© The Authors 2024

1. Introduction

The Brans–Dicke (BD) theory (Brans & Dicke 1961; Dicke 1962) of modified gravity (MG; Clifton et al. 2012) features a scalar field in addition to the spacetime metric. It effectively promotes Newton’s constant in general relativity (GR) to a dynamical gravitational strength. This scalar–tensor theory was introduced in the 1960s to implement Mach’s principle in GR. Since then it has become perhaps the most famous alternative to Einstein’s theory and a root for the development of other theories, such as the more general Horndeski theory (Horndeski 1974). BD features an additional parameter ω and reduces to GR as ω → ∞, so BD can be understood as a generalization of GR. In cosmology, replacing GR with BD in the standard Λ cold dark matter ((GR-)ΛCDM) model gives the alternative BD-ΛCDM model, where the effective gravitational strength evolves with time and modifies the expansion history and growth of structure.

As summarized in Table 1, observations in the 21st century have placed very strong constraints on ω. Already in 2003, Shapiro time delay measurements of the Cassini satellite on its way to Saturn constrained ω ≳ 10⁴, and recent strong-field tests from timing rapidly rotating neutron stars (pulsars) even bound ω ≳ 10⁵. On cosmological scales, however, the story is different. Roughly speaking, constraints have tightened from only ω ≳ 10² with stage-II survey data in the 2000s to ω ≳ 10³ with stage-III surveys in the 2010s. In the next decade, Alonso et al. (2017) expects cosmological constraints to improve by another order of magnitude with upcoming stage-IV surveys like Euclid (Laureijs et al. 2011), DESI (DESI Collaboration 2016), the Vera C. Rubin Observatory (LSST Science Collaboration 2009), SKA (Yahya et al. 2015) and next-generation cosmic microwave background (CMB) experiments (Abazajian et al. 2016). Fisher forecasts found that if ω = 800, then Euclid should be able to constrain this to 3% using galaxy clustering and weak lensing data (Frusciante et al. 2024). Likewise, Ballardini et al. (2019) forecast that upcoming clustering and weak lensing data in combination with BOSS and CMB observations have the potential to reach ω ≳ 10⁴ in the most optimistic case. With these seemingly ever-tightening competing constraints, BD has become a benchmark theory for testing (deviations from) GR. BD also survived the burial of theories with gravitational wave speeds v ≠ c from GW170817 (Ezquiaga & Zumalacárregui 2017), and Solà Peracaula et al. (2019, 2020) show that it alleviates the H₀ and S₈ tensions (Perivolaropoulos & Skara 2022).

In fact, several scalar–tensor theories reduce to BD on large (cosmological) scales, where gradients of the scalar field are suppressed (Avilez & Skordis 2014; Joudaki et al. 2022). This makes BD interesting not only on its own, but also as the large-scale limit of more general theories, so constraints on it have wider implications. Such theories typically feature screening mechanisms that hide the scalar field in dense small-scale regions, letting GR take over there. In other words, the true theory of gravity could be similar to BD on the largest scales, similar to GR on the smallest scales, and transition between them over intermediate scales. The very tight small-scale (astrophysical) constraints could therefore be irrelevant from a cosmological perspective.

In this work, we aim to develop a tool for predicting the nonlinear BD-ΛCDM matter power spectrum P_BD(k ≲ 5 h/Mpc,  z ≲ 3 | θ) across parameter space θ to around 1% precision, for use in constraining the theory with upcoming stage-IV large-scale structure survey data. There are several possible approaches to this, all of which rely on N-body simulations to account for the nonlinear evolution.

Halo modeling tools like HMCODE (Mead et al. 2015, 2021) and HALOFIT (Smith et al. 2003; Takahashi et al. 2012) rely on dark matter’s clustering in halos to effectively reverse-engineer the nonlinear matter distribution from halo statistics. They can be integrated directly with linear Einstein-Boltzmann solvers. The HMCODE framework has been extended to account for a wide range of dark energy and MG models, massive neutrinos and baryonic physics (Mead et al. 2016). Another related method is the halo model reaction approach of Cataneo et al. (2019), which combines the halo model and perturbation theory to model corrections coming from nonstandard physics. This is implemented in the publicly available REACT code (Bose et al. 2020; Giblin et al. 2019). The halo modeling approach was taken in Joudaki et al. (2022), who used N-body simulations to create a fitting formula in HMCODE to predict the nonlinear BD-ΛCDM power spectrum. As far as we know, this is the only nonlinear prediction for BD so far.

Emulators are trained from a set of (time-consuming) N-body simulations for a limited number of cosmological parameters, and then rely on machine learning to (quickly) interpolate in parameter space and output the nonlinear matter power spectrum for arbitrary cosmologies. There are already sophisticated emulators for the GR-ΛCDM power spectrum, such as EUCLIDEMULATOR2 (Euclid Collaboration 2021), COSMICEMU (Heitmann et al. 2016; Lawrence et al. 2017; Moran et al. 2023), and BACCO (Angulo et al. 2021; Aricò et al. 2021, 2022). Recent works, including Winther et al. (2019), Brando et al. (2022), Fiorini et al. (2023), SESAME (Mauland et al. 2024), and E-MANTIS (Sáez-Casares et al. 2024), have produced emulators for selected MG theories. It is also possible to combine halo modeling with emulation and emulate the ingredients of the halo model (Ruan et al. 2024). At the time of writing, no dedicated BD-ΛCDM power spectrum emulator exists.

Simulation rescaling (Angulo & White 2010) is a technique where data from a single N-body simulation is rescaled to produce output for a different cosmology. This technique was in fact used in training the BACCO emulator.

We demonstrate a hybrid technique, somewhat reminiscent of both rescaling and emulation, suited to extending existing predictors of GR-ΛCDM to BD-ΛCDM. By carefully selecting the cosmological parameters θ_BD and θ_GR, we map a given BD universe to a corresponding GR universe such that their cosmological evolutions are as similar as possible. This makes it possible to predict the nonlinear matter power spectrum boost $B=P_{\mathrm{BD}}/P_{\mathrm{GR}}\approx P_{\mathrm{BD}}^{\mathrm{lin}}/P_{\mathrm{GR}}^{\mathrm{lin}}$ using linear power spectra from cheap Einstein-Boltzmann solvers, which we verified by comparing it to the boost obtained from N-body simulations. In turn, we can predict the full BD-ΛCDM power spectrum P_BD = B ⋅ P_GR by combining it with any existing high-quality emulator for P_GR. This “trick” exploits the effort that has already gone into creating sophisticated GR emulators, and avoids duplicating this work for every alternative model to be explored. It could be a viable cheap route for constraining MG theories with upcoming surveys (Frusciante et al. 2024; Casas et al. 2017).

The rest of this paper is structured as follows. Section 2 reviews the background cosmology and perturbation theory of BD (and thus GR). Section 3 describes our pipeline for the Einstein-Boltzmann solver and N-body simulations. Section 4 justifies the parameterization we use to compute the boost. Section 5 shows the resulting linear and nonlinear boost and compares them to older predictions. Section 6 concludes.

