To constrain models beyond LambdaCDM the development of the Euclid analysis pipeline requires simulations that capture the nonlinear phenomenology of such models. We present an overview of numerical methods and Nbody simulation codes developed to study the nonlinear regime of structure formation in alternative dark energy and modified gravity theories. We review a variety of numerical techniques and approximations employed in cosmological Nbody simulations to model the complex phenomenology of scenarios beyond LambdaCDM. This includes discussions on solving nonlinear field equations accounting for fifth forces and implementing screening mechanisms. Furthermore we conduct a code comparison exercise to assess the reliability and convergence of different simulation codes across a range of models. Our analysis demonstrates a high degree of agreement among the outputs of different simulation codes providing confidence in current numerical methods for modelling cosmic structure formation beyond LambdaCDM. We highlight recent advances made in simulating the nonlinear scales of structure formation which are essential for leveraging the full scientific potential of the forthcoming observational data from the Euclid mission.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Euclid preparation. Simulations and nonlinearities beyond LambdaCDM. 1. Numerical methods and validation
{'ABSTRACT': 'To constrain models beyond Λ CDM, the development of the Euclid analysis pipeline requires simulations that capture the nonlinear phenomenology of such models. We present an overview of numerical methods and N -body simulation codes developed to study the nonlinear regime of structure formation in alternative dark energy and modified gravity theories. We review a variety of numerical techniques and approximations employed in cosmological N -body simulations to model the complex phenomenology of scenarios beyond Λ CDM. This includes discussions on solving nonlinear field equations, accounting for fifth forces, and implementing screening mechanisms. Furthermore, we conduct a code comparison exercise to assess the reliability and convergence of di ff erent simulation codes across a range of models. Our analysis demonstrates a high degree of agreement among the outputs of di ff erent simulation codes, providing confidence in current numerical methods for modelling cosmic structure formation beyond Λ CDM. We highlight recent advances made in simulating the nonlinear scales of structure formation, which are essential for leveraging the full scientific potential of the forthcoming observational data from the Euclid mission.  \nKey words. cosmological N -body simulations - dark matter - dark energy', '1. Introduction': "Significant progress in cosmological observations is expected in the upcoming years, in particular from the Euclid survey (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2024; Euclid Collaboration: Scaramella et al. 2022; Nesseris et al. 2022; Martinelli et al. 2021; Euclid Collaboration: Castro et al. 2023), Vera Rubin Observatory's Legacy Survey of Space and Time (LSST, Ivezi'c et al. 2019), the Roman Space Telescope (Spergel et al. 2015) and the Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration: Aghamousa et al. 2016). These surveys will o ff er precision observations to high redshifts, allowing us to study the evolution of the Universe with unprecedented accuracy and potentially uncover the nature of dark matter and dark energy (DE). Gaining a deeper understanding of the nature of DE and addressing the long-standing question of whether the cosmological constant ( Λ ) is responsible for the late-time accelerated expansion of the Universe is indeed one of the primary goals of the Euclid survey (Amendola et al. 2018).  \nThe Euclid space telescope was launched on July 1, 2023, and is going to observe billions of galaxies out to redshift z ≈ 2, covering more than a third of the sky in optical and near-infrared wavelengths. Euclid will deliver precise measurements of the shapes and redshifts of galaxies (Euclid Collaboration: Bretonnière et al. 2022, 2023; Euclid Collaboration: Merlin et al. 2023; Euclid Collaboration: Desprez et al. 2020; Euclid Collaboration: Ilbert et al. 2021), from which we will measure weak gravitational lensing (Euclid Collaboration: Ajani et al. 2023) and galaxy clustering (Euclid Collaboration: Adam et al. 2019). These primary probes can be used to rigorously investigate different cosmological scenarios, in particular those related to DE that go beyond the Λ -Cold-Dark-Matter ( Λ CDM) concordance model.  \nAlthough the Λ CDM model is generally very successful in matching observations, the true identities of CDM and the cosmological constant Λ remain unknown. Additionally, some tensions have persisted in recent years, most notably the Hubble tension (see Di Valentino et al. 2021a, for a summary and references) where local measurements of the Hubble parameter today, H 0, appear to disagree with those inferred from high-redshift observations by around 5 σ . Further examples are the S 8 tension (see Di Valentino et al. 2021b, for a summary and references) and some anomalies found in measurements of the cosmic microwave background (Abdalla et al. 2022). The presence of these tensions may hint at a breakdown of the Λ CDM model and further motivates the exploration of alternative scenarios.  \nOver the past few years, cosmologists have explored di ff erent possibilities to account for the late-time accelerating expansion of the Universe (Tsujikawa 2010; Clifton et al. 2012; Joyce et al. 2016) either by introducing a new field, referred to as the DE field or by proposing a modified theory of gravity (MG). A wide range of MG or DE models is equivalent to adding a new light scalar degree of freedom to the theory of General Relativity (GR).  \nIn these theories, the scalar degree of freedom exhibits time evolution, sometimes accompanied by spatial fluctuations within the cosmic horizon. Even in the absence of such fluctuations, the background evolution may be di ff erent from Λ CDM, leading to modifications in structure formation. Significant spatial fluctuations in these models may arise due to various factors, including a low characteristic speed of sound in the theory (Gleyzes et al. 2014; Hassani et al. 2019), or as a result of the non-minimal coupling of the scalar field to matter or gravity (see Amendola 2004, for an example). MG and DE theories featuring a coupling of the scalar field to matter can further a ff ect perturbations at subhorizon scales by mediating a fifth force. If the coupling is universal and includes baryons, a screening mechanism is essential to evade the precise constraints of local experiments (Will 2014). Screening mechanisms are typically achieved through nonlinear phenomena in such theories. If, on the other hand, the coupling to matter is non-universal and is confined entirely to the dark sector, local experiments have no constraining power, and cosmological observations provide the main constraints.  \nGiven the diversity of possible DE or MG scenarios, a large information gain is expected from nonlinear scales in the cosmological large-scale structure. These scales must be studied using N -body simulations that capture the essential aspects of the DE or MG models under consideration. This usually means that at least one additional equation needs to be solved for the extra degree of freedom. In many cases, this leads to a di ffi cult nonlinear problem that could require special techniques or approximations that need to be developed. This makes N -body simulations for models of DE or MG a challenging task.  \nIn this paper, we first review the main features of the di ff erent classes of DE and MG models that have been proposed over the past years (see also Amendola et al. 2018, and Frusciante et al. in prep. for a more comprehensive and detailed overview). For each of them, we then discuss the numerical methods implemented within a selection of existing N -body codes (summarised in Table 1). Focusing on MG models with a universal coupling, we then compare the results of di ff erent N -body implementations for two well-studied theories, namely the Hu-Sawicki f ( R ) gravity (Hu & Sawicki 2007) and the 'normal branch' of the DvaliGabadadze-Porrati braneworld model (nDGP, Dvali et al. 2000;  \nSchmidt 2009a). We choose simulation parameters following the code comparison paper by Winther et al. (2015) [W15 hereafter], allowing us to validate a number of new codes against existing results.  \nThis article is part of a series that collectively explores simulations and nonlinearities beyond the Λ CDMmodel:  \n- 1. Numerical methods and validation (this work).\n- 2. Results from non-standard simulations (Rácz et al. in prep.).\n- 3. Cosmological constraints on non-standard cosmologies from simulated Euclid probes (D'Amico et al. in prep.).\n- 4. Constraints on f ( R ) models from the photometric primary probes (Koyama et al. in prep.).  \nFor further details, see our companion papers. The purpose of this first article in the series is to serve as a reference for models beyond Λ CDMand their existing implementations in various codes. This paper is structured as follows. In Sect. 2 we give a broad overview of di ff erent numerical approaches to treat the additional physics of models beyond Λ CDM. In Sect. 3 we discuss a number of di ff erent codes that implement those approaches and carry out a validation exercise, comparing several recently developed codes with the existing state of the art. We conclude in Sect. 4. In an Appendix, we discuss some performance considerations.", '2.1. Non-standard background evolution': "A wide range of models beyond the simplest cosmological constant scenario are based on an additional scalar degree of freedom - e.g. a classical scalar field ϕ - that evolves dynamically in the expanding Universe and whose background energy density ρϕ provides the source for the observed DE abundance. To induce cosmic acceleration and to match existing constraints on the background expansion history, the equation-of-state parameter w of such an additional field must be su ffi ciently negative at recent epochs but is poorly constrained at earlier times, which allows for models where w also evolves dynamically as long as it converges to values close to w ≈ -1 in the late Universe. For these models, the DE component modifies the background expansion history of the Universe, which is encoded by the general expression of the Hubble function,  \nH 2 ( a ) H 2 0 = Ω m a -3 + Ω r a -4 + Ω k a -2 + Ω DEe -3 R a 1 1 + w ( a ' ) a ' d a ' , (1)  \nwhere the equation-of-state parameter of DE can be obtained by solving the background field equations - including the evolution of the additional scalar degree of freedom ϕ - or can be parameterised. A common parameterisation suggested by Chevallier, Polarski & Linder (CPL, Chevallier & Polarski 2001; Linder 2003) is based on the desired evolution of w at low redshifts,  \nw ( a ) = w 0 + wa (1 -a ) . (2)  \nAlternatively, one can set the desired relative abundance of DE at late ( Ω DE = 1 -Ω m) and early ( Ω EDE) epochs as in the Early Dark Energy (EDE, Wetterich 2004) parameterisation,  \nw ( a ) = w 0 1 + b ln(1 / a ) , b = 3 w 0 ln 1 -Ω EDE Ω EDE + ln 1 -Ω m Ω m . (3)  \nThe modified expansion history expressed by Eq. (1) will indirectly a ff ect the evolution of matter density fluctuations and modify the formation process of collapsed structures by changing the Hubble friction term in the equation for linear matter perturbations, which in Newtonian gauge and in Fourier space for sub-horizon scales reads:  \n¨ δ m + 2 H ˙ δ m = 4 π G ( ρ m δ m + ρ DE δ DE) , (4)  \nwhere δ m and δ DE are the density contrasts of matter and DE perturbations, respectively, G is Newton's constant, and a dot represents a derivative with respect to cosmic time.  \nBesides the richer background dynamics that is endowed by an evolving field, whenever DE is promoted from a cosmological constant to a dynamical degree of freedom, the model also acquires an additional layer of complexity: the presence and evolution of DE fluctuations around the mean-field configuration. This corresponds to the situation where δ DE in Eq. (4) is nonnegligible, whereas in Λ CDM it would vanish identically at all scales. Like any other density perturbations, inhomogeneities in the DE would then contribute to the peculiar gravitational potential that governs the evolution of matter perturbations and thus the formation of cosmic structures as shown by Eq. (4).  \nHowever, in many of the simplest scalar-field scenarios, such perturbations are negligible at sub-horizon scales because the speed of sound c s of the scalar field is naturally close to the speed of light. Ignoring them for the purpose of numerical simulations, the only modification of N -body algorithms required to simulate these DE models is given by an appropriate calculation of the cosmic expansion rate. The most common approach amounts to tabulating the specific expansion rate of the universe for the model to be simulated according to Eq. (1) and replacing the standard analytical calculation of the Hubble function within the N -body algorithm with an interpolated value from the tabulated solution that is provided to the code as an input. This approach is implemented by most of the simulation codes employed within the Euclid Collaboration to perform cosmological simulations in homogeneous DE models beyond Λ CDM.", '2.2. Linearised DE perturbations': 'Although a wide range of DE models are characterised by negligible DE fluctuations as discussed above, some specific scenarios may not fulfil such a condition at all scales and / or at all times, either because they feature a lower value of the DE speed of sound, allowing DE perturbations to grow on scales above the associated Jeans length that then falls inside the cosmological horizon, or because additional interactions - besides gravity - can induce the growth of such perturbations. The former case corresponds to the class of clustering DE models, while the latter is known as coupled DE.', '2.2.1. Clustering DE': 'The clustering DE models are characterised by two timedependent variables: the speed of sound c s and the equation of state parameter w . In these theories, the DE component clusters on scales larger than the associated sound-horizon / Jeans scale, λ s = H / c s and DE perturbations decay quickly below the soundhorizon scale. For a su ffi ciently small speed of sound, we may even expect nonlinear DE structures to form. At a fundamental level, clustering DE models are analogous to the k-essence type of theories, so that the action reads (Armendariz-Picon et al.', '2001)': 'interacting species in the so-called Einstein frame,  \nS = Z d 4 x √ -g " c 4 R 16 π G + P ( X , ϕ ) + L m # , (5)  \nwhere P is a general function of the kinetic term X ≡ -1 2 ∇ µϕ ∇ µ ϕ and the scalar field ϕ , and L m is the matter Lagrangian. For a given P ( X , ϕ ), the speed of sound and the equation of state are given by (Armendariz-Picon et al. 2000)  \nw = P P -2 XP , X , c 2 s = P , X 2 XP , XX + P , X , (6)  \nwhere the subscript \' , X \' denotes the partial derivative with respect to X . We therefore need to specify the function P ( X , ϕ ) to derive the equations of motion for the k-essence scalar field. However, since there are many possible choices, we can instead employ the e ff ective field theory (EFT) approach to model the dynamics of k-essence DE. The EFT framework, although not a fundamental theory, o ff ers several advantages (Gleyzes et al. 2014; Cheung et al. 2008), such as being a description of a wide range of theories within some scales. The EFT is a perturbative approach based on the assumption that the scalar field perturbations remain small over the scales of interest. It is worth noting that the regime of nonlinear matter clustering is accessible to the EFT framework as long as the scalar field perturbations remain small. The k-essence theories or clustering DE models are implemented in several N -body and Einstein-Boltzmann codes. In CLASS (Lesgourgues 2011) and CAMB (Lewis et al. 2000), these theories are implemented using the fluid picture. In hi\\_class (Zumalacárregui et al. 2017), the EFT equations are implemented and can be controlled using the EFT parameter α K ≡ 3(1 + w ) c -2 s within the code. On the other hand, in k-evolution (Hassani et al. 2019, 2020), which is an N -body code based on gevolution (Adamek et al. 2016), nonlinear equations for clustering DE are implemented as an independent component, and the k-essence field for small c s can form nonlinear structures. In some N -body codes, for example in gevolution (Adamek et al. 2016), clustering DE is implemented through a linear solution from an Einstein-Boltzmann solver. This is a good assumption for large speeds of sound, but for small ones, this method does not allow for the response of DE to the nonlinear matter structures.', '2.2.2. Coupled quintessence': 'Moving to the case of coupled DE models, the interaction can be formulated at a fundamental level by introducing a direct coupling between the scalar field and the spatial curvature R in the so-called Jordan frame (see e.g. Pettorino & Baccigalupi 2008), so that the action reads  \nS = Z d 4 x √ -g " c 4 f ( ϕ, R ) 16 π G -1 2 Z ( ϕ ) ∇ µ ϕ ∇ µϕ -V ( ϕ ) + L m # , (7)  \nwhere f ( ϕ, R ) is a function that couples the scalar field to the curvature, Z ( ϕ ) is a function that allows for non-standard kinetic terms, V ( ϕ ) is the scalar field self-interaction potential, and the matter Lagrangian contains at least one cold species characterised by some rest mass m 0.  \nAlternatively, the interaction can be formulated by including source terms in the covariant conservation equations of the  \nArticle number, page 4 of 20  \n∇ µ T (c) µν = -β (c)( ϕ ) M Pl T (c) ∇ νϕ , (8)  \n∇ µ T (b) µν = -β (b)( ϕ ) M Pl T (b) ∇ νϕ , (9)  \n∇ µ T ( ϕ ) µν = 1 M Pl h β (c)( ϕ ) T (c) + β (b)( ϕ ) T (b) i ∇ νϕ , (10)  \nwhere T (Y) µν is the stress-energy tensor of a given species Y, T (Y) is its trace, β (Y)( ϕ ) is the coupling function of species Y, the labels c , b , ϕ refer to the dark matter, baryon, and scalar field species, respectively, and M Pl ≡ ( ℏ c ) 1 / 2 (8 π G ) -1 / 2 is the reduced Planck mass.  \nIn the case of a universal coupling (that is, if β (b) = β (c)), the two approaches can be related to one another through a Weyl transformation of the metric (see again Pettorino & Baccigalupi 2008), and are therefore equivalent. On the other hand, the possibility to leave the baryonic component of the Universe only minimally coupled evades Solar System constraints (see e.g. Will 2014) on the deviations from standard gravity thereby avoiding the need for screening mechanisms. This is the case of Coupled Quintessence models (Wetterich 1995; Amendola 2000), where the direct coupling between the scalar field and massive (nonbaryonic) particles can support stable perturbations of the DE field at sub-horizon scales (see e.g. Amendola 2004). In general, such perturbations may even become nonlinear in the presence of a su ffi ciently strong coupling (as in the case of Growing Neutrino Quintessence models, see e.g. Amendola et al. 2008; Mota et al. 2008; Baldi 2012b; Ayaita et al. 2016). Nonetheless, a large class of widely studied coupled DE models is known to feature scalar perturbations of the order of the standard Newtonian gravitational potential ( δϕ ∼ Φ N, see again Amendola 2004, for an extended derivation), thereby remaining in the linear regime at all times and scales of cosmological interest. This allows us to linearise the corresponding field equations and derive modified equations of motions for massive particles, including the contribution of the additional force arising from the direct coupling with the scalar field (see e.g. Baldi 2011).  \nWhile in the former case the interaction will be universal (i.e. involving all matter species with the same strength), which goes under the name of Extended Quintessence, the latter approach allows for non-universal couplings that may selectively involve individual species, for example by separately choosing the coupling functions for baryons and dark matter.  \nIn fact, a general feature of coupled DE models is the existence of a \'fifth force\' mediated by the scalar field. The new force can be expressed as an additional acceleration experienced by a massive coupled particle, which in comoving coordinates will be given by  \na Y , 5 th = -β (Y)( ϕ ) ∇ δϕ, (11)  \nwhere Y identifies a coupled matter species, and δϕ is the scalar field fluctuation. This extra acceleration term is added to the standard Newtonian acceleration acting on all massive particles,  \na N = -H v -∇ Φ N , (12)  \nwhere v is the peculiar particle velocity in comoving coordinates, and Φ N is the peculiar Newtonian potential obeying the standard Poisson equation  \n∇ 2 Φ N = 4 π Ga 2 X Y ρ Y δ Y , (13)  \nwhere the sum runs over all clustering species in the Universe.  \nTherefore, solving for the dynamical evolution of massive coupled particles requires solving for the scalar field perturbation δϕ entering in Eq. (11), which in the most general case follows a nonlinear elliptic equation,  \n∇ 2 δϕ = F ( δϕ ) + X Y 8 π Ga 2 β (Y)( ϕ ) δ Y , (14)  \nwith F a function of the scalar field fluctuation δϕ , and where the sum runs over all the coupled matter species with their respective couplings β (Y)( ϕ ).  \nFor the particular case of a coupled DE model with a nonuniversal interaction (Damour et al. 1990) involving only dark matter and leaving baryons uncoupled (i.e. β (b) = 0) the function F ( δϕ ) in Eq. (14) is negligible compared to the term associated with matter density perturbations (see Amendola 2004, for a derivation) and can be safely discarded. As a result, the scalarfield equation reduces to  \n∇ 2 δϕ ≈ 8 π Ga 2 β (c)( ϕ ) δ c = 2 β (c)( ϕ ) ∇ 2 Φ c , (15)  \nwhere Φ c is the Newtonian potential generated by the distribution of the coupled dark matter particles, that is  \n∇ 2 Φ c = 4 π Ga 2 ρ c δ c . (16)  \nTherefore, the solution for the scalar field perturbations will be directly proportional to the potential Φ c according to the relation  \nδϕ ≈ 2 β (c) Φ c . (17)  \nThe acceleration equation for a coupled particle can then be rewritten as  \na c = -H v c -∇ Φ b -∇ GLYPH<16> 1 + 2 β 2 (c) ( ϕ ) GLYPH<17> Φ c , (18)  \nassuming here for simplicity that other clustering species (such as massive neutrinos) give a negligible contribution to the total Newtonian potential such that Φ N = Φ c + Φ b. This modified acceleration equation introduces a further modification to be implemented in N -body simulation codes for Coupled Quintessence cosmologies besides the specific expansion history of each particular model. This often requires substantial modifications in the gravity solvers of conventional N -body codes, as the algorithms need to evolve coupled and uncoupled massive particles (typically dark matter and baryons, respectively) with di ff erent equations and should therefore treat these components separately. Even under the approximation of a purely collisionless treatment (i.e., ignoring the hydrodynamical and astrophysical processes that a ff ect standard baryonic matter leading to the formation of stars and galaxies) that is often employed for largevolume simulations targeted at galaxy surveys such as Euclid , both coupled and uncoupled matter species must be included in the simulation to provide a consistent representation of the dynamics at all scales: as baryons and dark matter evolve di ff erently, assuming that all matter is dark would lead to an overestimation of the e ff ects of the coupling and, thus, biased results. This approach is implemented in the C-Gadget code (Baldi et al. 2010; Baldi 2012a) which is employed for Coupled DE simulations performed within the Euclid Collaboration.  \nDistinguishing between coupled and uncoupled particle types in simulations of Coupled Quintessence is also crucial for proper treatment of two other e ff ects that characterise these cosmological models beyond the fifth force described by Eq. (11). The first is the mass variation of coupled particles due to the  \nexchange of rest-frame energy with the DE scalar field, which arises as a direct consequence of the modified continuity equations (8) and of the assumption of particle number conservation. More specifically, the mass of coupled particles evolves as a result of the evolution of the background scalar field according to  \nm Y( a ) = m 0 exp -Z ϕ ( a ) ϕ 0 β (Y)( ϕ ) d ϕ M Pl ! . (19)  \nSuch a mass variation, which involves only particle species with a non-vanishing coupling to the scalar field, must be taken into account in N -body algorithms by changing the mass of individual simulation particles at every time step. This is normally done by tabulating the mass as a function of scale factor a by numerically integrating Eq. (19) along with the background dynamics of the scalar field ϕ , and interpolating from that table as the simulation progresses.  \nThe second e ff ect is an additional force (on top of the fifth force) acting on coupled particles as a consequence of momentum conservation due to the particles\' mass variation described by Eq. (19), which takes the form of a velocity-dependent extra acceleration behaving either as a friction or as a drag, depending on the relative signs of the coupling function β ( ϕ ) and of the scalar field velocity ˙ ϕ (see e.g. Baldi et al. 2010),  \na Y , v = -β (Y)( ϕ ) M Pl ˙ ϕ v . (20)  \nSuch a velocity-dependent acceleration is responsible for a very rich phenomenology characterising Coupled Quintessence models, especially on highly nonlinear scales (see e.g. Baldi et al. 2010; Baldi 2012a; Li & Barrow 2011; Baldi 2023), and must be included in N -body simulations as well for a fully consistent treatment of these scenarios. This is done by adding the extra acceleration described in Eq. (20) to the total acceleration (i.e. Newtonian plus fifth force) of all coupled particles in each time step,  \na Y = a N + a Y , 5 th + a Y , v . (21)  \nThe relevant quantities β ( ϕ ) and ˙ ϕ can again be interpolated from a table obtained by integrating the background dynamics of the system. This is the approach implemented in the C-Gadget code that has been used to run Coupled Quintessence simulations within the Euclid Collaboration.', '2.2.3. Momentum exchange and dark scattering': "A further example of interacting DE cosmologies characterised by scalar-field perturbations that always remain linear is given by models of pure momentum exchange (see e.g. Pourtsidou et al. 2013; Skordis et al. 2015) between the DE field and massive particles like dark matter or baryons. A limiting case is given by the Dark Scattering scenario (Simpson 2010) where the momentum transfer between the two components is modelled as the elastic scattering of massive particles moving through a homogeneous DE fluid with equation of state w . This results in an extra force acting on the moving massive particles which is proportional to their comoving velocity, similar to the velocity-dependent force described by Eq. (20) for Coupled Quintessence models. However, the origin of this force is completely di ff erent in this case, as it does not originate from the mass variation of particles but rather from the momentum transfer with the DE field. As a result, the scattering acceleration can be expressed as  \na s = -(1 + w ) 3 H 2 Ω DE 8 π G ξ v , (22)  \nwhere the parameter ξ is defined as  \nξ ≡ σ m , (23)  \nwith σ denoting the scattering cross section and m the typical mass of the scattering particle species.  \nAlthough Dark Scattering represents a limiting case of the more general class of pure momentum-exchange models between matter and DE (also known as 'Type 3' models in the classification of Skordis et al. 2015), for which further modifications to the standard particle dynamics are expected besides the drag force of Eq. (22), recent works (Palma & Candlish 2023) have shown that such additional modifications are generally subleading with respect to the drag force so that their e ff ect on structure formation can be neglected. This ensures that the current implementation of Dark Scattering within the simulations used in the Euclid Collaboration can be considered representative of the general class of momentum-exchange cosmologies.  \nThis type of interaction can be implemented in N -body algorithms (see e.g. Baldi & Simpson 2015, 2017) in a very similar way as the velocity-dependent acceleration in Coupled Quintessence scenarios, as the factors entering Eq. (22) are all either constants or background quantities that can be interpolated at every timestep from tabulated data. This is the approach implemented in the C-Gadget code that has been used to run the DAKAR and DAKAR2 simulations (Baldi & Simpson 2017, Rácz et al. in prep.).", '2.3. Nonlinear scalar field perturbations': "In models where a scalar field couples to matter universally, or at least to baryons in a relevant way, some mechanism to suppress the coupling is required to satisfy the stringent local tests of gravity. This is commonly referred to as 'screening'. Screening mechanisms are achieved by nonlinearity in the scalar-field equation coupled to matter. The equation determining the evolution of the scalar field is typically a wave equation of the form  \n□ ϕ = S ( ϕ, ∇ µϕ, ∇ µ ∇ νϕ, ρ m) . (24)  \nSeveral approximations are often used to solve these nonlinear equations. The most common one is the quasi-static approximation. The scalar field can be split into a background part, ¯ ϕ , and a perturbation, δϕ , as ϕ = ¯ ϕ + δϕ . The quasi-static approximation amounts to ignoring the time dependence of the scalar field perturbation, i.e. assuming ˙ ϕ ≃ ˙ ¯ ϕ . The partial differential equation (PDE) of the field perturbation, which in its original form may have been of the hyperbolic or parabolic type, is therefore cast into an elliptic form so that the scalar field solution at any given time depends solely on the matter configuration at that time. This is a good approximation whenever the speed of sound of the scalar field is small (Sawicki & Bellini 2015), which is the case for the MG models considered here. Non-quasistatic cosmological simulations have been conducted for several MG models using di ff erent techniques, such as the explicit leap-frog method and the implicit Newton-Gauss-Seidel method (Llinares & Mota 2013; Bose et al. 2015; Winther & Ferreira 2015b).  \nHere, □ ≡ ∇ µ ∇ µ represents the d'Alembertian operator and S is a nonlinear function that depends on the matter density, the scalar field, and its derivatives. Various methods have been developed to solve this nonlinear scalar field equation in N -body simulations, where the nonlinear density ρ m is modelled by collisionless particles (see W15, for more details).  \nThe scalar-field solution is required for the computation of the total gravitational potential Φ that acts on the matter particles,  \n∇ 2 Φ = 4 π Ga 2 GLYPH<16> δρ m + δρ e ff ( ϕ ) GLYPH<17> , (25)  \nwhere the e ff ective density depends on the scalar field. There are two common ways of solving for this total gravitational force. The first option is to solve first for ϕ and then use this solution to compute the source term in Eq. (25) and solve for Φ using a standard Poisson solver to get the total force ∇ Φ . The other option is to apply the total force ∇ Φ N + ∇ ϕ to the particles.  \nUnder the quasi-static approximation, the scalar-field equation assumes the same form as the usual Poisson equation,  \n∇ 2 ϕ = S ( ϕ, ∇ i ϕ, ∇ i ∇ j ϕ, δρ m) . (26)  \nThe main di ff erence is that the scalar-field equation is generally nonlinear. This nonlinear behaviour implies that conventional techniques, such as using Fourier analysis, cannot be used to solve the equation. Numerous approaches have been developed to address this challenge, and we refer the reader to W15 for details. For computational methods that aim to accurately solve a nonlinear equation on refined grids, the approach typically involves discretising the equation in a suitable way and employing an iterative algorithm, such as the Newton-Raphson method, to successively refine solutions based on an initial guess. To speed up convergence, many of these methods incorporate socalled 'multigrid' acceleration techniques which we quickly review here.", '2.3.1. Nonlinear multigrid algorithm': 'A generic way to solve nonlinear elliptic PDEs is to couple the multigrid algorithm to the Newton-Raphson method,  \nu new = u old -L ( u old ) ∂ L /∂ u old , (27)  \nwhere u is the discretised field, L is the di ff erential operator (which is a Laplacian for Newtonian gravity) and the superscripts refer to the new or old estimate of the solution in one Newton-Raphson iteration. The Newton-Raphson method produces linear equations for the correction terms, which are solved by the Full-Approximation-Storage Multigrid algorithm (Brandt 1977; Wesseling 2004; Guillet & Teyssier 2011). For a review of these methods applied to MG simulations, see e.g. Li (2018), Llinares (2018) and W15.  \nA simple sketch of the algorithm goes as follows. One starts with a guess for the solution on a grid; this could be anything from a constant value across the grid to using the solution from the previous timestep in the simulation. One then loops a few times over all cells in the grid, updating the solution using Eq. (27). This solution is then restricted to a grid with half the resolution, the solution is updated again, and this process is repeated recursively up to the coarsest grid (one with only 2 3 cells). The solution is then interpolated to the finer grid, updated once more, and this is done recursively until one reaches the finest grid we started with. One such cycle is called a V-cycle, and one repeats such V-cycles until convergence is achieved. The advantage of having this stack of coarser grids is that it helps to accelerate the convergence of the largest modes in the solution.', '2.3.2. Screening with nonlinearity in potentials': 'In models where screening is achieved by nonlinearity in a potential or coupling function, the equation for the scalar field becomes  \n∇ 2 ϕ = 4 π Ga 2 β ( ϕ ) δρ m + V ( ϕ ) , (28)  \nwhere β ( ϕ ) and V ( ϕ ) are nonlinear functions of ϕ . A typical example is f ( R ) gravity. In this class of models, the value of the scalar field changes by orders of magnitude. To enhance numerical stability, a common technique involves redefining the scalar field in terms of a new variable. The redefinition to choose depends on the specific model under consideration. It is typically chosen to prevent the occurrence of unphysical values of the scalar field during Newton-Raphson iterations. For example, for f ( R ) models the scalar field ϕ = f , R will be driven towards zero in high-density regions, but at the same time f , R cannot cross zero, as the potential becomes singular in this scenario. To avoid this issue, a commonly used field redefinition is u ≡ ln[ f , R / ¯ f , R ( a )] (Oyaizu 2008). However, this transformation introduces additional nonlinearity and in some models, such as Hu-Sawicki f ( R ) models, this transformation is not necessary and might even lead to considerable performance losses in a simulation.  \nIn some cases, this can be avoided. For example, Bose et al. (2017) noticed that for Hu-Sawicki f ( R ) gravity with n = 1 (Hu & Sawicki 2007), when making the change of variable u = p -f , R , the field equation could be recast as a depressed cubic equation,  \nu 3 + pu + q = 0 , (29)  \nwhich possesses analytical solutions (Ruan et al. 2022). Although the Gauss-Seidel smoothing procedure is still needed (because p depends on the values of the field u in neighbouring cells), this removes the Newton-Raphson part and expensive exponential / logarithmic operations from the method of Oyaizu (2008), therefore leading to significant performance gains. Ruan et al. (2022) also generalised this improved relaxation approach to the cases of n = 0 (which strictly speaking is not a variation of the Hu-Sawicki model) and n = 2.  \nFor other models like the symmetron, which has a Higgslike potential, the scalar field is free to cross zero, and no field redefinition is needed (apart from a simple rescaling).', '2.3.3. Screening with nonlinearity in kinetic terms': 'In another class of models, nonlinearity emerges within the kinetic term. For example, in models with the Vainshtein screening mechanism (Vainshtein 1972), the equation exhibits nonlinearity in the second derivatives of the scalar field,  \n∇ 2 ϕ = 4 π Ga 2 β ( a ) δρ m + g ( ∇ i ϕ, ∇ i ∇ j ϕ ) , (30)  \nwhere β ( a ) is a time-dependent coupling function. The simplest example here is the DGP model which was first simulated by Schmidt (2009b). In such cases, the operator-splitting trick (Chan & Scoccimarro 2009) can be employed. This approach can simplify the equations, avoiding potential issues associated with imaginary square roots, and improving code performance. This trick is particularly useful for the DGP braneworld models and other Vainshtein screening models, such as Cubic and Quartic Galileons (see Sect. 4.2.2 of W15, for more information).', '2.3.4. Approximate treatments of screening': "Some models allow linearisation of the nonlinear equation using some approximation. One approach (Khoury & Wyman 2009; Winther & Ferreira 2015a; see also the appendix of Schmidt 2009a) is to introduce the screening factor for the matter density perturbation  \n∇ 2 ϕ = c 2 m 2 a 2 ℏ 2 ϕ + 4 π Ga 2 δρ m ϵ screen( Φ N , |∇ Φ N | , ∇ 2 Φ N) , (31)  \nwhere the screening function depends on the Newtonian potential Φ N. This type of parameterised modified gravity is referred to as 'type 1' in Table 1. One specific method, developed in Brando et al. (2023), starts from linearising Klein-Gordon's equation. In this formalism, one solves the Poisson equation in Fourier space,  \n-k 2 Φ = 4 π G e ff ( a , k ) a 2 δρ m , , (32)  \nwhere the function G e ff ( a , k ) approximates the e ff ective Newton's constant introduced by the screening e ff ect of the scalar field on small scales. This function is given by  \nG e ff ( a , k ) = G + ∆ G e ff ( a , k ) = G + GLYPH<16> G lin e ff ( a ) -G GLYPH<17> M ( a , k ) , (33)  \nwith G lin e ff ( a ) being the asymptotic linear e ff ective Newton's constant that depends only on time, and M ( a , k ) is a function that approximately captures the nonlinear corrections introduced by the scalar field on small scales. This function allows Eq. (33) to transition from G e ff ( a , k ) → G lin e ff ( a ) on large scales to G e ff ( a , k ) → G on small scales. This type of parameterised modified gravity is referred to as 'type 2' in Table 1. A procedure to fix M ( a , k ) is described in Brando et al. (2023), which has the advantage of avoiding additional parameters to tune the screening e ffi ciency.  \nOne can also choose to parameterise the nonlinear contribution using an e ff ective Newton's constant at both small and large scales. If the modifications of gravity are encoded in a scaledependent function ∆ G e ff ( a , k ) as in Eq. (33), then we can propose a similar equation in real space,  \n˜ G e ff ( a , r ) = G + ∆ ˜ G e ff ( a , r ) , (34)  \nwhere ∆ ˜ G e ff ( a , r ) is the Fourier transform of ∆ G e ff ( a , k ). In practice, an additional approximation is made, namely ∆ ˜ G e ff ( a , r ) ≈ ∆ G e ff ( a , k → 1 / r ).  \nThis approach allows the encoding of nonlinear contributions over the whole range of scales modelled by N -body algorithms through real-space equations, for instance, the Tree Particle-Mesh (TreePM) method implemented in codes like Gadget4 . Provided the parameterisation is e ff ective, this is expected to increase the accuracy of the estimation of the nonlinear e ff ects.  \nSeveral parameterisations have been proposed for this kind of approach, either with additional tuning parameters, such as in Lombriser (2016), or based on local small-scale environmental properties to avoid the need for any extra parameters, such as in Winther & Ferreira (2015a).", '3. Additional code validation': 'A code comparison for simulations that solve nonlinear scalar field perturbations is presented in W15. Since then, various new codes have been developed. In this section, we show a comparison of the predictions for the power spectrum and the halo mass  \nTable 1: Summary table of the N -body codes implementing various extensions to the standard Λ CDM cosmology that have been used to produce simulations employed in Euclid pre-launch analysis, validation, and forecasting.  \nfunction using the simulations of W15 as a reference and starting from the same initial conditions. These were generated using second-order Lagrangian perturbation theory in a Λ CDM cosmology with Ω m = 0 . 269, ΩΛ = 0 . 731, h = 0 . 704, n s = 0 . 966 and σ 8 = 0 . 801. The simulations have N p = 512 3 particles of mass M p ≃ 8 . 756 × 10 9 h -1 M ⊙ in a box of size B = 250 h -1 Mpc and start at redshift z = 49. As in W15, we compare simulations for f ( R ) and nDGP models. In these models (Schmidt 2009a), the background expansion history is closely approximated by that of Λ CDM. Furthermore, the e ff ect of modified gravity can be ignored at z = 49, thus it is justified to use the initial conditions of Λ CDM. The measurements of the power spectrum and mass function are performed by the pipeline developed in the second article of this series, Rácz et al. (in prep.). Based on two models only, our comparison does not encompass the full diversity of numerical methods discussed in the previous section. In many cases some validation of the various implementations can be found in the corresponding references.', '3.1. Summary of codes used in the validation': 'Table 1 shows an overview of the simulation codes considered in this section and provides a quick reference of their capabilities  \nand limitations. For each of them, a short summary is presented here. In the Appendix we comment on the trade-o ff between accuracy and computational cost of the implementations.', '3.1.1. MG-Arepo': "First presented in Arnold et al. (2019), this code is based on the moving-mesh N -body and hydrodynamical simulation code Arepo (Springel 2010; Weinberger et al. 2020), which uses a TreePM algorithm to calculate gravitational forces. The additional modified gravity force (fifth force) is calculated with a relaxation solver (Bose et al. 2017) that is accelerated by the multigrid method and uses adaptive mesh refinement (AMR). It currently also supports simulations for the nDGP model (Hernández-Aguayo et al. 2021), as well as massive neutrinos implemented using the δ f method (Elbers et al. 2021).  \nTo solve the modified gravity equations, the density field is projected onto the AMR grid, constructed in such a way that each cell on the highest refinement level contains at most one particle (except if a pre-set maximum refinement level is reached; the cell size at this level is of the order of the smoothing length of the standard gravity solver). Once the field equation is solved to  \nobtain the scalar field configuration, the modified gravity force can be computed from its gradient using finite di ff erencing.  \nSince MG-Arepo computes the forces on all particles simultaneously, and the modified gravity field equations are generally highly nonlinear (with a poor convergence rate of the relaxation algorithm), this is computationally expensive compared to Arepo 's Newtonian gravity solver. However, the maximum acceleration of the modified gravity force is smaller than that of Newtonian gravity, mainly because the latter occurs in regions with high density where screening occurs. This allows the modified gravity solver to run using larger time steps (Arnold et al. 2019), resulting in significantly reduced computational cost. Together with Arepo 's e ffi cient MPI parallelisation and lean memory footprint, these have made it possible to run the large number of f ( R ) simulations used in various recent works, such as the FORGE -BRIDGE (Arnold et al. 2022; Harnois-Déraps et al. 2023; Ruan et al. 2024) simulation suite of 200 f ( R ) and nDGP models. This has allowed accurate emulators of various physical quantities or observables to be constructed.  \nAnother highlight of MG-Arepo is its capability to run realistic galaxy formation simulations in a cosmological box (Arnold et al. 2019; Hernández-Aguayo et al. 2021), thanks to the use of the Illustris-TNG subgrid physics model (Weinberger et al. 2017; Pillepich et al. 2018). More recently, it has been used for larger-box hydrodynamical simulations with a realistic recalibrated Illustris-TNG model (Mitchell et al. 2022), enabling the study of galaxy clusters in modified gravity.  \nThe MG-Arepo simulations used in this paper were run using a residual criterion of ϵ = 10 -2 and a maximum refinement level ( MaxAMRLevel ) of 10 for the nDGP simulations and 18 for the f ( R ) models with a gravitational softening of 0 . 01 h -1 Mpc.", '3.1.2. MG-Gadget': 'MG-Gadget (Puchwein et al. 2013) is a modified version of the TreePM N -body code Gadget-3 (which in turn is based on the public code Gadget-2 , see Springel 2005) implementing an AMR solver for the scalar degree of freedom f , R characterising the widely-studied Hu-Sawicki f ( R ) gravity model. In MG-Gadget , the same tree structure that is employed to solve for standard Newtonian gravity is also used as an adaptive grid to solve for the scalar field configuration through an iterative Newton-Gauss-Seidel (NGS) relaxation scheme (see Sect. 2.3) with the Full-Approximation-Storage Multigrid method (see Sect. 2.3.1) and with the field redefinition u ≡ ln[ f , R / ¯ f , R ( a )] (see Sect. 2.3.2). MG-Gadget also allows MG simulations to be run with massive neutrinos (see e.g. Baldi et al. 2014; Giocoli et al. 2018) using the neutrino particle method (see Adamek et al. 2023, for a review on numerical methods for massive neutrino simulations).  \nFor the simulations presented here, the relative tree opening criterion was used with an acceleration relative error threshold of 0.0025, and a uniform grid with 512 3 cells was employed to compute long-range Newtonian forces. Concerning the MG field solver, a residual tolerance of ϵ = 10 -2 was set for the V-cycle iteration and a maximum refinement level of 18 was used for the AMR grid, corresponding to a spatial resolution of 1 h -1 kpc at the finest grid level, compared to a gravitational softening of 18 h -1 kpc, following the setup adopted in W15.', '3.1.3. ECOSMOG': 'ECOSMOG (Li et al. 2012) is a generic modified gravity simulation code based on the publicly-available N -body and hydrodynamical simulation code RAMSES (Teyssier 2002). Originally developed for f ( R ) gravity, this code takes advantage of the adaptive mesh refinement of RAMSES to achieve the high resolution needed to solve the scalar field and hence the fifth force in highdensity regions. The nonlinear f ( R ) field equation is solved with the standard Gauss-Seidel approach as first applied by Oyaizu (2008), but it was later replaced by the more e ffi cient algorithm of Bose et al. (2017). The code has since been extended for simulations for the generalised chameleon (Brax et al. 2013), symmetron and dilaton (Brax et al. 2012), nDGP (Li et al. 2013b), cubic Galileon (Barreira et al. 2013, 2015), quartic Galileon (Li et al. 2013a), vector Galileon (Becker et al. 2020) and nonlocal gravity (Barreira et al. 2014).  \nFor the ECOSMOG simulations used in this paper, we have used a domain grid (the uniform mesh that covers the whole simulation domain) with 2 9 = 512 cells per dimension, and the cells are hierarchically refined if they contain 8 or more e ff ective 1 particles. The highest refined levels have e ff ectively 2 16 cells, leading to a force resolution of about 0 . 0075 h -1 Mpc.', '3.1.4. ISIS': 'The ISIS code (Llinares et al. 2014), like the ECOSMOG code above, is based on RAMSES (Teyssier 2002). It contains a scalar field solver that can be used to simulate generic MG models with nonlinear equation of motion and has been used to simulate models such as f ( R ) gravity, the symmetron model, nDGP and disformal coupled models (Gronke et al. 2014; Winther & Ferreira 2015a; Winther et al. 2015; Hagala et al. 2016; Llinares et al. 2020). It also allows for hydrodynamical simulations that have been used to study the interplay between baryonic physics and modified gravity (Hammami et al. 2015). Furthermore, the code has the capability to go beyond the quasistatic limit and study the full time dependence of the scalar field (Llinares & Mota 2013; Llinares & Mota 2014; Hagala et al. 2017). The scalar field solver used in the code is a Gauss-Seidel relaxation method with multigrid acceleration, very similar to the one in ECOSMOG described above.  \nFor the ISIS simulations presented in this paper, we have used a domain grid (the uniform mesh that covers the whole simulation domain) with 2 9 = 512 cells per dimension, and the cells are hierarchically refined if they contain 8 or more e ff ective particles.', '3.1.5. PySCo': 'PySCo 2 is a particle-mesh (PM) code written in Python and accelerated with Numba which currently supports Newtonian and f ( R ) gravity (parameterised as in Hu & Sawicki 2007, with n = 1 or n = 2). While multiple flavours of solvers based on Fast Fourier Transforms (FFT) are available, for the present paper, we use a multigrid solver to propose something di ff erent from other codes in this comparison project (other codes that also use multigrid are AMR-based). PySCo uses a triangular-shaped  \ncloud mass assignment scheme and solves the linear Poisson equation using multigrid V-cycles with a tolerance threshold of the residual of 10 -3 , and two F-cycles (such cycles go through the mesh levels more often than V-cycles, resulting in a higher convergence rate of the residuals at the cost of increased runtime, for details on multigrid cycles see also Ruan et al. 2022) to solve the additional field in f ( R ) gravity with the nonlinear multigrid method described in Sect. 2.3.1 and Eq. (29). Our convergence threshold is very conservative since we do not intend to conduct a convergence study in this paper (a less conservative threshold could still give reasonable results at much lower computational cost). Furthermore, to resolve the small scales we use a coarse grid with 2048 3 cells, resulting in roughly 500 time steps to complete the simulations.', '3.1.6. MG-COLA': 'MG-COLA simulations were performed using the FourierMultigrid Library FML . 3 These simulations are based on the COmoving Lagrangian Acceleration (COLA) method (Tassev et al. 2013), which combines Lagrangian perturbation theory with the PM method to reduce the number of time steps that are required to recover clustering on large scales. The MG-COLA N -body solver in the FML library contains implementations of various DE and MG models like the DGP model, the symmetron model, f ( R ) gravity and the Jordan-Brans-Dicke model (see Winther et al. 2017, for more details). For the PM part, we used N 1 / 3 mesh = 5 N 1 / 3 p , i.e. a mesh discretisation five times smaller than the mean particle separation, and the total of 150 time-steps linearly spaced along the scale factor to achieve a good agreement of the mass function with AMR simulations. We also ran lowresolution simulations with N 1 / 3 mesh = 3 N 1 / 3 p with 100 time steps, and checked that the nonlinear enhancement of the power spectrum from these low-resolution simulations agrees well with the one from the high-resolution ones. Screening is included using approximate treatments described in Sect. 2.3.4. For f ( R ) gravity models, the FML code has one parameter to tune the strength of chameleon screening called screening efficiency . For the model with ¯ f , R = 10 -6 , we used the default value screening efficiency=1 while we used screening efficiency=2 for ¯ f , R = 10 -5 . For nDGP, we used Gaussian smoothing with a smoothing radius of 1 h -1 Mpc to compute the density field for screening and did not use an option to enforce the linear solution at small wavenumbers k . We also ran simulations based on the screening approximation using G e ff ( a , k ) given by Eq. (33) for nDGP, and in these simulations, we have used the same COLA settings as in the other screening approximation implementation. The biggest advantage of this screening implementation is that it does not require any additional tuning parameter related to screening, that is, the function G e ff ( a , k ) is completely defined by the theoretical model one wants to simulate.', '3.1.7. PANDA': "PANDA is an extension of the TreePM code Gadget4 . It introduces modifications at large scales in the PM part with an effective Newton's constant according to Eq. (33), while the force at small scales is modified in the tree part with Eq. (34). In this respect, PANDA implements an approximate solver for the extra force induced by di ff erent possible MG theories. On the other hand, di ff erently from other approximate methods, the dynamics of matter particles under the e ff ect of the (dominant) standard  \ngravity force is treated with the full TreePM solver of Gadget4 of which it retains the accuracy in modeling the nonlinear density field. For nDGP, the functional form of ˜ G e ff ( a , r ) is described in Lombriser (2016), while for f ( R ) it is introduced as in Winther & Ferreira (2015a). In the latter case, no additional parameters are needed to describe the screening. For the first nDGP case instead, in addition to the theoretically defined parameters N 0 = 1, B = 1 / [3 β ( a )], b = 2 and a = 3, we consider the screening scale k th = 0 . 4 h Mpc -1 defined by the parameterisation  \nr V GLYPH<27> 2 3 3 Ω m H 2 r 2 c β 2 ! 1 / 3 r th . (35)  \nwhere r V is the Vainshtein radius of the model for a spherically symmetric density perturbation (see Lombriser 2016, for more details).", '3.2. Comparison of the power spectrum': "Before discussing the e ff ect of extra degrees of freedom, we compare the di ff erent codes within a Λ CDM cosmology. For reference and only for the case of Λ CDM, we include results from the code PKDGRAV3 that was also used for the Flagship simulations in Euclid (Potter et al. 2017; Euclid Collaboration: Castander et al. 2024). Figure 1 shows the matter power spectra at two di ff erent redshifts, z = 1 (left panels) and z = 0 . 667 (right panels). The lower panels show the relative di ff erence to the simulation carried out with MG-Arepo that we use as a reference throughout. Since ISIS and ECOSMOG have both been developed from the RAMSES code, the agreement between these two codes is better than 1%. These AMR codes tend to underestimate the power on small scales when compared to tree-based codes, mainly due to the mesh refinement criterion used in our simulations. Fine-tuning some precision parameters, we expect that a better agreement can be achieved, cf. Schneider et al. (2016). As we can see, the codes that perform closely to our benchmark at small scales are the ones that also share the same tree-PM gravity solver, PANDA and MG-Gadget . PKDGRAV3 , which uses a tree structure and the fast multipole method to compute gravitational forces, is in excellent agreement with MG-Gadget up to k ∼ 1 h Mpc -1 , yielding slightly more power at higher k .  \nFor modified gravity models, we compare the ratio between the power spectrum computed for that model and the one found in the Λ CDM reference cosmology. Taking such ratios largely cancels out the di ff erences seen in Λ CDM between the codes, making the di ff erences due to the treatments of the extra scalar field more apparent. Following W15, we consider two MG models that exhibit two di ff erent types of growth dependence and  \nAlso shown in the figure are the results from COLA , a pure PM code that uses a fixed grid to solve the Poisson equation. The fixed resolution of the PM grid explains why COLA consistently su ff ers from low force resolution at k ≳ 1 h Mpc -1 . All our simulations use the exact same initial particle data, such that our comparisons should not be contaminated by cosmic variance. However, for COLA simulations, we also reconstructed the initial density field from the initial particle distribution to estimate the displacement fields required for the 2LPT calculations in those simulations. While PySCo incorporates a full N -body solver, its small-scale accuracy is constrained by resolution limitations that arise from the absence of AMR. Moreover, within the scope of this code comparison project, PySCo employs a multigrid algorithm, resulting in reduced small-scale clustering compared to FFT (on a regular grid), thereby explaining the observed deficiency in power at wavenumbers k ≳ 1 h Mpc -1 .  \n∆  \nFig. 1: Matter power spectra from simulations carried out with di ff erent codes (di ff erent line styles) for a Λ CDM cosmology. The left panels show the spectra at redshift z = 1 whereas the right panels show the spectra at redshift z = 0 . 667. The bottom panels show the relative di ff erence with respect to the simulation carried out with MG-Arepo .  \n<!-- image -->  \nscreening mechanisms. The first one is the so-called Hu-Sawicki f ( R ) which exhibits a scale-dependent growth factor at linear order, and a screening mechanism realised by nonlinearities in the potential, see Sect. 2.3.2. The strength of the modification of gravity is characterised by the background value of the scalar degree of freedom, and we study the cases ¯ f , R = 10 -5 and ¯ f , R = 10 -6 , labelled 'F5' and 'F6', respectively, the latter being closer to GR. The second one is nDGP, a braneworld theory defined in a five-dimensional spacetime. The growth factor in this theory is scale independent, and screening is realised through nonlinearities in the kinetic terms, as discussed in Sect. 2.3.3. Here, the strength of the modification of gravity is characterised by the value of the cross-over scale, and we study the cases r c = 1 . 2 H -1 0 and r c = 5 . 6 H -1 0 , labelled 'N1.2' and 'N5.6', respectively, the latter being closer to GR.  \nFigure 2 shows the relative change (with respect to Λ CDM) of the matter power spectrum due to modifications of gravity, sometimes called 'boost', for the two f ( R ) scenarios at redshifts z = 1 (left panels) and z = 0 . 667 (right panels). Since f ( R ) already exhibits a scale-dependent growth factor at linear order, we can see that on scales of k ≈ 0 . 1 h Mpc -1 this e ff ect is already present, and gets enhanced at large wavenumbers due to nonlinearities of the density field. As we can see the codes MG-Arepo , MG-Gadget , ECOSMOG , ISIS , and PySCo all roughly agree to better than one percent down to scales of k ≈ 10 h Mpc -1 , as expected. These minor discrepancies can be caused by the refinement criterion used in the latter code or by slight variations in the redshift of the particle snapshot output. The present results are largely in agreement with a similar analysis performed in W15.  \nIn the same plots, we can also see how approximation schemes to introduce screening from nonlinearities in the potential perform in contrast to the exact solutions. The results that use these methods are showcased by the examples of COLA and PANDA , where each uses a di ff erent scheme to approximate the e ff ects of the dynamics of the scalar field in dense environments.  \nAs expected, the two codes do not exhibit the same level of agreement down to scales of k ≈ 10 h Mpc -1 as their counterparts using full solvers. However, in both cases, we can see that the deviations from the MG-Arepo reference results are limited to 2% even in the most extreme regime of departure from GR. Through closer inspection, we can see small di ff erences in the agreement between the approximate schemes and full solvers at di ff erent values of k . These deviations are caused by di ff erent implementations of the approximate MG solvers in the two codes. In fact, while COLA is an approximate method that uses a PM algorithm to solve the Poisson equation on a fixed grid with the use of the screening approximation to linearise the Poisson equation, PANDA is a new implementation that exploits the TreePM structure of the baseline Gadget-4 code to solve for the small-scale particle dynamics by incorporating the MG e ff ects (including the screening) through a scale-dependent Newton's constant in both real and Fourier space.  \nFigure 3 shows results for the matter power spectra in the two nDGP scenarios we consider, characterised by the two values of the cross-over scale, r c = 1 . 2 H -1 0 (N1.2) and r c = 5 . 6 H -1 0 (N5.6). As before, we compare simulation results at two different redshifts, z = 1 (left panels) and z = 0 . 667 (right panels). Since the linear growth function has a scale-independent enhancement compared to Λ CDM, we see a constant amplification of power at small values of k . The Vainshtein mechanism suppresses the deviation from Λ CDM on small scales, leading to a diminishing boost at large k . The agreement between different codes at k < 1 h Mpc -1 is better than 1% for all codes, including the approximate simulations with COLA and PANDA . At k > 1 h Mpc -1 , the AMR-based codes, ECOSMOG and ISIS , show larger deviations at the level of 2%, where the suppression of the deviation from Λ CDM is slightly underestimated compared with MG-Arepo . This could be caused by the di ff erence in the Λ CDMpower spectrum shown in Fig. 1, even though this di ff erence in the baseline code is largely cancelled out in the boost. On the other hand, despite their approximations in the MG solvers,  \nFig. 2: Amplification factor of the matter power spectra from simulations of two f ( R ) scenarios, ¯ f , R = 10 -5 (F5, top row ) and ¯ f , R = 10 -6 (F6, bottom row ), relative to the reference Λ CDM cosmology. The amplification factor BP = P / P Λ CDM is measured for di ff erent codes (di ff erent line styles) at two di ff erent redshifts, z = 1 ( left panels ) and z = 0 . 667 ( right panels ). The bottom panel of each plot shows the relative agreement of the individual measurements, using MG-Arepo as a common reference.  \n<!-- image -->  \nCOLA and PANDA agree with MG-Arepo at the level of 1% even on these scales.", '3.3. Comparison of the halo mass function': 'To gain further insight into the nonlinear dynamics of the di ff erent models and their respective implementations we also compare the cumulative halo mass functions measured in our simulations. For this purpose, halo catalogues are obtained with the Rockstar halo finder by running a pipeline described in the second paper of this series (Rácz et al. in prep.). To establish the baseline for the comparison, in the spirit of the previous section, Fig. 4 shows a comparison of the cumulative halo mass function for the di ff erent codes in a Λ CDM cosmology. Since MG-Arepo and MG-Gadget are based on the same TreePM gravity solver, their agreement is excellent. They also agree very well with results from PKDGRAV3 which, as we like to remind the reader, are only available for the case of Λ CDM and are shown for reference here. On the other hand, ECOSMOG and ISIS are based on the AMR method. They agree with each other, but these simulations underestimate the abundance of low-mass ha-  \nbelow 10 13 h -1 M ⊙ . COLA uses a fixed-grid PM method, thus it also underestimates the abundance of low-mass halos. With relatively high resolution in these simulations ( N 1 / 3 mesh = 5 N 1 / 3 p ), COLA agrees well with ECOSMOG and ISIS . We note that PySCo exhibits a similar behaviour to COLA , ECOSMOG and ISIS , albeit with a lower amplitude. This discrepancy can be attributed to the fact that PySCo has the least small-scale clustering (see Fig. 1), resulting in a smaller number of halos within the simulation.  \nThe chameleon screening in f ( R ) depends on the halo mass such that high-mass halos are typically screened. The critical mass above which screening is e ff ective is determined by the parameter ¯ fR . As we can see from Fig. 5, the ratio of the halo mass function between f ( R ) and Λ CDM is enhanced for lowmass halos but approaches unity at the high-mass end where all halos are e ff ectively screened. The agreement between full simulations ( ECOSMOG and ISIS , MG-Arepo , MG-Gadget ) is around 4% for all halo masses. COLA uses an approximation for screening based on the thin-shell condition for a spherically symmetric object. Although this captures an overall e ff ect of screening, it fails for low-mass halos in F6 and high-mass halos in F5, leading to larger deviations.  \nFig. 3: Amplification factor of the matter power spectra from simulations of two nDGP scenarios, r c = 1 . 2 H -1 0 (N1.2, top row ) and r c = 5 . 6 H -1 0 (N5.6, bottom row ), relative to the reference Λ CDM cosmology. The amplification factor BP = P / P Λ CDM is measured for di ff erent codes (di ff erent line styles) at two di ff erent redshifts, z = 1 ( left panels ) and z = 0 . 667 ( right panels ). The bottom panel of each plot shows the relative agreement of the individual measurements, using MG-Arepo as a common reference.  \n<!-- image -->  \nIn the case of the Vainshtein mechanism, there is no halomass dependence on the screening. Despite screening being effective inside dark matter halos, these halos still feel enhanced gravitational attraction. This increases the merger rate and ultimately leads to larger enhancements of the halo mass function for halos of larger masses. In nDGP, shown in Fig. 6, all codes agree within 4%, with a sub-per cent agreement seen in the intermediate mass range 10 12 h -1 M ⊙ < M < 10 13 h -1 M ⊙ . ECOSMOG shows a large deviation for M ≲ 5 × 10 11 h -1 M ⊙ . For these masses, there are less than 50 dark matter particles assigned to each halo, which indicates that the results could be a ff ected by the refinement criteria of the simulations. The deviations are larger for the most massive halos, but the number of these halos is low and the halo mass function therefore becomes very noisy in this regime.', '4. Conclusions': "In this work, we have presented a comprehensive review of numerical methods for cosmological N -body simulations in scenarios extending beyond the standard Λ CDM model. Our ex-  \nloration spanned a variety of alternative DE and MG theories, highlighting the critical role of N -body simulations in connecting theoretical models with observational data. Through the detailed examination of numerical solvers and approximations tailored to these extended theories, we have showcased the state of the art of modelling the nonlinear scales of cosmic structure formation under a wide range of cosmological scenarios. Our code comparison exercise, based on the simulations from W15 and extended by incorporating new codes and approximation techniques, has demonstrated a fair consensus among di ff erent numerical implementations. This validation is particularly important for the Euclid mission, as the forthcoming observational data will require precise nonlinear modelling to constrain cosmological parameters e ff ectively.  \nThis article is part of a series that explores simulations and nonlinearities in models beyond Λ CDM. The simulation codes that we have considered in this article are used to generate simulation products in the companion paper by Rácz et al. (in preparation) and are crucial for generating the simulations needed to create the nonlinear modelling in the companion paper by Koyama et al. (in preparation), which forecasts constraints for  \nFig. 4: Cumulative halo mass function from simulations carried out with di ff erent codes (di ff erent line styles) for a Λ CDMcosmology. The left panels show the halo mass function at redshift z = 1 whereas the right panels show the halo mass function at redshift z = 0 . 667. The bottom panels show the relative di ff erence with respect to the simulation carried out with MG-Arepo .  \n<!-- image -->  \n/circledot  \n/circledot  \nf ( R ) gravity from photometric primary probes of the Euclid mission. The validation checks performed in this paper are therefore critical for being able to trust these results.  \nThe main outcomes of this work can be summarised as follows:  \n- ⋆ N -body simulation codes have been developed for a wide range of extended cosmological scenarios, ranging from simple DE models to more complex interacting scalar field models and MG theories, to non-standard dark matter and initial conditions; the availability of such codes will be a crucial asset for the future developments of large galaxy surveys such as Euclid ; we have provided a concise yet comprehensive overview of several such codes, their main features, implementation methods, assumptions and approximations.\n- ⋆ As a result of our validation e ff ort, we found agreement in the power spectrum boost at ≲ 1% up to k ≲ 1 h Mpc -1 and at ≲ 3% up to k ≲ 10 h Mpc -1 among all the codes implementing full field solvers ( MG-Arepo, MG-Gadget, ECOSMOG, ISIS, PySCo ), while approximate methods ( PANDA, MG-COLA ) display slightly larger deviations not exceeding 3% up to k ≲ 7 h Mpc -1 .\n- ⋆ As a matter of fact, among these simulation codes the ones involving algorithms for the solution of nonlinear di ff erential equations of some additional degrees of freedom, as e.g. for the case of MG theories, are the most challenging in terms of implementation and numerical convergence; we have therefore performed a thorough validation e ff ort of these methods through a code comparison study, extending the approach adopted in W15 to more recent and diverse algorithms.\n- ⋆ The halo mass function shows larger deviations among the codes, also due to larger di ff erences in the outcomes for the baseline Λ CDM simulations; nonetheless, all codes agree within less than 5% on the relative change of the halo mass function except for the very largest and the smallest mass ranges where poor statistics and insu ffi cient resolution, respectively, may impact the results.  \n- ⋆ Wealso compared the computational requirements of the different codes by measuring the CPU time needed to complete a reference MG simulation starting from identical initial conditions; we found that while full field solvers generally imply a substantial increase - up to a factor of ten - of the total CPU time relative to a Λ CDMsimulation, approximate solvers are not significantly more demanding for MG simulations compared to standard Λ CDMruns.  \nLooking forward, the continued evolution of simulation techniques will be paramount in leveraging the full potential of upcoming large-scale structure surveys such as Euclid . N -body simulations therefore continue to set a solid foundation for the next generation of cosmological inquiries. By persistently pushing the boundaries of computational astrophysics, we are poised to uncover the underlying physics driving the accelerated expansion of the Universe, thereby opening new windows onto the fundamental nature of DE, dark matter, and gravity itself.  \nAcknowledgements. The lead authors thank A. Schneider and R.E. Smith for their diligent work as internal referees and P. Schneider for carefully proofreading the manuscript. The work of JA is supported by the Swiss National Science Foundation. During part of this work, AMCLB was supported by a fellowship of PSL University hosted by the Paris Observatory. This project was provided with computer and storage resources by GENCI at CINES thanks to the grant 2023A0150402287 on the supercomputer Adastra's GENOA partition. GR's research was supported by an appointment to the NASA Postdoctoral Program administered by Oak Ridge Associated Universities under contract with NASA. GR was supported by JPL, which is run under contract by the California Institute of Technology for NASA (80NM0018D0004). The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid , in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft, funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d'Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space O ffi ce (SSO), and  \nFig. 5: Increase of the cumulative halo mass function from simulations for the two f ( R ) scenarios, ¯ f , R = 10 -5 (F5, top row ) and ¯ f , R = 10 -6 (F6, bottom row ), relative to the reference Λ CDM cosmology. The enhancement factor Bn = n / n Λ CDM is measured for di ff erent codes (di ff erent line styles) at two di ff erent redshifts, z = 1 ( left panels ) and z = 0 . 667 ( right panels ). The bottom panel of each plot shows the relative agreement of the individual measurements, using MG-Arepo as a common reference.  \n<!-- image -->  \n/circledot  \n/circledot  \nthe United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site ( http://www.euclid-ec.org ).", 'References': 'Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, Journal of High Energy Astrophysics, 34, 49  \nAdamek, J., Angulo, R. E., Arnold, C., et al. 2023, JCAP, 06, 035  \nAdamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016, Nature Phys., 12, 346  \nAdamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016, JCAP, 07, 053  \nAmendola, L. 2000, Phys. Rev. D, 62, 043511  \nAmendola, L. 2004, Phys. Rev. D, 69, 103524  \n- Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Rev. Rel., 21, 2\n- Amendola, L., Baldi, M., & Wetterich, C. 2008, Phys. Rev. D, 78, 023015\n- Armendariz-Picon, C., Mukhanov, V., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 85, 4438  \nArmendariz-Picon, C., Mukhanov, V., & Steinhardt, P. J. 2001, Phys. Rev. D, 63, 103510  \n- Arnold, C., Leo, M., & Li, B. 2019, Nature Astron., 3, 945\n- Arnold, C., Li, B., Giblin, B., Harnois-Déraps, J., & Cai, Y.-C. 2022, MNRAS, 515, 4161\n- Ayaita, Y., Baldi, M., Führer, F., Puchwein, E., & Wetterich, C. 2016, Phys. Rev. D, 93, 063511  \nBaldi, M. 2011, MNRAS, 411, 1077  \nBaldi, M. 2012a, MNRAS, 422, 1028  \nBaldi, M. 2012b, ASP Conf. Ser., 453, 155  \n- Baldi, M. 2023, MNRAS, 521, 613\n- Baldi, M., Pettorino, V., Robbers, G., & Springel, V. 2010, MNRAS, 403, 1684\n- Baldi, M. & Simpson, F. 2015, MNRAS, 449, 2239\n- Baldi, M. & Simpson, F. 2017, MNRAS, 465, 653\n- Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75\n- Barreira, A., Bose, S., & Li, B. 2015, JCAP, 12, 059\n- Barreira, A., Li, B., Hellwing, W. A., Baugh, C. M., & Pascoli, S. 2013, JCAP, 10, 027\n- Barreira, A., Li, B., Hellwing, W. A., Baugh, C. M., & Pascoli, S. 2014, JCAP, 09, 031\n- Becker, C., Arnold, C., Li, B., & Heisenberg, L. 2020, JCAP, 10, 055\n- Bose, S., Hellwing, W. A., & Li, B. 2015, JCAP, 02, 034\n- Bose, S., Li, B., Barreira, A., et al. 2017, JCAP, 02, 050\n- Brando, G., Fiorini, B., Koyama, K., & Winther, H. A. 2022, JCAP, 09, 051\n- Brando, G., Koyama, K., & Winther, H. A. 2023, JCAP, 06, 045\n- Brandt, A. 1977, Mathematics of Computation, 31, 333\n- Brax, P., Davis, A.-C., Li, B., Winther, H. A., & Zhao, G.-B. 2012, JCAP, 10, 002\n- Brax, P., Davis, A.-C., Li, B., Winther, H. A., & Zhao, G.-B. 2013, JCAP, 04, 029\n- Chan, K. C. & Scoccimarro, R. 2009, Phys. Rev. D, 80, 104005\n- Cheung, C., Fitzpatrick, A. L., Kaplan, J., Senatore, L., & Creminelli, P. 2008, JHEP, 3, 014\n- Chevallier, M. & Polarski, D. 2001, Int. J. Mod. Phys. D, 10, 213\n- Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1\n- Damour, T., Gibbons, G. W., & Gundlach, C. 1990, Phys. Rev. Lett., 64, 123  \nFig. 6: Increase of the cumulative halo mass function from simulations for the two nDGP scenarios, r c = 1 . 2 H -1 0 (N1.2, top row ) and r c = 5 . 6 H -1 0 (N5.6, bottom row ), relative to the reference Λ CDM cosmology. The enhancement factor Bn = n / n Λ CDM is measured for di ff erent codes (di ff erent line styles) at two di ff erent redshifts, z = 1 ( left panels ) and z = 0 . 667 ( right panels ). The bottom panel of each plot shows the relative agreement of the individual measurements, using MG-Arepo as a common reference.  \n<!-- image -->  \n- DESI Collaboration: Aghamousa, A., Aguilar, J., Ahlen, S., et al. 2016, arXiv:1611.00036\n- Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021a, Astropart. Phys., 131, 102605\n- Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021b, Astropart. Phys., 131, 102604\n- Dvali, G., Gabadadze, G., & Porrati, M. 2000, Phys. Lett. B, 485, 208\n- Elbers, W., Frenk, C. S., Jenkins, A., Li, B., & Pascoli, S. 2021, MNRAS, 507, 2614\n- Euclid Collaboration: Adam, R., Vannier, M., Maurogordato, S., et al. 2019, A&A, 627, A23\n- Euclid Collaboration: Ajani, V., Baldi, M., Barthelemy, A., et al. 2023, A&A, 675, A120\n- Euclid Collaboration: Bretonnière, H., Huertas-Company, M., Boucaud, A., et al. 2022, A&A, 657, A90\n- Euclid Collaboration: Bretonnière, H., Kuchner, U., Huertas-Company, M., et al. 2023, A&A, 671, A102\n- Euclid Collaboration: Castander, F. J., Fosalba, P., Stadel, J., et al. 2024 [ arXiv:2405.13495 ]\n- Euclid Collaboration: Castro, T., Fumagalli, A., Angulo, R. E., et al. 2023, A&A, 671, A100\n- Euclid Collaboration: Desprez, G., Paltani, S., Coupon, J., et al. 2020, A&A, 644, A31\n- Euclid Collaboration: Ilbert, O., de la Torre, S., Martinet, N., et al. 2021, A&A, 647, A117\n- Euclid Collaboration: Mellier, Y. et al. 2024 [ arXiv:2405.13491 ]  \nEuclid Collaboration: Merlin, E., Castellano, M., Bretonnière, H., et al. 2023, A&A, 671, A101  \n- Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112\n- Fiorini, B., Koyama, K., & Baker, T. 2023, JCAP, 12, 045\n- Giocoli, C., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813\n- Gleyzes, J., Langlois, D., & Vernizzi, F. 2014, Int. J. Mod. Phys. D, 23, 1443010\n- Gordon, J., de Aguiar, B. F., Rebouças, J. a., et al. 2024 [ arXiv:2404.12344 ]\n- Gronke, M. B., Llinares, C., & Mota, D. F. 2014, A&A, 562, A9\n- Guillet, T. & Teyssier, R. 2011, J. Comp. Phys., 230, 4756\n- Hagala, R., Llinares, C., & Mota, D. F. 2016, A&A, 585, A37\n- Hagala, R., Llinares, C., & Mota, D. F. 2017, Phys. Rev. Lett., 118, 101301\n- Hammami, A., Llinares, C., Mota, D. F., & Winther, H. A. 2015, MNRAS, 449, 3635\n- Harnois-Déraps, J., Hernandez-Aguayo, C., Cuesta-Lazaro, C., et al. 2023, MNRAS, 525, 6336\n- Hassani, F., Adamek, J., Kunz, M., & Vernizzi, F. 2019, JCAP, 12, 011\n- Hassani, F., L\'Huillier, B., Shafieloo, A., Kunz, M., & Adamek, J. 2020, JCAP, 04, 039\n- Hernández-Aguayo, C., Arnold, C., Li, B., & Baugh, C. M. 2021, MNRAS, 503, 3867\n- Hu, W. & Sawicki, I. 2007, Phys. Rev. D, 76, 064004\n- Ivezi\'c, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111\n- Joyce, A., Lombriser, L., & Schmidt, F. 2016, Ann. Rev. Nucl. Part. Sci., 66, 95 Khoury, J. & Wyman, M. 2009, Phys. Rev. D, 80, 064023\n- Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193\n- Lesgourgues, J. 2011, arXiv:1104.2932  \nLewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473  \n- Li, B. 2018, Simulating Large-Scale Structure for Models of Cosmic Acceleration, 2514-3433 (IOP Publishing)\n- Li, B., Barreira, A., Baugh, C. M., et al. 2013a, JCAP, 11, 012\n- Li, B. & Barrow, J. D. 2011, Phys. Rev. D, 83, 024007\n- Li, B., Zhao, G.-B., & Koyama, K. 2013b, JCAP, 05, 023\n- Li, B., Zhao, G.-B., Teyssier, R., & Koyama, K. 2012, JCAP, 01, 051\n- Linder, E. V. 2003, Phys. Rev. Lett., 90, 091301\n- Llinares, C. 2018, Int. J. Mod. Phys. D, 27, 1848003\n- Llinares, C., Hagala, R., & Mota, D. F. 2020, MNRAS, 491, 1868\n- Llinares, C. & Mota, D. 2013, Phys. Rev. Lett., 110, 161101\n- Llinares, C. & Mota, D. F. 2014, Phys. Rev. D, 89, 084023  \nLlinares, C., Mota, D. F., & Winther, H. A. 2014, A&A, 562, A78  \n- Lombriser, L. 2016, JCAP, 11, 039\n- Martinelli, M., Martins, C. J. A. P., Nesseris, S., et al. 2021, A&A, 654, A148  \nMauland, R., Winther, H. A., & Ruan, C.-Z. 2024, A&A, 685, A156  \n- Mitchell, M. A., Arnold, C., & Li, B. 2022, MNRAS, 514, 3349\n- Mota, D. F., Pettorino, V., Robbers, G., & Wetterich, C. 2008, Phys. Lett. B, 663, 160\n- Nesseris, S., Sapone, D., Martinelli, M., et al. 2022, A&A, 660, A67 Oyaizu, H. 2008, Phys. Rev. D, 78, 123523\n- Palma, D. & Candlish, G. N. 2023, MNRAS, 526, 1904\n- Pettorino, V. & Baccigalupi, C. 2008, Phys. Rev. D, 77, 103003\n- Pillepich, A. et al. 2018, MNRAS, 473, 4077\n- Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2  \nPourtsidou, A., Skordis, C., & Copeland, E. J. 2013, Phys. Rev. D, 88, 083505  \nPuchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348  \n- Ramachandra, N., Valogiannis, G., Ishak, M., & Heitmann, K. 2021, Phys. Rev. D, 103, 123525  \nRuan, C.-Z., Cuesta-Lazaro, C., Eggemeier, A., et al. 2024, MNRAS, 527, 2490 Ruan, C.-Z., Hernández-Aguayo, C., Li, B., et al. 2022, JCAP, 05, 018  \n- Sáez-Casares, I. n., Rasera, Y., & Li, B. 2023 [ arXiv:2303.08899 ]\n- Sawicki, I. & Bellini, E. 2015, Phys. Rev. D, 92, 084061\n- Schmidt, F. 2009a, Phys. Rev. D, 80, 123003\n- Schmidt, F. 2009b, Phys. Rev. D, 80, 043001\n- Schneider, A., Teyssier, R., Potter, D., et al. 2016, JCAP, 04, 047\n- Simpson, F. 2010, Phys. Rev. D, 82, 083505\n- Skordis, C., Pourtsidou, A., & Copeland, E. J. 2015, Phys. Rev. D, 91, 083537\n- Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757\n- Springel, V. 2005, MNRAS, 364, 1105\n- Springel, V. 2010, MNRAS, 401, 791\n- Tassev, S., Zaldarriaga, M., & Eisenstein, D. 2013, JCAP, 06, 036\n- Teyssier, R. 2002, A&A, 385, 337\n- Tsujikawa, S. 2010, Lect. Notes Phys., 800, 99\n- Vainshtein, A. I. 1972, Phys. Lett. B, 39, 393\n- Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291\n- Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32\n- Wesseling, P. 2004, An Introduction to Multigrid Methods (Philadelphia: R.T. Edwards)\n- Wetterich, C. 1995, A&A, 301, 321\n- Wetterich, C. 2004, Phys. Lett. B, 594, 17\n- Will, C. M. 2014, Living Rev. Rel., 17, 4\n- Winther, H. A. & Ferreira, P. G. 2015a, Phys. Rev. D, 91, 123507\n- Winther, H. A. & Ferreira, P. G. 2015b, Phys. Rev. D, 92, 064005\n- Winther, H. A., Koyama, K., Manera, M., Wright, B. S., & Zhao, G.-B. 2017, JCAP, 08, 006\n- Winther, H. A., Schmidt, F., Barreira, A., et al. 2015, MNRAS, 454, 4208\n- Wright, B. S., Winther, H. A., & Koyama, K. 2017, JCAP, 10, 054\n- Zumalacárregui, M., Bellini, E., Sawicki, I., Lesgourgues, J., & Ferreira, P. G. 2017, JCAP, 08, 019\n- 1 Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland\n- 3 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93 / 2, 40129 Bologna, Italy\n- 2 Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK\n- 4 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93 / 3, 40129 Bologna, Italy\n- 6 Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Muhlenberg 1, D-14476 Potsdam-Golm, Germany\n- 5 INFN-Sezione di Bologna, Viale Berti Pichat 6 / 2, 40127 Bologna, Italy\n- 7 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s / n, 08193 Barcelona, Spain\n- 8 Institut de Ciencies de l\'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s / n Cerdanyola del Vallés, 08193 Barcelona, Spain\n- 10 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway\n- 9 Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France\n- 11 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA\n- 13 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK\n- 12 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany\n- 14 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK\n- 16 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy\n- 15 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy\n- 17 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK\n- 19 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany\n- 18 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK\n- 20 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France\n- 22 Institut de Physique Théorique, CEA, CNRS, Université ParisSaclay 91191 Gif-sur-Yvette Cedex, France\n- 21 Université St Joseph; Faculty of Sciences, Beirut, Lebanon\n- 23 School of Mathematics and Physics, University of Surrey, Guildford, Surrey, GU2 7XH, UK\n- 25 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy\n- 24 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy\n- 26 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy\n- 28 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy\n- 27 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy\n- 29 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy\n- 31 INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 30 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 32 Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126, Napoli, Italy\n- 34 INFN section of Naples, Via Cinthia 6, 80126, Napoli, Italy\n- 33 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy\n- 35 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal\n- 37 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France\n- 36 Faculdade de Ciências da Universidade do Porto, Rua do Campo de Alegre, 4150-007 Porto, Portugal\n- 38 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy\n- 40 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy\n- 39 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy\n- 41 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c / o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy\n- 43 Port d\'Informació Científica, Campus UAB, C. Albareda s / n, 08193 Bellaterra (Barcelona), Spain\n- 42 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain  \n- 44 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany\n- 45 Institut d\'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain\n- 46 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6 / 2, 40127 Bologna, Italy\n- 47 Instituto de Astrofísica de Canarias, Calle Vía Láctea s / n, 38204, San Cristóbal de La Laguna, Tenerife, Spain\n- 48 European Space Agency / ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy\n- 49 ESAC / ESA, Camino Bajo del Castillo, s / n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain\n- 50 Université Claude Bernard Lyon 1, CNRS / IN2P3, IP2I Lyon, UMR 5822, Villeurbanne, F-69100, France\n- 51 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland\n- 52 Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, 08028 Barcelona, Spain\n- 53 Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig de Lluís Companys 23, 08010 Barcelona, Spain\n- 54 UCB Lyon 1, CNRS / IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France\n- 55 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal\n- 56 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal\n- 57 Department of Astronomy, University of Geneva, ch. d\'Ecogia 16, 1290 Versoix, Switzerland\n- 58 Université Paris-Saclay, CNRS, Institut d\'astrophysique spatiale, 91405, Orsay, France\n- 59 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy\n- 60 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy\n- 61 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France\n- 62 FRACTAL S.L.N.E., calle Tulipán 2, Portal 13 1A, 28231, Las Rozas de Madrid, Spain\n- 63 INAF-Osservatorio Astronomico di Padova, Via dell\'Osservatorio 5, 35122 Padova, Italy\n- 64 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany\n- 65 Universitäts-Sternwarte München, Fakultät für Physik, LudwigMaximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany\n- 66 Dipartimento di Fisica "Aldo Pontremoli", Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy\n- 67 Felix Hormuth Engineering, Goethestr. 17, 69181 Leimen, Germany\n- 68 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark\n- 69 Cosmic Dawn Center (DAWN), Denmark\n- 70 Université Paris-Saclay, CNRS / IN2P3, IJCLab, 91405 Orsay, France\n- 71 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany\n- 72 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA\n- 73 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 74 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland\n- 75 Aix-Marseille Université, CNRS / IN2P3, CPPM, Marseille, France\n- 76 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, CH1211 Genève 4, Switzerland\n- 77 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland  \n- 78 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland\n- 80 Centre de Calcul de l\'IN2P3 / CNRS, 21 avenue Pierre de Coubertin 69627 Villeurbanne Cedex, France\n- 79 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands\n- 81 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany\n- 83 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France\n- 82 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93 / 2, 40129 Bologna, Italy\n- 84 University of Applied Sciences and Arts of Northwestern Switzerland, School of Engineering, 5210 Windisch, Switzerland\n- 86 Institut d\'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France\n- 85 Institut d\'Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris, France\n- 87 Institut de Física d\'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain\n- 89 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, Denmark\n- 88 European Space Agency / ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands\n- 90 Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada\n- 92 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada\n- 91 Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada\n- 93 Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy\n- 95 Institute of Space Science, Str. Atomistilor, nr. 409 M˘agurele, Ilfov, 077125, Romania\n- 94 Centre National d\'Etudes Spatiales - Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France\n- 96 Dipartimento di Fisica e Astronomia "G. Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy\n- 98 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25 / 8, 6020 Innsbruck, Austria\n- 97 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile\n- 99 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain\n- 101 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain\n- 100 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal\n- 102 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands\n- 104 Dipartimento di Fisica, Università degli studi di Genova, and INFN-Sezione di Genova, via Dodecaneso 33, 16146, Genova, Italy\n- 103 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy\n- 105 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA\n- 107 Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA), Loc. Lignan 39, I-11020, Nus (Aosta Valley), Italy\n- 106 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy\n- 108 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK\n- 110 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France\n- 109 School of Physics and Astronomy, Cardi ff University, The Parade, Cardi ff , CF24 3AA, UK  \n- 111 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy\n- 113 CERCA / ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA\n- 112 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain\n- 114 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy\n- 116 Departamento de Astrofísica, Universidad de La Laguna, 38206, La Laguna, Tenerife, Spain\n- 115 Departamento de Física Fundamental. Universidad de Salamanca. Plaza de la Merced s / n. 37008 Salamanca, Spain\n- 117 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy\n- 119 Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan\n- 118 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy\n- 120 Ludwig-Maximilians-University, Schellingstrasse 4, 80799 Munich, Germany\n- 122 Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy\n- 121 Max-Planck-Institut für Physik, Boltzmannstr. 8, 85748 Garching, Germany\n- 123 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA\n- 125 Université Côte d\'Azur, Observatoire de la Côte d\'Azur, CNRS, Laboratoire Lagrange, Bd de l\'Observatoire, CS 34229, 06304 Nice cedex 4, France\n- 124 Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands\n- 126 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA\n- 128 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary\'s University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada\n- 127 Department of Physics & Astronomy, University of California Irvine, Irvine CA 92697, USA\n- 129 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain\n- 131 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK\n- 130 Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica, Universidad de La Laguna (ULL), 38200, La Laguna, Tenerife, Spain\n- 132 CEA Saclay, DFR / IRFU, Service d\'Astrophysique, Bat. 709, 91191 Gif-sur-Yvette, France\n- 134 Instituto de Astrofísica de Canarias, c / Via Lactea s / n, La Laguna E-38200, Spain. Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco Sanchez, La Laguna, E-38200, Spain\n- 133 Department of Computer Science, Aalto University, PO Box 15400, Espoo, FI-00 076, Finland\n- 135 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany\n- 137 Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of Turku, Finland\n- 136 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000, Grenoble, France\n- 138 Serco for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 140 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, Victoria 3122, Australia\n- 139 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia\n- 141 School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK\n- 142 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town, 7535, South Africa\n- 143 ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil\n- 145 Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm, SE-106 91, Sweden\n- 144 IRFU, CEA, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France\n- 146 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK\n- 148 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy\n- 147 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125, Firenze, Italy\n- 149 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal\n- 151 INFN, Sezione di Roma 2, Via della Ricerca Scientifica 1, Roma, Italy\n- 150 Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, Roma, Italy\n- 152 HE Space for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 154 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy, and INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 153 Aurora Technology for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 155 Theoretical astrophysics, Department of Physics and Astronomy, Uppsala University, Box 515, 751 20 Uppsala, Sweden\n- 157 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK\n- 156 Department of Physics, Royal Holloway, University of London, TW20 0EX, UK\n- 158 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA\n- 160 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark\n- 159 Cosmic Dawn Center (DAWN)\n- 161 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA\n- 162 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA  \nFig. A.1: Comparison of the computational costs of simulations run under di ff erent gravity models (labels on the x -axis) and with di ff erent codes (markers and colours as described in the legend).  \n<!-- image -->', 'Appendix A: Computational cost': 'One of the obstacles in producing accurate predictions in the nonlinear regime of structure formation is presented by the computational cost of running large N -body simulations. This issue becomes even more pronounced in the case of MG simulations, due to the additional computational cost associated with the solution of the Klein-Gordon equation for the scalar field describing the additional degree of freedom in these theories. To provide some insights into the trade-o ff s between accuracy and time-tosolution, we attempt here a comparative analysis of the computational cost of the simulations run for this paper. Given that these simulations were run on di ff erent machines and with di ff erent parallelisation settings, we cannot conduct a precise assessment of their computational cost. Instead, we limit this analysis to an order-of-magnitude comparison of the simulations run with the various codes and models discussed in this paper.  \nFor this exercise, we use the information on the wall-clock time, T real, and the number of cores, N cores, as recorded in the log files of the simulation runs. This information was available for all simulations except for the ones run with ISIS . We estimate the computational cost of the simulations, C , as the product of the wall-clock time and the number of cores:  \nC ≡ N cores T real . (A.1)  \nThe performance of multigrid codes depends on the convergence criterion chosen. For PySCo , we chose an extremely con-  \nThe estimates of the computational cost for the simulations are compared in Fig. A.1. We can see that the cost of Λ CDM simulations of Tree-PM and AMR codes is C ∼ 10 3 CPUh, while for MG-COLA and PySCo the cost is about one order of magnitude lower at C ∼ 10 2 CPUh. For f ( R ) gravity models instead, the computational cost increases significantly (approximately by a factor of ten) for the codes that solve the full Klein-Gordon equation of the scalar field, namely MG-Arepo , MG-Gadget , ECOSMOG , and PySCo , while the overhead is smaller for the codes that adopt screening approximations, namely PANDA and MG-COLA . Finally, in nDGP gravity, only ECOSMOG has a significant overhead compared to Λ CDM, while MG-Arepo and the approximate codes, PANDA and MG-COLA , have just a small overhead.  \nservative approach (see Sect. 3.1.5) with a very low tolerance threshold, resulting in more V-cycles and almost double the CPU time needed to solve the linear Poisson equation compared to a more standard setup. Regarding the nonlinear solver, we use two F-cycles instead of a single one (which should in principle be enough, but it is not the goal of the present paper to provide a convergence study), therefore roughly doubling the CPU time needed for the f ( R ) gravity models.  \nWe stress that a thorough assessment of the e ffi ciency of the codes is beyond the scope of this paper and would have required a much more methodical e ff ort including (but not limited to)  \n- -running the simulations in a controlled environment,\n- -testing the scaling performance of each code.\n- -conducting convergence tests for the various hyperparameters,  \nIn fact, when focusing only on predictions of the amplification factors, it is possible to achieve a similar level of accuracy with lower force, mass or time resolution, since resolution effects mostly cancel out when taking ratios of quantities a ff ected by the same inaccuracies (Brando et al. 2022). This has been shown to be the case for MG-COLA simulations in Fiorini et al. (2023), where the use of a lower resolution allowed accurate predictions of power spectrum boosts in nDGP gravity with a theoretical gain of ∼ 300 with respect to the computational cost of the COLA simulations presented here.  \nSuch large speed-ups have paved the way for creating emulators for the nonlinear amplification of the power spectrum in models beyond Λ CDM in a cost-e ff ective way, i.e. without the need for supercomputers. This has already been done for some of the models we consider in this paper (see e.g. Ramachandra et al. 2021; Mauland et al. 2024; Fiorini et al. 2023). For instance, Fiorini et al. (2023) found that an emulator for the nDGP model can be constructed with as little as a few thousand CPUh worth of computational time. Likewise, Mauland et al. (2024) who presented a generic pipeline for using COLA to create such emulators, used f ( R ) gravity as an example and found similar numbers for the required computational time. Gordon et al. (2024) described a simulation setup that can also be used for emulating the full power spectrum (up to a reasonable high wavenumber k ∼ 1 h Mpc -1 ), requiring around ∼ 100 CPUh per simulation on a modern CPU. Emulators have also been constructed for beyondΛ CDM models using high-resolution direct simulations in the same way as has been done for Λ CDM. This approach generally gives more accurate emulators than those created with approximate methods, but this comes at a much higher cost. For example, both Sáez-Casares et al. (2023) and Arnold et al. (2022) have each presented a high-fidelity emulator for the f ( R ) model considered in this paper, but at a higher cost of about 3.5-4 million CPUh.'}