2. Brans–Dicke modified gravity

Brans–Dicke (BD) theory (Brans & Dicke 1961; Dicke 1962) generalizes general relativity (GR). In addition to the metric tensor g_μν, it features a scalar field ϕ with a constant dimensionless parameter ω. Here we use units where c = 1, but explicitly keep Newton’s constant G = 6.67 × 10⁻¹¹ Nm²/kg² to show how BD effectively modifies the gravitational constant of GR. In the Jordan frame, the total action of BD gravity minimally coupled to the matter action S_M is

$$\begin{array}{c}\hfill S[g_{\mu \nu },\varphi ]=\frac{1}{16\pi G}{\displaystyle \int {\mathrm{d}}^4}x\sqrt{-|g|}[\varphi R-\frac{\omega }{\varphi }{(\mathrm{\nabla }\varphi )}^2]+S_{\mathrm{M}},\end{array}$$(1)

where g_μν is the spacetime metric with determinant |g|< 0, Ricci scalar R and covariant derivative ∇_μ; and (∇ϕ)² = (∇^μϕ)(∇_μϕ) = g^μν(∇_μϕ)(∇_νϕ). From the principle of least action δS = 0, the (classical) equations of motion for g_μν and ϕ are the modified Einstein and Klein-Gordon field equations (Clifton et al. 2012)

$$\begin{array}{cc}\hfill G_{\mu \nu }& =\frac{8\pi G}{\varphi }{T}_{\mu \nu }+\frac{\omega }{{\varphi }^2}[({\mathrm{\nabla }}_{\mu }\varphi )({\mathrm{\nabla }}_{\nu }\varphi )-\frac{1}{2}{g}_{\mu \nu }{(\mathrm{\nabla }\varphi )}^2]\hfill \\ \hfill & +\frac{1}{\varphi }[{\mathrm{\nabla }}_{\mu }{\mathrm{\nabla }}_{\nu }\varphi -g_{\mu \nu }{\mathrm{\nabla }}^2\varphi ],\hfill \end{array}$$(2a)

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}^2\varphi & =\frac{8\pi G}{3+2\omega }T,\hfill \end{array}$$(2b)

where ∇² = g^μν∇_μ∇_ν is the four-dimensional Laplacian, and $T_{\mu \nu }=(-2/\sqrt{-g})(\delta S_{\mathrm{M}}/\delta g^{\mu \nu })$ is the energy-momentum tensor of matter with trace $T=T_{\mu }^{\mu }$. Effectively, ϕ promotes Newton’s constant G in the Einstein field equations (2a) to a dynamical field G/ϕ. In the following, we see in more detail how BD is similar to GR with such an effective gravitational strength.

GR is recovered in the limit ω → ∞: then ∇²ϕ = 0, whose solution ϕ = 1 everywhere recovers Einstein’s field equations G_μν = 8πG T_μν with Newton’s constant G.

2.1. Cosmology

We assume a spatially flat, homogeneous and isotropic Friedmann–Lemaître–Robertson–Walker background with small perturbations in the Newtonian gauge:

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-[1+2\mathrm{\Psi }(t,\mathit{x})]\mathrm{d}{t}^2+a^2(t)[1+2\mathrm{\Phi }(t,\mathit{x})]\mathrm{d}{\mathit{x}}^2.\end{array}$$(3)

We fill the universe with a perfect fluid with total energy-momentum tensor

$$\begin{array}{c}\hfill T^{\mu \nu }={\displaystyle \sum _{\mathrm{s}}}{T}_{(s)}^{\mu \nu }={\displaystyle \sum _{\mathrm{s}}}({\rho }_{\mathrm{s}}+P_{\mathrm{s}})u_{\mathrm{s}}^{\mu }{u}_{\mathrm{s}}^{\nu }+P_{\mathrm{s}}{g}^{\mu \nu },\end{array}$$(4)

with energy densities ρ_s, pressures P_s, and four-velocities $u_{\text{s}}^{\mu }$ from the species s:

• radiation r (from photons γ and massless¹ neutrinos ν) with equation of state ${\overline{P}}_{\mathrm{r}}={\overline{\rho }}_{\mathrm{r}}/3$,

• matter m (from baryons b and cold dark matter c) with equation of state ${\overline{P}}_{\mathrm{m}}=0$,

• dark energy Λ (from a cosmological constant) with equation of state ${\overline{P}}_{\mathrm{\Lambda }}=-{\overline{\rho }}_{\mathrm{\Lambda }}$.

The universe’s background and (linear) perturbations are described perturbatively up to first order around a homogeneous and isotropic universe.

2.2. Background cosmology

To zeroth order, the equation of motions (2) give the Friedmann and scalar field equations

$$\begin{array}{cc}\hfill H^2& =\frac{8\pi G}{3\overline{\varphi }}\overline{\rho }+\frac{\omega }{6}{\left(\frac{\dot{\overline{\varphi }}}{\overline{\varphi }}\right)}^2-H\frac{\dot{\overline{\varphi }}}{\overline{\varphi }},{(}_{\mu \nu }{=}_{00})\hfill \end{array}$$(5a)

$$\begin{array}{cc}\hfill \frac{1}{a^3}\frac{\mathrm{d}}{\mathrm{d}t}\left(\dot{\overline{\varphi }}{a}^3\right)& =\frac{8\pi G}{3+2\omega }(\overline{\rho }-3\overline{P}),\hfill \end{array}$$(5b)

where $\overline{\rho }={\overline{\rho }}_{\mathrm{r}}+{\overline{\rho }}_{\mathrm{m}}+{\overline{\rho }}_{\mathrm{\Lambda }}$, $\overline{P}={\overline{P}}_{\mathrm{r}}+{\overline{P}}_{\mathrm{m}}+{\overline{P}}_{\mathrm{\Lambda }}$, ˙ = d/dt and H = ȧ/a is the Hubble parameter. The Friedmann equation in BD is the same as in GR, only with additional effective energy from $\dot{\overline{\varphi }}$ and an effective gravitational strength $G_{\mathrm{H}}=G/\overline{\varphi }$.

According to the Boltzmann equation, the three species evolve as ${\overline{\rho }}_{\mathrm{r}}={\overline{\rho }}_{\mathrm{r}0}{a}^{-4}$, ${\overline{\rho }}_{\mathrm{m}}={\overline{\rho }}_{\mathrm{m}0}{a}^{-3}$ and ${\overline{\rho }}_{\mathrm{\Lambda }}={\overline{\rho }}_{\mathrm{\Lambda }0}{a}^0$, where ₀ indexes quantities by their values today. Defining the standard density parameters ${\mathrm{\Omega }}_{\mathrm{s}}={\overline{\rho }}_{\mathrm{s}}/(3H^2/8\pi G)$, physical densities today are measured by ${\omega }_{\mathrm{s}0}={\mathrm{\Omega }}_{\mathrm{s}0}{h}^2\propto {\overline{\rho }}_{\mathrm{s}0}$². Re-parameterizing the evolution in terms of ′=d/dlna = H⁻¹ ⋅ d/dt, the system (5) can be rewritten in the dimensionless form

$$\begin{array}{cc}\hfill E^2={\left(\frac{H}{H_0}\right)}^2& =\frac{1}{\overline{\varphi }}·\frac{{\displaystyle {\mathrm{\Omega }}_{\mathrm{r}0}{a}^{-4}+{\mathrm{\Omega }}_{\mathrm{m}0}{a}^{-3}+{\mathrm{\Omega }}_{\mathrm{\Lambda }0}}}{1-\frac{\omega }{6}{[{(ln\overline{\varphi })}^{\prime }]}^2+{(ln\overline{\varphi })}^{\prime }},\hfill \end{array}$$(6a)

$$\begin{array}{cc}\hfill ({\overline{\varphi }}^{\prime }E{a}^3{)}^{\prime }& =\frac{3}{3+2\omega }\frac{{\mathrm{\Omega }}_{\mathrm{m}0}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }0}{a}^{-3}}{E}.\hfill \end{array}$$(6b)

This system can be integrated from an initial scale factor a_ini, given $\overline{\varphi }(a_{\mathrm{ini}})$, $\overline{\varphi }\prime (a_{\mathrm{ini}})$, Ω_r0, Ω_m0, Ω_Λ0, and ω. But these are not all free parameters, as today’s Friedmann equation E₀ = 1 constrains one of them – say Ω_Λ0. To find a consistent solution, one can iterate over Ω_Λ0 and solve the system repeatedly until this constraint is satisfied.

The integration can be simplified by selecting a_ini in the radiation era. There H ∝ a⁻² and $\overline{\rho }-3\overline{P}\simeq {\overline{\rho }}_{\mathrm{r}}-3{\overline{P}}_{\mathrm{r}}=0$, so the Klein–Gordon equation (5b) gives $\dot{\overline{\varphi }}\propto a^{-3}$ with the solution $\overline{\varphi }=A/a+B$. The diverging solution is unphysical, leaving a frozen $\overline{\varphi }=\mathrm{const}$ with

$$\begin{array}{c}\hfill \dot{\overline{\varphi }}={\dot{\overline{\varphi }}}_{\mathrm{ini}}=0\text{(during radiation domination)}.\end{array}$$(7)

In the Poisson equation (9) below, we see that the effective gravitational strength felt by matter today depends on ω and ${\overline{\varphi }}_0$. We therefore used the shooting method to vary ${\overline{\varphi }}_{\mathrm{ini}}$ and integrate the system repeatedly until we hit the desired ${\overline{\varphi }}_0$.

2.3. Linear perturbations

Ignoring anisotropic relativistic stresses and polarization and working in the quasi-static limit $\{\dot{\mathrm{\Phi }},\dot{\mathrm{\Psi }},\dot{\delta \varphi }\}=H·\{\mathrm{\Phi },\mathrm{\Psi },\delta \varphi \}$ and sub-horizon limit k ≫ H (Orjuela-Quintana & Nesseris 2023), the linearly perturbed Boltzmann, Einstein and Klein–Gordon equations in k-space are (Solà Peracaula et al. 2019)

$$\begin{array}{cc}\hfill {({\delta }_{\mathrm{c}}+3\mathrm{\Phi })}^{\prime }& =\frac{k}{\mathit{aH}}{v}_{\mathrm{c}},\hfill \end{array}$$(8a)

$$\begin{array}{cc}\hfill {({\delta }_{\mathrm{b}}+3\mathrm{\Phi })}^{\prime }& =\frac{k}{\mathit{aH}}{v}_{\mathrm{b}},\hfill \end{array}$$(8b)

$$\begin{array}{cc}\hfill {({\delta }_{\gamma }+4\mathrm{\Phi })}^{\prime }& =\frac{4}{3}\frac{k}{\mathit{aH}}{v}_{\gamma },\hfill \end{array}$$(8c)

$$\begin{array}{cc}\hfill {(a{v}_{\mathrm{c}})}^{\prime }& =-\frac{k}{H}\mathrm{\Psi },\hfill \end{array}$$(8d)

$$\begin{array}{cc}\hfill {(a{v}_{\mathrm{b}})}^{\prime }& =-\frac{k}{H}\mathrm{\Psi }+\frac{4}{3}\frac{{\overline{\rho }}_{\gamma }}{{\overline{\rho }}_{\mathrm{b}}}(v_{\mathrm{b}}-v_{\gamma })a{\tau }^{\prime },\hfill \end{array}$$(8e)

$$\begin{array}{cc}\hfill {(a{v}_{\gamma })}^{\prime }& =-\frac{k}{H}\mathrm{\Psi }-\frac{k}{4H}{\delta }_{\gamma }+a{v}_{\gamma }+(v_{\gamma }-v_{\mathrm{b}})a{\tau }^{\prime },\hfill \end{array}$$(8f)

$$\begin{array}{cc}\hfill k^2(\mathrm{\Phi }+\frac{1}{2}\frac{\delta \varphi }{\overline{\varphi }})& =4\pi \frac{G}{\overline{\varphi }}{a}^2\delta \rho ,\hfill \end{array}$$(8g)

$$\begin{array}{cc}\hfill k^2\frac{\delta \varphi }{\overline{\varphi }}& =\frac{8\pi }{2\omega +3}\frac{G}{\overline{\varphi }}{a}^2(\delta \rho -3\delta P),\hfill \end{array}$$(8h)

$$\begin{array}{cc}\hfill \mathrm{\Psi }+\mathrm{\Phi }& =-\frac{\delta \varphi }{\overline{\varphi }}.\hfill \end{array}$$(8i)

Motion of matter is sensitive to Ψ through the geodesic equation. In matter domination, δP ≃ 0 and $\delta \rho \simeq {\overline{\rho }}_{\mathrm{m}}{\delta }_{\mathrm{m}}\equiv {\overline{\rho }}_{\mathrm{b}}{\delta }_{\mathrm{b}}+{\overline{\rho }}_{\mathrm{c}}{\delta }_{\mathrm{c}}$, and solving equations (8g), (8h) and (8i) for Ψ gives its Poisson equation

$$\begin{array}{c}\hfill -k^2\mathrm{\Psi }=4\pi \underset{G_{\mathrm{m}}}{\underbrace{\frac{2\omega +4}{2\omega +3}\frac{G}{\overline{\varphi }}}}\frac{{\overline{\rho }}_{\mathrm{m}0}}{a}{\delta }_{\mathrm{m}}=\frac{3}{2}\frac{G_{\mathrm{m}}}{G}\frac{{\mathrm{\Omega }}_{\mathrm{m}0}{H}_0^2}{a}{\delta }_{\mathrm{m}}.\end{array}$$(9)

The Poisson equation in BD is the same as in GR, only with an effective gravitational strength $G_{\mathrm{m}}=(2\omega +4)/(2\omega +3)(G/\overline{\varphi })=(2\omega +4)/(2\omega +3)G_{\mathrm{H}}$ felt by matter that depends on time through the scalar field in the background. Cavendish experiments today measure G_m0, so the scalar field today must be ${\overline{\varphi }}_0=(2\omega +4)/(2\omega +3)$ in “restricted models” with the correct Newtonian limit G_m0 = G. However, we also allowed for “unrestricted models” with arbitrary G_m0 and ${\overline{\varphi }}_0$, so we can constrain G_m0 with cosmology (Zahn & Zaldarriaga 2003; Ballardini et al. 2022).

To find a simplified growth equation for the total matter overdensity δ_m, baryon–photon interactions at late times can be neglected. Combining the Boltzmann equations (8a), (8b), (8d) and (8e) with the Poisson equation (9) then gives the growth equation

$$\begin{array}{cc}\hfill 0& ={\ddot{\delta }}_{\mathrm{m}}+2H{\dot{\delta }}_{\mathrm{m}}-4\pi G_{\mathrm{m}}{\overline{\rho }}_{\mathrm{m}}{\delta }_{\mathrm{m}}\hfill \\ \hfill & =H^2\{{\delta }_{\mathrm{m}}^{″}+[2+{(lnE)}^{\prime }]{\delta }_{\mathrm{m}}^{\prime }-\frac{3}{2}\frac{G_{\mathrm{m}}}{G}{\mathrm{\Omega }}_{\mathrm{m}}(a){\delta }_{\mathrm{m}}\}.\hfill \end{array}$$(10)

The general solution of this second-order differential equation is a linear combination δ_m(k, a) = C₊(k)D₊(a)+C₋(k)D₋(a) of a growing mode D₊(a) and a decaying mode D₋(a). The growing mode is called the scale-independent growth factor D₊(a), and, neglecting the decaying mode, relates D₊(a₂)/D₊(a₁) = δ_m(k, a₂)/δ_m(k, a₁) for any k.

On the other hand, deflection of light through weak lensing is sensitive to the Weyl potential (Ψ − Φ)/2. To find its source equation, eliminate δϕ from equations (8g), (8h), and (8i) to get

$$\begin{array}{c}\hfill -k^2\left(\frac{\mathrm{\Phi }-\mathrm{\Psi }}{2}\right)=4\pi \frac{G}{\overline{\varphi }}{\overline{\rho }}_{\mathrm{m}0}{a}^{-1}{\delta }_{\mathrm{m}}.\end{array}$$(11)

Here we see that it is the effective $G_{\gamma }=G/\overline{\varphi }$ that affects light trajectories. In general, we encountered three effective gravitational strengths G_γ = G_H ≠ G_m that affect light trajectories, the expansion and matter clustering, respectively. Although they generally take on two different values, we have G_γ = G_H ≈ G_m for (viable) models with appreciable ω.

The growth of structure is summarized by the matter power spectrum

$$\begin{array}{c}\hfill P(k,z)=P_{\mathrm{prim}}(k)·{|{\mathrm{\Delta }}_{\mathrm{m}}(k,z)|}^2,\end{array}$$(12a)

of the gauge-invariant density perturbation Δ_m ≈ δ_m amplified from a nearly scale-invariant (n_s ≲ 1) primordial power spectrum

$$\begin{array}{c}\hfill P_{\mathrm{prim}}(k)=\frac{2{\pi }^2}{k^3}·A_{\mathrm{s}}{\left(\frac{k}{k_{\mathrm{p}}}\right)}^{n_{\mathrm{s}}-1}\end{array}$$(12b)

with primordial amplitude A_s and spectral index n_s.

2.4. Nonlinear structure formation

To perform N-body simulations in BD we only need the geodesic equation

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathit{x}}{\mathrm{d}{t}^2}+2H\frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}=-\frac{1}{a^2}\mathbf{\nabla }\mathrm{\Psi },\end{array}$$(13)

and the Poisson equation

$$\begin{array}{c}\hfill {\mathbf{\nabla }}^2\mathrm{\Psi }=\frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m}0}\frac{G_{\mathrm{m}}}{G}{H}_0^2{a}^{-1}{\delta }_{\mathrm{m}}.\end{array}$$(14)

The geodesic equation is the same in BD and GR, and the Poisson equation is the same as the linear version (9). The only two modifications from a GR N-body simulation are the modified Hubble function H(a) and the effective gravitational strength G_m(a), making it trivial to extend any N-body code from GR to BD.

Figure 1 compares the background evolution in two BD-ΛCDM and GR-ΛCDM universes with cosmological parameters to be introduced in Sect. 4 and Table 2.

Fig. 1. Evolution of the effective gravitational parameter G, Hubble function E(a) = H(a)/H0, growth factor D(a), and growth rate f(a) = d lnD(a)/d lna from equations (5) and (10) with the scale factor a in the fiducial BD-ΛCDM and GR-ΛCDM cosmologies in Table 2 (equivalent to the transformation in Fig. 5a). The common time of matter-radiation equality is marked by ${\overline{\rho }}_{\mathrm{r}}={\overline{\rho }}_{\mathrm{m}}$.

3. Computational pipeline

To compute the matter power spectrum P(k, z | θ) in a BD-ΛCDM or GR-ΛCDM universe, we joined three codes into the pipeline in Fig. 2:

• HI_CLASS (Zumalacárregui et al. 2017; Bellini et al. 2020), based on CLASS (Blas et al. 2011), solves the background and perturbation Einstein-Boltzmann equations in GR-ΛCDM and several MG cosmologies, including BD-ΛCDM. The perturbed equations it solves are more sophisticated than our simplified presentation (8) and too lengthy to reproduce here. It outputs the linear matter power spectrum P_{HI_CLASS}(k,z). Bellini et al. (2018) showed that its BD matter power spectrum agrees to 0.1% with four other codes.

• FML³ is a particle-mesh N-body code that uses the COLA method (Tassev et al. 2013) in GR-ΛCDM and several MG cosmologies, including BD-ΛCDM. It supersedes the code MG-PICOLA (Winther et al. 2017), in turn based on L-PICOLA (Howlett et al. 2015). The code first backscales the linear power spectrum (see e.g. Fidler et al. (2017)) from HI_CLASS, and then seeds initial particle positions. It then solves the Poisson equation (14) with the time-dependent G_m(t) in BD and Newton’s constant G_m = G in GR, Finally, it outputs the (quasi-)nonlinear matter power spectrum P_COLA(k, z) by constructing the density field from the particle position snapshots. These results are used for our main analysis. We used a box with volume L³ = (384 h/Mpc)³ and N_part = N_cell = 512³ particles and cells, evolving from redshift z_init = 10 until today with N_step = 30 time steps. These simulations sacrifice accuracy at the deepest nonlinear scales for overall computational speed: the particle-mesh nature of the code does not resolve small-scale forces, and the COLA method maintains accuracy up to quasilinear scales with few time steps. This is an excellent trade-off in our scenario. To some extent, the error in P at smaller scales even cancels after we compute the boost P_BD/P_GR.

• RAMSES (Teyssier 2002) solves the Poisson equation with the standard N-body method with adaptive mesh refinement and time-stepping in GR-ΛCDM and (with our modifications) BD-ΛCDM. We modified its GR version⁴ to solve the same Poisson equation as FML using arbitrary splined functions for H(a) and G_m(a) over time, passed on from HI_CLASS. It starts from the same initial particle positions at the same redshift z_init as FML, uses the same box size L and number of particles N_part, adaptively refines a base grid with the same size N_cell, but uses adaptive time stepping that ignores FML’s N_step. The snapshots are analyzed in the same way as FML, except that we quadrupled the number of grid cells along each dimension due to the finer resolution. This also outputs the nonlinear matter power spectrum P_RAMSES(k,z), and is only used to validate FML’s results.

Fig. 2. Our BD- or GR-specific simulation pipeline begins from a set of input parameters θBD or θGR. First, they are passed to the Einstein-Boltzmann solver HI_CLASS, which solves the background and perturbations and outputs the linear power spectrum. Next, FML back-scales PHI_CLASS(k,z) from z = 0 to z = zinit, realizes initial conditions for an N-body simulation with the random seed s and evolves it forward with the COLA method on a uniform grid, outputting the nonlinear power spectrum. Finally, FML’s first snapshot can be passed on to RAMSES, which evolves the same initial conditions with the standard N-body method on an adaptive mesh and yields an independent nonlinear power spectrum.

All codes branch to solve the equations specific to BD-ΛCDM or GR-ΛCDM at all times. To compensate for slight code differences, such as time stepping and splining effects, we normalized P from FML and RAMSES to match that from HI_CLASS at the most linear common scale k = 0.02 h/Mpc.

3.1. Boost computation

To compute the matter power spectrum boost

$$\begin{array}{c}\hfill B(k_{\mathrm{BD}},z|{\mathit{\theta }}_{\mathrm{BD}})=\frac{P_{\mathrm{BD}}(k_{\mathrm{BD}},z|{\mathit{\theta }}_{\mathrm{BD}})}{P_{\mathrm{GR}}(k_{\mathrm{GR}},z|{\mathit{\theta }}_{\mathrm{GR}})},\end{array}$$(15)

our pipeline first runs one BD simulation with freely chosen parameters θ_BD, followed by one GR simulation with parameters θ_GR = θ_GR(θ_BD) decided by a transformation of the BD parameters. It then compares the power spectra at equal redshifts z, but generally different wavenumbers k_BD and k_GR = k_GR(k_BD), which we justify in Sect. 4.

To minimize cosmic variance within each (BD, GR) simulation pair, it uses a common initial condition seed s_BD = s_GR = s (for fixed θ_BD and θ_GR). But the seed s = s(θ_BD) is changed for every new set of parameters to avoid bias toward one particular configuration of the universe. To further reduce cosmic variance, we also used amplitude-fixed initial conditions (Angulo & Pontzen 2016; Villaescusa-Navarro et al. 2018).

3.2. Tests and convergence analysis

HI_CLASS and FML independently calculate E(a) = H(a)/H₀, Ω_Λ0, and σ₈ in GR, and also $\overline{\varphi }(a)$ and $\overline{\varphi }\prime (a)$ in BD. We verified that these agree within a small tolerance after every run.

Figure 3 shows a convergence analysis of the results from FML. The results change by less than 1% when the increasing the resolution of the simulation parameters, showing that the results from our setup have converged.

Fig. 3. Convergence of FML’s boost (15) at z = 0 in the fiducial cosmology in Table 2 (equivalent to the transformation in Fig. 5a), as computational parameters are varied from the fiducial Npart = 5123, Ncell = 5123, Nstep = 30, zinit = 10, and L = 384 Mpc/h.

To test the absolute accuracy of the particle-mesh simulations we ran using FML, we also compared it to results from the adaptive-mesh-refinement code RAMSES. These results are shown in Fig. 4 and again, the results are equal up to 1% error up to k ≲ 5 h/Mpc. Together, these tests show that the results of the boost from FML are both precise and accurate to 1% error. We only performed these tests for the fiducial cosmology in Table 2: we chose a small value ω = 100 in BD to maximize the deviation from GR, and assume the results can be trusted for other cosmological parameters.

Fig. 4. Comparison of the matter power spectrum (boost) obtained from FML using the COLA method on a uniform grid and RAMSES using the standard N-body method on an adaptive grid, in the fiducial cosmology in Table 2 (equivalent to the transformation in Fig. 5a). To compensate for small code differences, we have rescaled all spectra to match the linear HI_CLASS spectrum at their smallest common k. We cut RAMSES’ spectra at increasingly higher k ≤ {1/2, 1, 2}⋅kNyquist as the mesh refinement increases for z = {3, 1.5, 0}. Although P from FML and RAMSES disagree, that error largely cancels in their boosts B = PBD/PGR, which agree to 1% up to k ≲ 5 h/Mpc.

4. Model parameterization

What parameter and wavenumber transformations θ_GR(θ_BD) and k_GR(k_BD) should we use in the boost (15)? That is, how should we compare a BD universe to a GR universe? The answer may sound as obvious as the question sounds pedantic: that we should use the same parameters in BD and GR. But this is vague or impossible, as it depends on our selection of cosmological parameters, and all cosmological parameters cannot be equal in two different universes. For example, the background equations (5) forbid having the same ${\overline{\rho }}_{\mathrm{r}0}$, ${\overline{\rho }}_{\mathrm{m}0}$, ${\overline{\rho }}_{\mathrm{\Lambda }0}$, H₀ and ${\overline{\varphi }}_0$ today. For us this is a question of convenience. We can choose to compare two universes in the boost (15) such that this ratio behaves as simple and predictable as possible. At the end of the day, we want to multiply it by P_GR to predict P_BD, anyway.

Next, we discuss some transformation alternatives while testing them in Fig. 5.

Fig. 5. Matter power spectrum boost (15) with different parameter transformations θGR = θGR(θBD), where ($A_{\mathrm{s}}^{\mathrm{BD}}=A_{\mathrm{s}}^{\mathrm{GR}}$ or ${\sigma }_8^{\text{BD}}={\sigma }_8^{\text{GR}}$) and ($H_0^{\text{BD}}=H_0^{\text{GR}}$ or $H_0^{\mathrm{GR}}=H_0^{\mathrm{BD}}{\overline{\varphi }}_{\mathrm{ini}}^{1/2}$). Requiring equal σ8 (instead of equal As) flattens the boost in the nonlinear regime, while requiring equal ϕinih2 (instead of equal h) synchronizes the sound horizon and avoids rapid baryon acoustic oscillations.

4.1. Wavenumbers and Hubble parameter

Consider the perturbations (8) and forget about recombination for a moment, which enters through the optical depth τ. In the Boltzmann equations, k and H appear in the particular combination k/H. At early times, the same is true for the Einstein and Klein-Gordon equations after substituting the effective $G_{\mathrm{H}}=G/\overline{\varphi }=3H^2/8\pi \overline{\rho }$ from the Friedmann equation (5a) with the frozen scalar field (7). For large enough ω that sufficiently suppress the anisotropic stress $\delta \varphi /\overline{\varphi }$, perturbations in BD and GR with equal k/H (or $k/\sqrt{G_{\mathrm{H}}}$) evolve similarly.

In other words, BD resembles GR with a different expansion rate H due to a different effective gravitational strength. This difference can be factored out by comparing modes with equal k/H (or $k/\sqrt{G_{\mathrm{H}}}$). Similarly, (Zahn & Zaldarriaga 2003, equation (12)) shows that a mode k′ in a GR universe with gravitational constant G′ evolves like a mode $k=k\prime /\sqrt{G\prime /G}$ in another GR universe with Newton’s constant G. In other words, as they elegantly explained it, “gravity introduces no preferred scale, so the dynamics of the perturbations remains the same when scales are measured in units of the expansion time.” This also holds for BD’s scale-independent effective gravitational strength.

However, it breaks when recombination is activated and τ is included, as the faster expansion hinders (re)combination of H⁺ + e⁻ and dampens the small-scale acoustic peaks in the CMB spectrum. This is indeed how Zahn & Zaldarriaga (2003) shows how the CMB constrains G to percent-level. However, most matter is dark, so this effect is much smaller for the power spectrum of matter that concerns us here.

This motivates matching modes {k_BD, k_GR} such that k/H is the same at early times. To be compatible with N-body codes that factor H₀ out of the relevant equations, we can compare modes with equal k/H₀ as long as we make E_BD = E_GR at early times. The frozen scalar field (7) cuts the Friedmann equation (5a) off at $H_{\mathrm{BD}}^2=H_{\mathrm{GR}}^2/\overline{\varphi }$ provided we use the same physical densities ${\overline{\rho }}_{\mathrm{s}0}\propto {\omega }_{\mathrm{s}0}$ in BD and GR⁵ so we can accomplish this by transforming

$$\begin{array}{c}\hfill H_0^{\mathrm{GR}}=H_0^{\mathrm{BD}}{\overline{\varphi }}_{\mathrm{ini}}^{1/2}⟺H_0^{\mathrm{GR}}/\sqrt{[b]G_{\mathrm{GR}}^{\mathrm{ini}}}=H_0^{\mathrm{BD}}/\sqrt{[b]G_{\mathrm{BD}}^{\mathrm{ini}}}.\end{array}$$(16)

This changes the simulation box size L = 384 Mpc/h with H₀ = 100 h km/s/Mpc, so we can match modes at different scales. Our “scale scaling” $k_{\mathrm{GR}}=(H_0^{\mathrm{GR}}/H_0^{\mathrm{BD}})k_{\mathrm{BD}}={\overline{\varphi }}_{\mathrm{ini}}^{1/2}{k}_{\mathrm{BD}}$ contributes to our goal of “mapping” a BD universe onto a GR universe.

Another insightful way to understand this rescaling is to consider the sound horizon

$$\begin{array}{c}\hfill s={\displaystyle {\int }_0^t}{c}_{\mathrm{s}}\frac{\mathrm{d}t}{a}={\displaystyle {\int }_0^a}\sqrt{\frac{4{\overline{\rho }}_{\gamma }/3{\overline{\rho }}_{\mathrm{m}}}{1+4{\overline{\rho }}_{\gamma }/3{\overline{\rho }}_{\mathrm{m}}}}\frac{\mathrm{d}a}{H},\end{array}$$(17)

associated with baryonic acoustic oscillations (BAO). This length scale “freezes in” when the baryons and photons decouple at recombination, leaving a signature in the matter power spectrum at (multiples of) the characteristic wavenumber k_s ∝ 1/s. Indeed, the Hubble parameter transformation (16) ensures that k_s/h will be the same in the two universes. This synchronizes the phases of the oscillations in P_BD(k_BD) and P_GR(k_GR), avoiding oscillations in their ratio P_BD/P_GR. Figure 5 shows the suppressed oscillations when parameterized by ${\overline{\varphi }}_{\mathrm{ini}}{h}^2$ instead of h.

The key to understand is that matching modes by k/h combined with the parameter transformation (16) leverages the frozen scalar field (7) during radiation domination to make the initial evolution of perturbations as similar as possible in BD and GR. This only holds until matter domination, when the scalar field begins to move. But modes that have entered the horizon by this time simply grow (approximately) with the scale-independent growth factor. These are precisely the interesting smaller-scale modes beyond the peak of the matter power spectrum! Larger-scale modes are perfectly described by linear perturbation theory, anyway.

4.2. Power spectrum normalization

It is common to parameterize the primordial power spectrum (12b) with A_s. However, as we are interested in predicting the late-time matter power spectrum, it is more natural to ask for “equal power” in BD and GR today than at early times. The normalization of today’s power spectrum is conventionally done in terms of the parameter

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{R}}^2=\frac{1}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}{P}_0(k){\widehat{W}}_{\mathrm{R}}^2(k)k^2\mathrm{d}k,\end{array}$$(18)

where ${\widehat{W}}_{\mathrm{R}}(k)=3[sin(kR)-kRcos(kR)]/{(kR)}^3$ is the Fourier transform of the top-hat profile W_R(r) = Θ(R − r)/(4πR³/3) with radius R. It is common to define the normalization using the value of σ₈, i.e. σ_R with R = 8 Mpc/h. By pinning ${\sigma }_8^{\text{BD}}={\sigma }_8^{\text{GR}}$ today⁶, we get less primordial power $A_{\mathrm{s}}^{\mathrm{BD}}<A_{\mathrm{s}}^{\mathrm{GR}}$ to compensate for the increased structure growth in BD due to both less Hubble friction E_BD < E_GR and stronger gravity G_BD > G_GR (see Fig. 1). Figure 5 shows that using the same A_s gives a peak in the boost, while using the same σ₈ gives a flatter boost.

This shows that the parameterization in Fig. 5a with σ₈ and ${\overline{\varphi }}_{\mathrm{ini}}{h}^2$ gives the simplest and most predictable nonlinear boost. As an added bonus, it also deviates the least from its linear counterpart, and we show below that this holds through parameter space.

5. Resulting matter power spectrum boost

We now focus on the boost with the particular parameterization developed in Sect. 4, defaulting to the fiducial parameters in Table 2 (equivalent to the transformation in Fig. 5a). In particular, due to the cosmological constraints ω ≳ 1000 summarized in Table 1, we restricted ourselves to ω ≥ 100 and used the minimum ω = 100 in the fiducial cosmology.

5.1. Understanding the linear boost growth

The simplest way to understand the evolution of the linear scale-dependent boost is in terms of the scale-independent growth factor D₊(a) from equation (10). During matter domination, the power spectrum scales like $P\propto {\mathrm{\Delta }}_{\text{m}}^2\approx {\delta }_{\text{m}}^2\propto D_{+}^2$ at sub-horizon scales, so

$$\begin{array}{c}\hfill B(k,z)\propto {\left(\frac{D_{+}^{\mathrm{BD}}(a(z))}{D_{+}^{\mathrm{GR}}(a(z))}\right)}^2.\end{array}$$(19)

Fig. 6 compares this to the evolution of the full linear boost from early matter domination till today. As seen in Fig. 1, the scalar field in BD begins to move during matter domination, resulting in both less Hubble friction $E_{\text{BD}}^{\prime }<E_{\text{GR}}^{\prime }$ and stronger gravity G_BD > G_GR. This gives faster growth in BD!

Fig. 6. Evolution of the linear matter power spectrum boost from HI_CLASS compared to the scale-independent growth (19) in Fig. 1, starting from a = 10−3 during early matter domination. The boosts are calculated in the fiducial cosmology in Table 2 (equivalent to the transformation in Fig. 5a). Notice that the curves for z = 0 and z = 9 overlap, and that the boost turns around at the peak z = 1.5 in Fig. 1.

When dark energy takes over, Fig. 1 shows that the expansion accelerates more in BD than in GR with $E_{\text{BD}}^{\prime }>E_{\text{GR}}^{\prime }$, while the gravitational strengths G_BD → G_GR also become equal. This causes the boost to peak around z = 1.5 and then fall off.

5.2. Cosmological dependence of the nonlinear boost

We now study how the nonlinear boost changes with the cosmological parameters in Table 2. In Fig. 7, we vary one by one parameter away from its fiducial value while keeping the others fixed. Not only is the boost (divided by ${(h_{\mathrm{BD}}/h_{\mathrm{GR}})}^{-3}={\overline{\varphi }}_{\mathrm{ini}}^{3/2}$) very resistant toward parameter changes, but so is its resemblance $P_{\mathrm{BD}}/P_{\mathrm{GR}}\approx P_{\mathrm{BD}}^{\mathrm{lin}}/P_{\mathrm{GR}}^{\mathrm{lin}}$ of the linear boost within 1% up to k ≤ 1 h/Mpc, and 2% up to k ≤ 5 h/Mpc.

Fig. 7. Variation of the matter power spectrum boost (15) at z = 0 as each cosmological parameter is varied away from its fiducial value in Table 2 (equivalent to the transformation in Fig. 5a). The bottom right plot is special; showing the ratio of the nonlinear boost from FML/COLA to the linear boost from HI_CLASS for all boosts in the 7 other plots at more redshifts z = 0.0, 1.5, 3.0. It shows that B ≈ Blin to 1% for k ≤ 1 h/Mpc, and 2% for k ≤ 5 h/Mpc.

5.3. Linear prediction of the nonlinear boost

Conveniently, this lets us get away with predicting the nonlinear boost $B\approx B_{\mathrm{lin}}=P_{\mathrm{BD}}^{\mathrm{lin}}/P_{\mathrm{GR}}^{\mathrm{lin}}$ using linear theory and codes, like HI_CLASS. In turn, we can predict the nonlinear power spectrum P_BD = B_lin ⋅ P_GR from an existing nonlinear prediction of P_GR. In fact, existing P_GR emulators like EUCLIDEMULATOR2 only output the nonlinear correction factor $C=P_{\mathrm{GR}}/P_{\mathrm{GR}}^{\mathrm{lin}}$ and defer multiplication with $P_{\mathrm{GR}}^{\mathrm{lin}}$, so it even looks like we can cancel $P_{\mathrm{GR}}^{\mathrm{lin}}$ and simply calculate $P_{\mathrm{BD}}({\mathit{\theta }}_{\mathrm{BD}})=C({\mathit{\theta }}_{\mathrm{GR}}({\mathit{\theta }}_{\mathrm{BD}}))·P_{\mathrm{BD}}^{\mathrm{lin}}({\mathit{\theta }}_{\mathrm{BD}})$. However, we still need to calculate $P_{\mathrm{GR}}^{\mathrm{lin}}$ in order to calculate the σ₈ integral (18) from the perturbations to make the parameter transformation θ_GR(θ_BD).

We provide the program BD.PY⁷ that handles these subtleties and predicts the nonlinear BD spectrum from 2 runs of HI_CLASS for the linear spectra BD and GR spectra and 1 run of EUCLIDEMULATOR2 for the nonlinear GR spectrum. This is our main product, and can be used in future fitting to large-scale structure survey data.

5.4. Comparison with existing nonlinear prediction

To demonstrate our program BD.PY, we compared it to the prediction of Joudaki et al. (2022) that we mentioned in Sect. 1. They incorporated BD into HMCODE by tuning the expression for the virial overdensity to match a set of N-body simulations. For their purpose of using KiDS data to constrain the model, the accuracy aim was only a few percent up to k ≲ 1 h/Mpc and 5 − 10% up to k ≲ 10 h/Mpc. Fig. 8 shows that their prediction was (almost) within their tolerance, but produces too little power in BD when trialed with a tighter 1% tolerance.

Fig. 8. A nonlinear BD-ΛCDM power spectrum (boost) with ω = 100 obtained with our script BD.PY compared to that from an older fitting-formula modification to HMCODE in HI_CLASS from Joudaki et al. (2022), and COLA N-body results from FML.

6. Conclusion

Despite increasingly tight constraints from small-scale experiments, BD gravity remains interesting at cosmological scales, where constraints are weaker and it approximates more general theories of gravity with screening that falls back to GR at small scales. More precise data from upcoming stage-IV large-scale structure surveys, like Euclid, can increase the competitiveness of cosmological constraints, and shed new light on this situation to illuminate the path forward for constraints on gravity.

In fact, this work began as an effort to train a traditional machine learning emulator for the nonlinear BD matter power spectrum from a set of COLA N-body simulations. Instead, we found that it is possible to predict the nonlinear boost B = P_BD/P_GR using linear theory and codes if one transforms the cosmological parameters θ_BD → θ_GR in a particular way. In detail, by comparing the linear boost with the nonlinear boosts from both COLA and standard N-body simulations, we verified that

$$\begin{array}{c}\hfill \frac{P_{\mathrm{BD}}(k/h_{\mathrm{BD}},z)}{P_{\mathrm{GR}}(k/h_{\mathrm{GR}},z)}\stackrel{2\%}{\simeq }\frac{P_{\mathrm{BD}}^{\mathrm{lin}}(k/h_{\mathrm{BD}},z)}{P_{\mathrm{GR}}^{\mathrm{lin}}(k/h_{\mathrm{GR}},z)}\end{array}$$

$$\begin{array}{c}\hfill \text{for}\left\{\begin{array}{c}k\lesssim 5h/\mathrm{Mpc}\\ z\lesssim 3\end{array}\right\}\text{if}\left\{\begin{array}{c}{\omega }_{\mathrm{BD}}\gtrsim 100\\ {\sigma }_8^{\mathrm{GR}}={\sigma }_8^{\mathrm{BD}}\\ H_0^{\mathrm{GR}}=H_0^{\mathrm{BD}}{\overline{\varphi }}_{\mathrm{ini}}^{1/2}\\ {\overline{\rho }}_{m0}^{\mathrm{GR}}={\overline{\rho }}_{m0}^{\mathrm{BD}}\\ {\overline{\rho }}_{r0}^{\mathrm{GR}}={\overline{\rho }}_{r0}^{\mathrm{BD}}\end{array}\right\}.\end{array}$$

This “trick” exploits the boost’s dependence on cosmological parameters to map the nonlinear boost to the linear boost.

This has the advantage that the boost can be computed with cheap linear codes like HI_CLASS, instead of expensive nonlinear N-body methods. When paired with an existing nonlinear P_GR predictor, such as EUCLIDEMULATOR2, this paves a simple and efficient path for predicting P_BD and constraining BD theory with the precision of upcoming stage-IV large-scale structure surveys. We provide the program BD.PY for this.

A drawback of our approach is that the parameter transformation θ_BD → θ_GR requires the emulator P_GR(θ_GR) to cover a broader range of cosmological parameters than that motivated by GR alone. For example, for reasonable input parameters h_BD and $A_{\mathrm{s}}^{\mathrm{BD}}$, our transformation can request values of h_GR and $A_{\mathrm{s}}^{\mathrm{GR}}$ that are outside the parameter ranges covered by a GR emulator. To facilitate MG studies, we therefore encourage developers of base GR emulators to “think big” when setting parameter bounds. It would also be valuable to develop a traditional matter power spectrum emulator dedicated to BD, or improve the precision of previous halo modeling approaches.

We also encourage others to investigate whether such simplifying parameter mappings θ_MG → θ_GR exist for other MG theories. Our findings could be specific to BD gravity or apply to more theories, perhaps with a similar scale-independent nature. Even if one does not find $P_{\mathrm{MG}}({\mathit{\theta }}_{\mathrm{MG}})/P_{\mathrm{GR}}({\mathit{\theta }}_{\mathrm{GR}})\approx P_{\mathrm{MG}}^{\mathrm{lin}}({\mathit{\theta }}_{\mathrm{MG}})/P_{\mathrm{GR}}^{\mathrm{lin}}({\mathit{\theta }}_{\mathrm{GR}})$ “exactly,” investing thought in a clever parameter transformation may still pay off by simplifying the shape of P_MG/P_GR so that one can obtain it from $P_{\mathrm{MG}}^{\mathrm{lin}}/P_{\mathrm{GR}}^{\mathrm{lin}}$ through a fitting formula or ease the training of emulators. This reduces the effort needed to explore alternatives to ΛCDM with MG, so one does not have to duplicate the effort gone into GR for every such alternative.

Acknowledgments

We thank the Research Council of Norway for their support under grant no. 287772.

References

All Tables

All Figures

Fig. 1. Evolution of the effective gravitational parameter G, Hubble function E(a) = H(a)/H0, growth factor D(a), and growth rate f(a) = d lnD(a)/d lna from equations (5) and (10) with the scale factor a in the fiducial BD-ΛCDM and GR-ΛCDM cosmologies in Table 2 (equivalent to the transformation in Fig. 5a). The common time of matter-radiation equality is marked by ${\overline{\rho }}_{\mathrm{r}}={\overline{\rho }}_{\mathrm{m}}$.

In the text

Fig. 2. Our BD- or GR-specific simulation pipeline begins from a set of input parameters θBD or θGR. First, they are passed to the Einstein-Boltzmann solver HI_CLASS, which solves the background and perturbations and outputs the linear power spectrum. Next, FML back-scales PHI_CLASS(k,z) from z = 0 to z = zinit, realizes initial conditions for an N-body simulation with the random seed s and evolves it forward with the COLA method on a uniform grid, outputting the nonlinear power spectrum. Finally, FML’s first snapshot can be passed on to RAMSES, which evolves the same initial conditions with the standard N-body method on an adaptive mesh and yields an independent nonlinear power spectrum.

In the text

	Fig. 3. Convergence of FML’s boost (15) at z = 0 in the fiducial cosmology in Table 2 (equivalent to the transformation in Fig. 5a), as computational parameters are varied from the fiducial Npart = 5123, Ncell = 5123, Nstep = 30, zinit = 10, and L = 384 Mpc/h.
In the text

Fig. 4. Comparison of the matter power spectrum (boost) obtained from FML using the COLA method on a uniform grid and RAMSES using the standard N-body method on an adaptive grid, in the fiducial cosmology in Table 2 (equivalent to the transformation in Fig. 5a). To compensate for small code differences, we have rescaled all spectra to match the linear HI_CLASS spectrum at their smallest common k. We cut RAMSES’ spectra at increasingly higher k ≤ {1/2, 1, 2}⋅kNyquist as the mesh refinement increases for z = {3, 1.5, 0}. Although P from FML and RAMSES disagree, that error largely cancels in their boosts B = PBD/PGR, which agree to 1% up to k ≲ 5 h/Mpc.

In the text

Fig. 5. Matter power spectrum boost (15) with different parameter transformations θGR = θGR(θBD), where ($A_{\mathrm{s}}^{\mathrm{BD}}=A_{\mathrm{s}}^{\mathrm{GR}}$ or ${\sigma }_8^{\text{BD}}={\sigma }_8^{\text{GR}}$) and ($H_0^{\text{BD}}=H_0^{\text{GR}}$ or $H_0^{\mathrm{GR}}=H_0^{\mathrm{BD}}{\overline{\varphi }}_{\mathrm{ini}}^{1/2}$). Requiring equal σ8 (instead of equal As) flattens the boost in the nonlinear regime, while requiring equal ϕinih2 (instead of equal h) synchronizes the sound horizon and avoids rapid baryon acoustic oscillations.

In the text

Fig. 6. Evolution of the linear matter power spectrum boost from HI_CLASS compared to the scale-independent growth (19) in Fig. 1, starting from a = 10−3 during early matter domination. The boosts are calculated in the fiducial cosmology in Table 2 (equivalent to the transformation in Fig. 5a). Notice that the curves for z = 0 and z = 9 overlap, and that the boost turns around at the peak z = 1.5 in Fig. 1.

In the text

Fig. 7. Variation of the matter power spectrum boost (15) at z = 0 as each cosmological parameter is varied away from its fiducial value in Table 2 (equivalent to the transformation in Fig. 5a). The bottom right plot is special; showing the ratio of the nonlinear boost from FML/COLA to the linear boost from HI_CLASS for all boosts in the 7 other plots at more redshifts z = 0.0, 1.5, 3.0. It shows that B ≈ Blin to 1% for k ≤ 1 h/Mpc, and 2% for k ≤ 5 h/Mpc.

In the text

	Fig. 8. A nonlinear BD-ΛCDM power spectrum (boost) with ω = 100 obtained with our script BD.PY compared to that from an older fitting-formula modification to HMCODE in HI_CLASS from Joudaki et al. (2022), and COLA N-body results from FML.
In the text