11institutetext: Institute of Theoretical Astrophysics, University of Oslo,
P.O.Box 1029 Blindern, N-0315 Oslo, Norway
11email<EMAIL_ADDRESS>
# A simple prediction of the non-linear matter power spectrum in Brans-Dicke
gravity from linear theory

Herman Sletmoen 11 Hans A. Winther 11

Brans-Dicke (BD) was one of the first proposed scalar-tensor theories of
gravity, and effectively turns the gravitational constant of General
Relativity (GR) time-dependent. Constraints on the BD parameter $\omega$ serve
as a benchmark for testing GR, which is recovered in the limit
$\omega\rightarrow\infty$. Current small-scale astrophysical constraints
$\omega\gtrsim 10^{5}$ are much tighter than large-scale cosmological
constraints $\omega\gtrsim 10^{3}$, but these decouple if the true theory of
gravity features screening. On the largest cosmological scales BD approximates
the most general second order scalar-tensor (Horndeski) theory, so constraints
here have wider implications. These will improve with upcoming large-scale
structure and CMB surveys. To constrain BD with weak gravitational lensing,
one needs its non-linear matter power spectrum $P_{\mathrm{BD}}$. By comparing
the boost $B=P_{\mathrm{BD}}/P_{\mathrm{GR}}$ from linear theory and non-
linear $N$-body simulations, we show that the non-linear boost can simply be
predicted from linear theory if the $\mathrm{BD}$ and $\mathrm{GR}$ universes
are parametrized in a way that makes their early cosmological evolution and
quasi-linear power today similar. In particular, they need the same
$H_{0}/\sqrt{\smash[b]{G_{\rm eff}(a=0)}}$ and $\sigma_{8}$, where $G_{\rm
eff}$ are their (effective) gravitational strengths. Our prediction is $1\%$
accurate for $\omega\geq 100$, $z\leq 3$ and $k\leq 1\,h/\mathrm{Mpc}$, and
$2\%$ further up to $k\leq 5\,h/\mathrm{Mpc}$. It also holds for
$G_{\mathrm{BD}}$ that do not match Newton’s constant today, so one can study
GR with different gravitational constants $G_{\mathrm{GR}}$ by sending
$\omega\rightarrow\infty$. We provide a code that computes $B$ with the linear
Einstein-Boltzmann solver hi_class and multiplies it by the non-linear
$P_{\mathrm{GR}}$ from EuclidEmulator2 to predict $P_{\mathrm{BD}}$.

###### Key Words.:

cosmology: large-scale structure of Universe – cosmological parameters –
methods: numerical

##  1 Introduction

The Brans-Dicke (BD) theory (Brans & Dicke, 1961; Dicke, 1962) of modified
gravity (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
$\omega$ and reduces to GR as $\omega\rightarrow\infty$, so the former can be
understood as a generalization of the latter. In cosmology, replacing
$\mathrm{GR}$ by $\mathrm{BD}$ in the standard ($\mathrm{GR}$-)$\Lambda$CDM
model gives the alternative $\mathrm{BD}$-$\Lambda$CDM model, where the
effective gravitational strength evolves with time and modifies the expansion
history and growth of structure.

Experiment (dataset) |  Constraint $\downarrow$ (conf.) | Year
---|---|---
Pulsar timing (PSR J0337+1715) (Voisin et al., 2020) |  $\omega>140\,000$ ($95\%$) | 2020
Pulsar timing (PSR J0337+1715) (Archibald et al., 2018) |  $\omega>\phantom{0}73\,000$ ($95\%$) | 2018
Shapiro time delay (Cassini satellite) (Bertotti et al., 2003; Will, 2014) |  $\omega>\phantom{0}40\,000$ ($95\%$) | 2003
Pulsar timing (PSR J1738+0333) (Freire et al., 2012) |  $\omega>\phantom{0}25\,000$ ($68\%$) | 2012
CMB (Planck), BAO (BOSS, 6dFGS, SDSS) (Ooba et al., 2017) |  $\omega>\phantom{00}2\,000$ ($95\%$) | 2017
CMB (Planck), BAO (BOSS, 6dFGS, SDSS), Hubble (several) (Amirhashchi & Yadav, 2020) |  $\omega>\phantom{00}1\,560$ ($95\%$) | 2019
CMB (Planck), BAO (BOSS), SN (Pantheon), LSS (KiDS, 2dFLenS) (Joudaki et al., 2022) |  $\omega>\phantom{00}1\,540$ ($95\%$) | 2022
CMB (Planck) (Avilez & Skordis, 2014) |  $\omega>\phantom{000\,}890$ ($99\%$) | 2014
CMB (Planck), BAO (6dFGS, SDSS) (Umiltà et al., 2015) |  $\omega>\phantom{000\,}208$ ($95\%$) | 2015
CMB (WMAP, Planck), BAO (SDSS, 6dF) (Li et al., 2013) |  $\omega>\phantom{000\,}182$ ($95\%$) | 2013
CMB (WMAP, others), LSS (2dF) (Acquaviva et al., 2005) |  $\omega>\phantom{000\,}120$ ($95\%$) | 2005
CMB (WMAP, others), LSS (SDSS) (Wu et al., 2010; Wu & Chen, 2010) |  $\omega>\phantom{000\,0}98$ ($95\%$) | 2010
Table 1: Constraints on the Brans-Dicke parameter $\omega$ from various tests,
beginning with the Cassini bound from 2003. Sorted by $\omega$, the table
naturally divides into two “camps” from small-scale and large-scale
constraints. We exclude $\omega<0$ due to ghost instability. Note that Avilez
& Skordis (2014); Joudaki et al. (2022) report separate constraints with and
without fixing $G_{m0}=G$ today; we show only the tightest latter bounds. Note
also that Joudaki et al. (2022)’s bounds $\omega>\\{1540,160,160,350\\}$
depend on the prior used for $\omega$ in the Markov chains. See also Nagata et
al. (2004); Ooba et al. (2016, 2017); Ballardini et al. (2016, 2020); Umiltà
et al. (2015) for constraints on the gravitational strength and generalized
Brans-Dicke models.

As summarized in Table 1, observations in the $21$st century have placed very
strong constraints on $\omega$. Already in 2003, Shapiro time delay
measurements of the Cassini satellite on its way to Saturn constrained
$\omega\gtrsim 10^{4}$, and recent strong-field tests from timing of rapidly
rotating neutron stars (pulsars) even bounds $\omega\gtrsim 10^{5}$. On
cosmological scales, however, the story is different. Roughly speaking,
constraints have tightened from only $\omega\gtrsim 10^{2}$ with stage-II
survey data in the 2000s to $\omega\gtrsim 10^{3}$ 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 CMB experiments (Abazajian et al., 2016).
Fisher forecasts found that if $\omega=800$, then Euclid should be able to
constrain this to $3\%$ using galaxy clustering and weak lensing data
(Frusciante et al., 2023). Likewise, the forecast of Ballardini et al. (2019)
found that upcoming clustering and weak lensing data in combination with BOSS
and CMB observations has the potential to reach $\omega\gtrsim 10\,000$ 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\neq c$ from GW170817 (Ezquiaga & Zumalacárregui, 2017), and Solà Peracaula
et al. (2019); Sola et al. (2020) show that it alleviates the $H_{0}$ and
$S_{8}$ 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 like BD
on the largest scales, like 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 non-linear
$\mathrm{BD}$-$\Lambda$CDM matter power spectrum $P_{\mathrm{BD}}(k\lesssim
5\,h/\mathrm{Mpc},\,z\lesssim 3\,|\,\boldsymbol{\theta})$ across parameter
space $\boldsymbol{\theta}$ 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 non-linear evolution. For example:

  * •

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 non-linear matter distribution from
halo statistics, and can be integrated directly with linear Einstein-Boltzmann
solvers. The hmcode framework have been extended to take into account a wide
range of dark energy and modified gravity 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 both the halo model
and perturbation theory to model corrections coming from non-standard physics.
This is implemented in the publicly available ReACT code (Bose et al., 2020;
Giblin et al., 2019).

This approach was taken in Joudaki et al. (2022), who used $N$-body
simulations to create a fitting formula in hmcode to predict the non-linear
$\mathrm{BD}$-$\Lambda$CDM power-spectrum. As far as we know, this is the only
non-linear 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 non-linear matter
power spectrum for arbitrary cosmologies. There already exists sophisticated
emulators for the $\mathrm{GR}$-$\Lambda$CDM power spectrum, such as
EuclidEmulator2 (Euclid Collaboration et al., 2021), CosmicEmu (Heitmann et
al., 2016; Lawrence et al., 2017; Moran et al., 2023) and BACCO (Angulo et
al., 2021; Aricò et al., 2021a, b). Recent works like Winther et al. (2019),
Brando et al. (2022), Fiorini et al. (2023), Sesame (Mauland et al., 2023) and
e-MANTIS (Sáez-Casares et al., 2024) have produced emulators for selected
modified gravity 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 $\mathrm{BD}$-$\Lambda$CDM power spectrum
emulator exists.

  * •

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

We demonstrate a hybrid technique, somewhat reminiscent of both rescaling and
emulation, suited to extending existing predictors of
$\mathrm{GR}$-$\Lambda$CDM to $\mathrm{BD}$-$\Lambda$CDM. By carefully
selecting the cosmological parameters $\boldsymbol{\theta}_{\mathrm{BD}}$ and
$\boldsymbol{\theta}_{\mathrm{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 _non-linear_
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, as we will verify by
comparing it to the boost obtained from $N$-body simulations. In turn, we can
predict the full $\mathrm{BD}$-$\Lambda$CDM power spectrum
$P_{\mathrm{BD}}=B\cdot P_{\mathrm{GR}}$ by combining it with any existing
high-quality emulator for $P_{\mathrm{GR}}$. This “trick” exploits the effort
that has already gone into creating sophisticated $\mathrm{GR}$ emulators, and
avoids duplicating this work for every alternative model to be explored. It
could be a viable cheap route for constraining modified gravity theories with
upcoming surveys (Frusciante et al., 2023; Casas et al., 2017).

This paper is structured as follows. Section 1 is this introduction. 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 motivates the parametrization we use to
compute the boost. Section 5 shows the resulting linear and non-linear boost
and compares it 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_{\mu\nu}$, it features a
scalar field $\phi$ with a constant dimensionless parameter $\omega$. Here we
will use units where $c=1$, but explicitly keep Newton’s constant $G=6.67\cdot
10^{-11}\,\mathrm{Nm}^{2}/\mathrm{kg}^{2}$ to show how BD effectively modifies
the gravitational constant of GR. In the Jordan frame, the total action of
Brans-Dicke gravity minimally coupled to the matter action $S_{M}$ is

$S\left[g_{\mu\nu},\phi\right]=\frac{1}{16\pi
G}\int\mathrm{d}^{4}x\sqrt{-|g|}\,\bigg{[}\phi
R-\frac{\omega}{\phi}(\nabla\phi)^{2}\bigg{]}+S_{M},$ (1)

where $g_{\mu\nu}$ is the spacetime metric with determinant $|g|<0$, Ricci
scalar $R$ and covariant derivative $\nabla_{\mu}$; and
$(\nabla\phi)^{2}=(\nabla^{\mu}\phi)(\nabla_{\mu}\phi)=g^{\mu\nu}(\nabla_{\mu}\phi)(\nabla_{\nu}\phi)$.
From the principle of least variaction $\delta S=0$, the (classical) equations
of motion for $g_{\mu\nu}$ and $\phi$ are the modified Einstein and Klein-
Gordon field equations (Clifton et al., 2012)

$\displaystyle G_{\mu\nu}$ $\displaystyle=\frac{8\pi
G}{\phi}T_{\mu\nu}+\frac{\omega}{\phi^{2}}\bigg{[}(\nabla_{\mu}\phi)(\nabla_{\nu}\phi)-\frac{1}{2}g_{\mu\nu}(\nabla\phi)^{2}\bigg{]}$
$\displaystyle+\frac{1}{\phi}\bigg{[}\nabla_{\mu}\nabla_{\nu}\phi-
g_{\mu\nu}\nabla^{2}\phi\bigg{]},$ (2a) $\displaystyle\nabla^{2}\phi$
$\displaystyle=\frac{8\pi G}{3+2\omega}T,$ (2b)

where $\nabla^{2}=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the 4-dimensional
Laplacian, and $T_{\mu\nu}=(-2/\sqrt{-g})\,(\delta S_{M}/\delta g^{\mu\nu})$
is the energy-momentum tensor of matter with trace
$T=T^{\mu}_{\phantom{\mu}\mu}$. Notice that $\phi$ effectively promotes
Newton’s constant $G$ in the Einstein field equations (2a) to a dynamical
field $G/\phi$. In the following, we will see in more detail how BD is similar
to GR with such an effective gravitational strength.

GR is recovered in the limit $\omega\rightarrow\infty$: then
$\nabla^{2}\phi=0$, whose solution $\phi=1$ everywhere recovers Einstein’s
field equations $G_{\mu\nu}=8\pi G\,T_{\mu\nu}$ with Newton’s constant $G$.

###  2.1 Cosmology

We assume a spatially flat, homogeneous and isotropic Friedmann-Lemaître-
Robertson-Walker (FLRW) background with small perturbations in the Newtonian
gauge

$\mathrm{d}s^{2}=-\big{[}1+2\Psi(t,\mathbf{x})\big{]}\mathrm{d}t^{2}+a^{2}(t)\big{[}1+2\Phi(t,\mathbf{x})\big{]}\mathrm{d}\mathbf{x}^{2}.$
(3)

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

$T^{\mu\nu}=\sum_{s}T^{\mu\nu}_{(s)}=\sum_{s}(\rho_{s}+P_{s})u_{s}^{\mu}u_{s}^{\nu}+P_{s}g^{\mu\nu},$
(4)

with energy densities $\rho_{s}$, pressures $P_{s}$ and 4-velocities
$u^{\mu}_{s}$ from the species $s$:

  * •

radiation $r$ (from photons $\gamma$ and massless111Including _massive_
neutrinos is straightforward, but previous studies (such as Mauland et al.
(2023)) have shown that they have little effect on the matter power spectrum
_boost_. neutrinos $\nu$) with equation of state
$\bar{P}_{r}=\bar{\rho}_{r}/3$,

  * •

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

  * •

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

We then calculate perturbatively up to first order to study the universe’s
background and linear perturbations.

###  2.2 Background cosmology

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

$\displaystyle H^{2}$ $\displaystyle=\frac{8\pi
G}{3\bar{\phi}}\bar{\rho}+\frac{\omega}{6}\left(\frac{\dot{\bar{\phi}}}{\bar{\phi}}\right)^{2}-H\frac{\dot{\bar{\phi}}}{\bar{\phi}},\quad({}_{\mu\nu}={}_{00})$
(5a)
$\displaystyle\frac{1}{a^{3}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\dot{\bar{\phi}}a^{3}\right)$
$\displaystyle=\frac{8\pi G}{3+2\omega}(\bar{\rho}-3\bar{P}),$ (5b)

where $\bar{\rho}=\bar{\rho}_{r}+\bar{\rho}_{m}+\bar{\rho}_{\Lambda}$,
$\bar{P}=\bar{P}_{r}+\bar{P}_{m}+\bar{P}_{\Lambda}$,
$\dot{{}}=\mathrm{d}/\mathrm{d}t$ and $H=\dot{a}/a$ is the Hubble parameter.
Notice that the Friedmann equation in BD is the same as in GR, only with
additional effective energy from $\smash{\dot{\bar{\phi}}}$ and an effective
gravitational strength $G_{H}=G/\bar{\phi}$.

According to the Boltzmann equation, the three species evolve like
$\bar{\rho}_{r}=\bar{\rho}_{r0}a^{-4}$, $\bar{\rho}_{m}=\bar{\rho}_{m0}a^{-3}$
and $\bar{\rho}_{\Lambda}=\bar{\rho}_{\Lambda 0}a^{0}$, where 0 indexes
quantities by their values today. We define the standard density parameters
$\Omega_{s}=\bar{\rho}_{s}/(3H^{2}/8\pi G)$, such that
$\omega_{s0}=\Omega_{s0}h^{2}\propto\bar{\rho}_{s0}$ measures the physical
densities today.222Some authors instead define
$\Omega_{s}=\rho_{s}/(3H^{2}/8\pi(G/\bar{\phi}))$ in terms of the effective
$G/\bar{\phi}$. This is natural in BD because their sum is
$\sum_{s}\Omega_{s}=1$ at all times, if one defines an effective
$\Omega_{\phi}$. However, we opt for the GR definition
$\Omega_{s}=\rho_{s}/(3H^{2}/8\pi G)$ to minimize confusion and maintain
consistency with hi_class’ input parameter
$\texttt{Omega\\_cdm}=\rho_{s}/(3H_{0}^{2}/8\pi G)$. Re-parametrizing the
evolution in terms of ${}^{\prime}=\mathrm{d}/\mathrm{d}\ln
a=H^{-1}\cdot\mathrm{d}/\mathrm{d}t$, the system (5) can be rewritten in the
dimensionless form

$\displaystyle E^{2}=\left(\frac{H}{H_{0}}\right)^{2}$
$\displaystyle=\frac{1}{\bar{\phi}}\cdot\frac{\displaystyle\Omega_{r0}a^{-4}+\Omega_{m0}a^{-3}+\Omega_{\Lambda
0}}{1-\frac{\omega}{6}[(\ln\bar{\phi})^{\prime}]^{2}+(\ln\bar{\phi})^{\prime}},$
(6a) $\displaystyle\big{(}\bar{\phi}^{\prime}E\,a^{3}\big{)}^{\prime}$
$\displaystyle=\frac{3}{3+2\omega}\frac{\Omega_{m0}+4\,\Omega_{\Lambda
0}\,a^{-3}}{E}.$ (6b)

This system can be integrated from an initial scale factor $a_{\mathrm{ini}}$,
given $\bar{\phi}(a_{\mathrm{ini}})$, $\bar{\phi}^{\prime}(a_{\mathrm{ini}})$,
$\Omega_{r0}$, $\Omega_{m0}$, $\Omega_{\Lambda 0}$ and $\omega$. But these are
not all free parameters, as today’s Friedmann equation $E_{0}=1$ constrains
one of them – say $\Omega_{\Lambda 0}$. To find a consistent solution, one can
iterate over $\Omega_{\Lambda 0}$ and solve the system repeatedly until this
constraint is satisfied.

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

$\smash{\dot{\bar{\phi}}}=\smash{\dot{\bar{\phi}}_{\mathrm{ini}}}=0\qquad\text{(during
radiation domination)}.$ (7)

In the Poisson equation (9) we will see that the effective gravitational
strength felt by matter today depends on $\omega$ and $\bar{\phi}_{0}$. We
therefore use the shooting method to vary $\bar{\phi}_{\mathrm{ini}}$ and
integrate the system repeatedly until we hit the desired $\bar{\phi}_{0}$.

###  2.3 Linear perturbations

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

$\displaystyle(\delta_{c}+3\Phi)^{\prime}$
$\displaystyle=\phantom{-}\frac{k}{aH}v_{c},$ (8a)
$\displaystyle(\delta_{b}+3\Phi)^{\prime}$
$\displaystyle=\phantom{-}\frac{k}{aH}v_{b},$ (8b)
$\displaystyle(\delta_{\gamma}+4\Phi)^{\prime}$
$\displaystyle=\frac{4}{3}\frac{k}{aH}v_{\gamma},$ (8c)
$\displaystyle(av_{c})^{\prime}$ $\displaystyle=-\frac{k}{H}\Psi,$ (8d)
$\displaystyle(av_{b})^{\prime}$
$\displaystyle=-\frac{k}{H}\Psi+\frac{4}{3}\frac{\bar{\rho}_{\gamma}}{\bar{\rho}_{b}}\left(v_{b}-v_{\gamma}\right)a\tau^{\prime},$
(8e) $\displaystyle(av_{\gamma})^{\prime}$
$\displaystyle=-\frac{k}{H}\Psi-\frac{k}{4H}\delta_{\gamma}+av_{\gamma}+(v_{\gamma}-v_{b})a\tau^{\prime},$
(8f) $\displaystyle
k^{2}\Big{(}\Phi+\frac{1}{2}\frac{\delta\phi}{\bar{\phi}}\Big{)}$
$\displaystyle=4\pi\frac{G}{\bar{\phi}}a^{2}\delta\rho,$ (8g) $\displaystyle
k^{2}\frac{\delta\phi}{\bar{\phi}}$
$\displaystyle=\frac{8\pi}{2\omega+3}\frac{G}{\bar{\phi}}a^{2}(\delta\rho-3\,\delta
P),$ (8h) $\displaystyle\Psi+\Phi$
$\displaystyle=-\frac{\delta\phi}{\bar{\phi}}.$ (8i)

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

$-k^{2}\Psi=4\pi\underbrace{\frac{2\omega+4}{2\omega+3}\frac{G}{\bar{\phi}}}_{G_{m}}\frac{\bar{\rho}_{m0}}{a}\delta_{m}=\frac{3}{2}\frac{G_{m}}{G}\frac{\Omega_{m0}H_{0}^{2}}{a}\delta_{m}.$
(9)

Notice that the Poisson equation in BD is the same as in GR, only with an
effective gravitational strength
$G_{m}=(2\omega+4)/(2\omega+3)(G/\bar{\phi})=(2\omega+4)/(2\omega+3)G_{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 $\bar{\phi}_{0}=(2\omega+4)/(2\omega+3)$ in “restricted models”
with the correct Newtonian limit $G_{m0}=G$. However, we will also allow for
“unrestricted models” with arbitrary $G_{m0}$ and $\bar{\phi}_{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
$\delta_{m}$, we neglect baryon-photon interactions at late times and combine
the Boltzmann equations (8a), (8b), (8d) and (8e) with the Poisson equation
(9) to get

$\displaystyle 0$ $\displaystyle=\ddot{\delta}_{m}+2H\dot{\delta}_{m}-4\pi
G_{m}\bar{\rho}_{m}\delta_{m}$
$\displaystyle=H^{2}\bigg{\\{}\delta_{m}^{\prime\prime}+\Big{[}2+(\ln
E)^{\prime}\Big{]}\delta_{m}^{\prime}-\frac{3}{2}\frac{G_{m}}{G}\Omega_{m}(a)\delta_{m}\bigg{\\}}.$
(10)

The general solution of this second order differential equation is a linear
combination $\delta_{m}(k,a)=C_{+}(k)D_{+}(a)+C_{-}(k)D_{-}(a)$ of a growing
mode $D_{+}(a)$ and a decaying mode $D_{-}(a)$. The former is called the
_scale-independent growth factor_ $D_{+}(a)$, and neglecting the decaying mode
$D_{-}(a)$ relates
$D_{+}(a_{2})/D_{+}(a_{1})=\delta_{m}(k,a_{2})/\delta_{m}(k,a_{1})$ for any
$k$.

On the other hand, deflection of light through weak lensing is sensitive to
the Weyl potential $(\Psi-\Phi)/2$. To find its source equation, we eliminate
$\delta\phi$ from equations (8g), (8h) and (8i) to get

$-k^{2}\left(\frac{\Phi-\Psi}{2}\right)=4\pi\frac{G}{\bar{\phi}}\bar{\rho}_{m0}a^{-1}\delta_{m}.$
(11)

Here we see that it is the effective $G_{\gamma}=G/\bar{\phi}$ that affects
light trajectories. In general, we encountered three effective gravitational
strengths $G_{\gamma}=G_{H}\neq G_{m}$ that affect light trajectories, the
expansion and matter clustering, respectively. Although they generally take on
two different values, we have $G_{\gamma}=G_{H}\approx G_{m}$ for (viable)
models with appreciable $\omega$.

Figure 1: Evolution of the effective gravitational parameter $G$, Hubble
function $E(a)=H(a)/H_{0}$, growth factor $D(a)$ and growth rate
$f(a)=\mathrm{d}\,{\ln D(a)}/\mathrm{d}\,{\ln a}$ from equations (5) and (2.3)
with the scale factor $a$ in the fiducial $\mathrm{BD}$-$\Lambda$CDM and
$\mathrm{GR}$-$\Lambda$CDM cosmologies in Table 2. The common time of matter-
radiation equality is marked by $\bar{\rho}_{r}=\bar{\rho}_{m}$.

The growth of structure is summarized by the matter power spectrum

$P(k,z)=P_{\mathrm{prim}}(k)\cdot|\Delta_{m}(k,z)|^{2},$ (12a) of the gauge-
invariant density perturbation $\Delta_{m}\approx\delta_{m}$ amplified from a
nearly scale-invariant ($n_{s}\lesssim 1$) primordial power spectrum
$P_{\mathrm{prim}}(k)=\frac{2\pi^{2}}{k^{3}}\cdot
A_{s}\left(\frac{k}{k_{p}}\right)^{n_{s}-1}$ (12b)

with primordial amplitude $A_{s}$ and spectral index $n_{s}$.

###  2.4 Non-linear structure formation

To perform $N$-body simulations in BD we only need the two equations: the
geodesic equation

$\displaystyle\frac{d^{2}\boldsymbol{x}}{dt^{2}}+2H\frac{d\boldsymbol{x}}{dt}=-\frac{1}{a^{2}}\boldsymbol{\nabla}\Psi,$
(13)

and the Poisson equation

$\displaystyle\boldsymbol{\nabla}^{2}\Psi=\frac{3}{2}\Omega_{m0}\frac{G_{m}}{G}H_{0}^{2}a^{-1}\delta_{m}.$
(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)/G$, making it trivial to implement in any
$N$-body code for GR.

Figure 1 compares the background evolution in two $\mathrm{BD}$-$\Lambda$CDM
and $\mathrm{GR}$-$\Lambda$CDM universes with cosmological parameters to be
introduced in Sect. 4 and Table 2.

##  3 Computational pipeline

ptpthi_class: solve background; solve perturbations. fml: solve background;
seed particles with $s$; evolve with COLA.ramses: solve background; read first
snapshot; evolve without COLA.Particle snapshots
$\boldsymbol{x}_{i}(z_{j})$Particle snapshots $\boldsymbol{x}_{i}(z_{j})$Model
parameters $\boldsymbol{\theta}$ and initial condition seed $s$Linear
$P_{\textsc{hi\\_class}}(k,z)$Quasi-linear $P_{\mathrm{COLA}}(k,z)$Non-linear
$P_{\mathrm{RAMSES}}(k,z)$ Figure 2: Our $\mathrm{BD}$\- or
$\mathrm{GR}$-specific simulation pipeline begins from a set of input
parameters $\boldsymbol{\theta}_{\mathrm{BD}}$ or
$\boldsymbol{\theta}_{\mathrm{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 backscales
$P_{\textsc{hi\\_class}}(k,z)$ from $z=0$ to $z=z_{\mathrm{init}}$, 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 non-
linear 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 non-linear power
spectrum.

To compute the matter power spectrum $P(k,z\,|\,\boldsymbol{\theta})$ in a
$\mathrm{BD}$-$\Lambda$CDM or $\mathrm{GR}$-$\Lambda$CDM universe, we join
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 $\mathrm{GR}$-$\Lambda$CDM and several modified gravity
cosmologies, including $\mathrm{BD}$-$\Lambda$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_{\textsc{hi\\_class}}(k,z)$. Bellini et al. (2018) showed that its Brans-
Dicke matter power spectrum agrees to $0.1\%$ with four other codes.

  * •

fml333The fml library is available at https://github.com/HAWinther/FML/, and
the COLA $N$-body code in the subdirectory FML/COLASolver. is a particle-mesh
N-body code that uses the COLA method (Tassev et al., 2013) in
$\mathrm{GR}$-$\Lambda$CDM and several modified gravity cosmologies, including
$\mathrm{BD}$-$\Lambda$CDM. It supersedes the code mg-picola (Winther et al.,
2017), in turn based on l-picola (Howlett et al., 2015). The code starts by
generating initial particle positions from hi_class’ linear power spectrum
today using the backscaling method (see e.g. Fidler et al. (2017)). It then
solves the Poisson equation (14) with the time-dependent $G_{m}(t)$ in
$\mathrm{BD}$ and Newton’s constant $G_{m}=G$ in $\mathrm{GR}$, Finally, it
outputs the (quasi-)non-linear matter power spectrum $P_{\mathrm{COLA}}(k,z)$
by constructing the density field from the particle position snapshots. These
results are used for our main analysis. We use a box with volume
$L^{3}=(384\,h/\mathrm{Mpc})^{3}$ and
$N_{\mathrm{part}}=N_{\mathrm{cell}}=512^{3}$ particles and cells, evolving
from redshift $z_{\mathrm{init}}=10$ till today with $N_{\mathrm{step}}=30$
time steps. These simulations sacrifice accuracy at the deepest non-linear
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 quasi-linear scales with few time steps. This is an excellent trade-off in
our scenario. To some extent, the error in $P$ at smaller scales will even
cancel after we compute the _boost_ $P_{\mathrm{BD}}/P_{\mathrm{GR}}$.

  * •

ramses (Teyssier, 2002) solves the Poisson equation with the standard $N$-body
method with adaptive mesh refinement and time-stepping in
$\mathrm{GR}$-$\Lambda$CDM and (with our modifications)
$\mathrm{BD}$-$\Lambda$CDM. We have modified its $\mathrm{GR}$ version444This
ramses patch is available at https://github.com/HAWinther/JordanBransDicke. 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’ output. It starts
from the same initial particle positions at the same redshift
$z_{\mathrm{init}}$ as fml, uses the same box size $L$ and number of particles
$N_{\mathrm{part}}$, adaptively refines a base grid with the same size
$N_{\mathrm{cell}}$, but uses adaptive time stepping that ignores fml’s
$N_{\mathrm{step}}$. The snapshots are analyzed in the same way as fml, except
that we quadruple the number of grid cells along each dimension due to the
finer resolution. This also outputs the non-linear matter power spectrum
$P_{\textsc{ramses}}(k,z)$, and is only used to validate fml’s results.

All codes branch to solve the equations specific to $\mathrm{BD}$-$\Lambda$CDM
or $\mathrm{GR}$-$\Lambda$CDM at all times. To compensate for slight code
differences, such as time stepping and splining effects, we normalize $P$ from
fml and ramses to match that from hi_class at the most linear common scale
$k=0.02\,h/\mathrm{Mpc}$.

To compute the matter power spectrum boost

$B(k_{\mathrm{BD}},z\,|\,\boldsymbol{\theta}_{\mathrm{BD}})=\frac{P_{\mathrm{BD}}(k_{\mathrm{BD}},z\,|\,\boldsymbol{\theta}_{\mathrm{BD}})}{P_{\mathrm{GR}}(k_{\mathrm{GR}},z\,|\,\boldsymbol{\theta}_{\mathrm{GR}})},$
(15)

we first run one $\mathrm{BD}$ simulation with freely chosen parameters
$\boldsymbol{\theta}_{\mathrm{BD}}$, followed by one $\mathrm{GR}$ simulation
with parameters
$\boldsymbol{\theta}_{\mathrm{GR}}=\boldsymbol{\theta}_{\mathrm{GR}}(\boldsymbol{\theta}_{\mathrm{BD}})$
decided by a transformation of the former. We then compare the power spectra
at equal redshifts $z$, but generally different wavenumbers $k_{\mathrm{BD}}$
and $k_{\mathrm{GR}}=k_{\mathrm{GR}}(k_{\mathrm{BD}})$, as we will motivate in
Sect. 4.

We minimize cosmic variance _within_ each $(\mathrm{BD},\mathrm{GR})$
simulation pair by using a common initial condition seed
$s_{\mathrm{BD}}=s_{\mathrm{GR}}=s$, but change
$s=s(\boldsymbol{\theta}_{\mathrm{BD}})$ for every new set of parameters to
avoid bias towards one particular configuration of the universe. To further
reduce cosmic variance, we also use amplitude-fixed initial conditions (Angulo
& Pontzen, 2016; Villaescusa-Navarro et al., 2018).

###  3.1 Tests and convergence analysis

hi_class and fml independently calculate $E(a)=H(a)/H_{0}$, $\Omega_{\Lambda
0}$ and $\sigma_{8}$ in $\mathrm{GR}$, and also $\bar{\phi}(a)$ and
$\bar{\phi}^{\prime}(a)$ in $\mathrm{BD}$. We check that these agree within a
small tolerance after every run.

Figure 3: Convergence of fml’s boost (15) at $z=0$ in the fiducial cosmology
in Table 2, as computational parameters are varied from the fiducial
$N_{\mathrm{part}}=512^{3}$, $N_{\mathrm{cell}}=512^{3}$,
$N_{\mathrm{step}}=30$, $z_{\mathrm{init}}=10$ and $L=384\,\mathrm{Mpc}/h$.

Figure 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. 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\leq\\{1/2,1,2\\}\cdot k_{\mathrm{Nyquist}}$
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=P_{\mathrm{BD}}/P_{\mathrm{GR}}$, which agree to $1\%$ up to $k\lesssim
5\,h/\mathrm{Mpc}$.

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

To test the absolute accuracy of the particle-mesh simulations we run using
fml we also compare 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\lesssim 5\,h/\mathrm{Mpc}$. Together, these tests show
that the results of the boost from fml are both precise and accurate to $1\%$
error. We only perform these tests for the fiducial cosmology in Table 2,
where we have chosen a small value $\omega=100$ in $\mathrm{BD}$ to maximize
the deviation from GR, and assume the results can be trusted for other
cosmological parameters.

##  4 Model parametrization

What parameter and wavenumber transformations
$\boldsymbol{\theta}_{\mathrm{GR}}(\boldsymbol{\theta}_{\mathrm{BD}})$ and
$k_{\mathrm{GR}}(k_{\mathrm{BD}})$ should we use in the boost (15)? That is,
_how_ should we compare a $\mathrm{BD}$ universe to a $\mathrm{GR}$ universe?
It may sound as obvious as the question sounds overpedantic 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 $\bar{\rho}_{r0}$,
$\bar{\rho}_{m0}$, $\bar{\rho}_{\Lambda 0}$, $H_{0}$ and $\bar{\phi}_{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_{\mathrm{GR}}$ to predict $P_{\mathrm{BD}}$, anyway.

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

(a) Equal $\sigma_{8}=0.80$ and $\phi_{\mathrm{ini}}h^{2}=0.67$

(b) Equal $A_{s}=2\cdot 10^{-9}$ and $\phi_{\mathrm{ini}}h^{2}=0.67$ (c) Equal
$\sigma_{8}=0.80$ and $h=0.7$

(d) Equal $A_{s}=2\cdot 10^{-9}$ and $h=0.7$

Figure 5: Matter power spectrum boost (15) with different parameter
transformations
$\boldsymbol{\theta}_{\mathrm{GR}}=\boldsymbol{\theta}_{\mathrm{GR}}(\boldsymbol{\theta}_{\mathrm{BD}})$,
where ($A_{s}^{\mathrm{BD}}=A_{s}^{\mathrm{GR}}$ or
$\sigma_{8}^{\mathrm{BD}}=\sigma_{8}^{\mathrm{GR}}$) and
($H_{0}^{\mathrm{BD}}=H_{0}^{\mathrm{GR}}$ or
$H_{0}^{\mathrm{GR}}=H_{0}^{\mathrm{BD}}\bar{\phi}_{\mathrm{ini}}^{1/2}$).
Requiring equal $\sigma_{8}$ (instead of equal $A_{s}$) flattens the boost in
the non-linear regime, while requiring equal $\phi_{\mathrm{ini}}h^{2}$
(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 $\tau$. 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_{H}=G/\bar{\phi}=3H^{2}/8\pi\bar{\rho}$ from the Friedmann
equation (5a) with the frozen scalar field (7). For large enough $\omega$ that
sufficiently suppress the anisotropic stress $\delta\phi/\bar{\phi}$,
perturbations in BD and GR with equal $k/H$ (or $k/\sqrt{G_{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_{H}}$). Similarly,
(Zahn & Zaldarriaga, 2003, equation (12)) shows that a mode $k^{\prime}$ in a
GR universe with gravitational constant $G^{\prime}$ evolves like a mode
$k=k^{\prime}/\sqrt{G^{\prime}/G}$ in another GR universe with Newton’s
constant $G$. As explained there: “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 is also true for BD’s scale-independent
effective gravitational strength.

This breaks when recombination is activated, 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 _matter_ power spectrum that concerns us here.

This motivates matching modes $\\{k_{\mathrm{BD}},k_{\mathrm{GR}}\\}$ such
that $k/H$ is the same at early times. To be compatible with $N$-body codes
that factor $H_{0}$ out of the relevant equations, we can compare modes with
equal $k/H_{0}$ as long as we make $E_{\mathrm{BD}}=E_{\mathrm{GR}}$ at early
times. The frozen scalar field (7) cuts the Friedmann equation (5a) off at
$H_{\mathrm{BD}}^{2}=H_{\mathrm{GR}}^{2}/\bar{\phi}$ provided we use the same
physical densities $\bar{\rho}_{s0}\propto\omega_{s0}$ in BD and
GR,555Strictly speaking, we cannot have
$\bar{\rho}_{s0}^{\mathrm{GR}}=\bar{\rho}_{s0}^{\mathrm{BD}}$ for all species
$s$ because $\bar{\rho}_{\Lambda 0}$ is constrained by the Friedmann equation,
but we have $\bar{\rho}_{\mathrm{GR}}\approx\bar{\rho}_{\mathrm{BD}}$ at
_early times_ , when $\bar{\rho}_{r}\gg\bar{\rho}_{m}\gg\bar{\rho}_{\Lambda}$.
so we can accomplish this by transforming

$H_{0}^{\mathrm{GR}}=H^{\mathrm{BD}}_{0}\bar{\phi}_{\mathrm{ini}}^{1/2}\quad\iff\quad
H_{0}^{\mathrm{GR}}/\sqrt{\smash[b]{G_{\mathrm{GR}}^{\mathrm{ini}}}}=H_{0}^{\mathrm{BD}}/\sqrt{\smash[b]{G_{\mathrm{BD}}^{\mathrm{ini}}}}.$
(16)

This changes the simulation box size $L=384\,\mathrm{Mpc}/h$ with
$H_{0}=100\,h\,\mathrm{km}/\mathrm{s}/\mathrm{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}}=\bar{\phi}_{\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

$s=\int_{0}^{t}c_{s}\frac{\mathrm{d}t}{a}=\int_{0}^{a}\sqrt{\frac{4\bar{\rho}_{\gamma}/3\bar{\rho}_{m}}{1+4\bar{\rho}_{\gamma}/3\bar{\rho}_{m}}}\frac{\mathrm{d}a}{H},$
(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}\propto 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_{\mathrm{BD}}(k_{\mathrm{BD}})$ and $P_{\mathrm{GR}}(k_{\mathrm{GR}})$,
avoiding oscillations in their ratio $P_{\mathrm{BD}}/P_{\mathrm{GR}}$. Figure
5 shows the suppressed oscillations when parametrized by
$\bar{\phi}_{\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 parametrize 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 $\mathrm{BD}$ and
$\mathrm{GR}$ _today_ than at early times. The normalization of today’s power
spectrum is conventionally done in terms of the parameter

$\sigma_{R}^{2}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}P_{0}(k)\hat{W}_{R}^{2}(k)k^{2}\mathrm{d}k,$
(18)

where $\hat{W}_{R}(k)=3\,[\sin(kR)-kR\cos(kR)]/(kR)^{3}$ is the Fourier
transform of the top-hat profile $W_{R}(r)=\Theta(R-r)/(4\pi R^{3}/3)$ with
radius $R$. It is common to define the normalization using the value of
$\sigma_{8}$, i.e. $\sigma_{R}$ with $R=8\,\mathrm{Mpc}/h$. By pinning
$\sigma_{8}^{\mathrm{BD}}=\sigma_{8}^{\mathrm{GR}}$ today,666In principle, we
must shoot to hit a given $\sigma_{8}(A_{s})$. But $A_{s}$ scales the whole
linear power spectrum (12a) and is detached from the linear perturbations, and
the $\sigma_{8}$ integral (18) dies off at non-linear scales. To avoid
shooting, we can guess any value $A_{s}^{\prime}$ of $A_{s}$, find
$(\sigma_{8}^{\prime})^{2}\propto A_{s}^{\prime}$ linearly, and just proceed
with $A_{s}=(\sigma_{8}/\sigma_{8}^{\prime})^{2}A_{s}^{\prime}$. we get less
primordial power $A_{s}^{\mathrm{BD}}<A_{s}^{\mathrm{GR}}$ to compensate for
the increased structure growth in $\mathrm{BD}$ due to both less Hubble
friction $E_{\mathrm{BD}}<E_{\mathrm{GR}}$ and stronger gravity
$G_{\mathrm{BD}}>G_{\mathrm{GR}}$ (see Fig. 1). Figure 5 shows that using the
same $A_{s}$ gives a peak in the boost, while using the same $\sigma_{8}$
gives a flatter boost.

This shows that the parametrization in Fig. 5(a) with $\sigma_{8}$ and
$\bar{\phi}_{\mathrm{ini}}h^{2}$ gives the simplest and most predictable non-
linear boost. As an added bonus, it also deviates the least from its linear
counterpart, and we will soon see that this holds through parameter space.

Table 2: Independent cosmological parameters $\boldsymbol{\theta}$ for $\mathrm{BD}$-$\Lambda$CDM and $\mathrm{GR}$-$\Lambda$CDM and their fiducial values. In $\mathrm{GR}$, $\omega=\infty$, $G_{m0}/G=1$ and $\bar{\phi}_{\mathrm{ini}}=1$ implicitly. Listed parameters are varied, while the photon temperature today $T_{\gamma 0}=2.7255\,\mathrm{K}$, effective massless neutrino number $N_{\mathrm{eff}}=3.0$ and primordial power spectrum pivot scale $k_{\mathrm{pivot}}=0.05\,\mathrm{Mpc}$ are fixed, and we always assume a flat universe. Derived parameters like the early scalar field $\bar{\phi}_{\mathrm{ini}}$, Hubble constant $h$, primordial power spectrum amplitude $A_{s}$ and dark energy density parameter $\omega_{\Lambda 0}$ can differ in BD and GR universes with equal independent parameters. Parameter | Fiducial value
---|---
$\omega$ | $100$
$G_{m0}/G=\bar{\phi}^{-1}(2\omega+4)/(2\omega+3)$ | $1$
$h\,\bar{\phi}_{\mathrm{ini}}^{1/2}$ | $0.674$
$\sigma_{8}$ | $0.80$
$\omega_{m0}=\Omega_{m0}h^{2}$ | $0.15$
$\omega_{b0}=\Omega_{b0}h^{2}$ | $0.02$
$n_{s}$ | $1.0$

##  5 Resulting matter power spectrum boost

We now focus on the boost with the particular parametrization developed in
Sect. 4, defaulting to the fiducial parameters in Table 2. In particular, due
to the cosmological constraints $\omega\gtrsim 1000$ summarized in Table 1, we
restrict ourselves to $\omega\geq 100$ and use the minimum $\omega=100$ in the
fiducial cosmology.

###  5.1 Understanding the linear boost evolution

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 (2.3), giving

$B(k,z)=B_{\mathrm{prim}}(k)\cdot\left(\frac{D_{+}^{\mathrm{BD}}(a(z))}{D_{+}^{\mathrm{GR}}(a(z))}\right)^{2},$
(19)

where $B_{\mathrm{prim}}$ is the boost of the primordial spectra (12b). Figure
6 shows that the two are indeed roughly consistent:

Figure 6: Evolution of the linear matter power spectrum boost from hi_class
compared to the scale-independent growth (19) obtained from Fig. 1. The boosts
are divided by that of the primordial spectra (12b) and calculated in the
fiducial cosmology in Table 2. Notice that the curves for $z=0$ and $z=9$
overlap, and that the boost is highest at $z=1.5$ due to the peak in Fig. 1.

  1. 1.

At the earliest times, structure has not had time to form, so $P\rightarrow
P_{\mathrm{prim}}$ and $B\rightarrow B_{\mathrm{prim}}$.

  2. 2.

During radiation domination, the frozen scalar field (7) and Hubble parameter
transformation (16) make $E_{\mathrm{BD}}=E_{\mathrm{GR}}$, as seen in Fig. 1.
This creates equal Hubble frictions $(\ln E_{\mathrm{BD}})^{\prime}=(\ln
E_{\mathrm{GR}})^{\prime}$ in the growth equation (2.3). Only stronger gravity
$G_{\mathrm{BD}}>G_{\mathrm{GR}}$ leads to slightly faster growth in
$\mathrm{BD}$ in a nevertheless stagnant period.

  3. 3.

During matter domination, Fig. 1 shows that the scalar field in $\mathrm{BD}$
begins to move, resulting in both less Hubble friction
$E_{\mathrm{BD}}^{\prime}<E_{\mathrm{GR}}^{\prime}$ and stronger gravity
$G_{\mathrm{BD}}>G_{\mathrm{GR}}$ at the same time. This leads to faster
growth in $\mathrm{BD}$!

  4. 4.

During dark energy domination, Fig. 1 shows that the expansion accelerates
more in $\mathrm{BD}$ than in $\mathrm{GR}$ with
$E_{\mathrm{BD}}^{\prime}>E_{\mathrm{GR}}^{\prime}$, while the gravitational
strengths $G_{\mathrm{BD}}\rightarrow G_{\mathrm{GR}}$ also become equal. This
causes the boost to peak around $z=1.5$ and then fall off.

###  5.2 Cosmological dependence of the non-linear boost

Figure 7: Variation of the matter power spectrum boost (15) at $z=0$ as each
of the cosmological parameters in Table 2 are varied individually away from
their fiducial values. The bottom right plot is special; showing the ratio of
the non-linear 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\approx B_{\mathrm{lin}}$ to $1\%$ for $k\leq 1\,h/\mathrm{Mpc}$ and $2\%$
for $k\leq 5\,h/\mathrm{Mpc}$.

We now study how the non-linear 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}=\bar{\phi}^{3/2}$) very resistant
towards 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\leq 1\,h/\mathrm{Mpc}$ and $2\%$ up to $k\leq
5\,h/\mathrm{Mpc}$.

###  5.3 Linear prediction of the non-linear boost

Conveniently, this lets us get away with predicting the _non-linear_ 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
non-linear power spectrum $P_{\mathrm{BD}}=B_{\mathrm{lin}}\cdot
P_{\mathrm{GR}}$ from an existing non-linear prediction of $P_{\mathrm{GR}}$.
In fact, existing $P_{\mathrm{GR}}$ emulators like EuclidEmulator2 only output
the non-linear 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}}(\boldsymbol{\theta}_{\mathrm{BD}})=C(\boldsymbol{\theta}_{\mathrm{GR}}(\boldsymbol{\theta}_{\mathrm{BD}}))\cdot
P_{\mathrm{BD}}^{\mathrm{lin}}(\boldsymbol{\theta}_{\mathrm{BD}})$. However,
we still need to calculate $P_{\mathrm{GR}}^{\mathrm{lin}}$ in order to
calculate the $\sigma_{8}$ integral (18) from the perturbations to make the
parameter transformation
$\boldsymbol{\theta}_{\mathrm{GR}}(\boldsymbol{\theta}_{\mathrm{BD}})$.

We provide the program bd.py777The bd.py code is available at
https://github.com/hersle/jbd. that handles these subtleties and predicts the
non-linear $\mathrm{BD}$ spectrum from $2$ runs of hi_class for the linear
spectra $\mathrm{BD}$ and $\mathrm{GR}$ spectra and $1$ run of EuclidEmulator2
for the non-linear $\mathrm{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 non-linear prediction

To demonstrate our program bd.py, we compare it to the prediction of Joudaki
et al. (2022) that we mentioned in Sect. 1. They incorporate $\mathrm{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\lesssim
1\,h/\mathrm{Mpc}$ and $5$-$10\%$ up to $k\lesssim 10\,h/\mathrm{Mpc}$. Figure
8 shows that their prediction was (almost) within their tolerance, but
produces too little power in $\mathrm{BD}$ when trialed with a tighter $1\%$
tolerance.

Figure 8: A non-linear $\mathrm{BD}$-$\Lambda$CDM power spectrum (boost) with
$\omega=1000$ 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).

##  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 non-linear BD matter power spectrum from a set of COLA
$N$-body simulations. Instead, we found that it is possible to predict the
non-linear boost $B=P_{\mathrm{BD}}/P_{\mathrm{GR}}$ using linear theory and
codes _if_ one makes a particular transformation
$\boldsymbol{\theta}_{\mathrm{BD}}\rightarrow\boldsymbol{\theta}_{\mathrm{GR}}$
of the cosmological parameters. In detail, by comparing the linear boost with
the non-linear boosts from both COLA and standard $N$-body simulations, we
verified that

$\displaystyle\frac{P_{\mathrm{BD}}(k/h_{\mathrm{BD}},z)}{P_{\mathrm{GR}}(k/h_{\mathrm{GR}},z)}\overset{2\%}{\simeq}\frac{P_{\mathrm{BD}}^{\mathrm{lin}}(k/h_{\mathrm{BD}},z)}{P_{\mathrm{GR}}^{\mathrm{lin}}(k/h_{\mathrm{GR}},z)}$
$\displaystyle\text{for}\quad\begin{Bmatrix}[l]k\lesssim 5\,h/\mathrm{Mpc}\\\
z\lesssim
3\end{Bmatrix}\quad\text{if}\quad\begin{Bmatrix}[l]\omega_{\mathrm{BD}}\gtrsim
100\\\ \sigma_{8}^{\mathrm{GR}}=\sigma_{8}^{\mathrm{BD}}\\\
H_{0}^{\mathrm{GR}}=H_{0}^{\mathrm{BD}}\smash{\bar{\phi}_{\mathrm{ini}}^{1/2}}\\\
\bar{\rho}_{m0}^{\mathrm{GR}}=\bar{\rho}_{m0}^{\mathrm{BD}}\\\
\bar{\rho}_{r0}^{\mathrm{GR}}=\bar{\rho}_{r0}^{\mathrm{BD}}\end{Bmatrix}.$

This “trick” exploits the boost’s dependence on cosmological parameters to map
the non-linear 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 non-linear $N$-body methods. When paired
with an existing non-linear $P_{\mathrm{GR}}$ predictor, such as
EuclidEmulator2, this paves a simple and efficient path for predicting
$P_{\mathrm{BD}}$ and constrain Brans-Dicke theory with the precision of
upcoming stage-IV large-scale structure surveys.We provide the program bd.py
for this at https://github.com/hersle/jbd.

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

We also encourage others to investigate whether such simplifying parameter
mappings
$\boldsymbol{\theta}_{\mathrm{MG}}\rightarrow\boldsymbol{\theta}_{\mathrm{GR}}$
exist for other modified gravity ($\mathrm{MG}$) theories. Our findings could
be specific to Brans-Dicke gravity, or apply to more theories, perhaps with a
similar scale-independent nature. Even if one does not find
$P_{\mathrm{MG}}(\boldsymbol{\theta}_{\mathrm{MG}})/P_{\mathrm{GR}}(\boldsymbol{\theta}_{\mathrm{GR}})\approx
P_{\mathrm{MG}}^{\mathrm{lin}}(\boldsymbol{\theta}_{\mathrm{MG}})/P_{\mathrm{GR}}^{\mathrm{lin}}(\boldsymbol{\theta}_{\mathrm{GR}})$
“exactly”, investing thought in a clever parameter transformation may still
pay off by simplifying the shape of $P_{\mathrm{MG}}/P_{\mathrm{GR}}$, so 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 $\Lambda$CDM with modified gravity, so one
does not have to duplicate the effort gone into general relativity for every
such alternative.

###### Acknowledgements.

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

## References

  * Abazajian et al. (2016) Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, arXiv e-prints, arXiv:1610.02743 [arXiv:1610.02743]
  * Acquaviva et al. (2005) Acquaviva, V., Baccigalupi, C., Leach, S. M., Liddle, A. R., & Perrotta, F. 2005, Phys. Rev. D, 71, 104025 [arXiv:astro-ph/0412052]
  * Alonso et al. (2017) Alonso, D., Bellini, E., Ferreira, P. G., & Zumalacárregui, M. 2017, Phys. Rev. D, 95, 063502 [arXiv:1610.09290]
  * Amirhashchi & Yadav (2020) Amirhashchi, H. & Yadav, A. K. 2020, Physics of the Dark Universe, 30, 100711 [arXiv:1908.04735]
  * Angulo & Pontzen (2016) Angulo, R. E. & Pontzen, A. 2016, Monthly Notices of the Royal Astronomical Society, 462, L1 [arXiv:1603.05253]
  * Angulo & White (2010) Angulo, R. E. & White, S. D. M. 2010, MNRAS, 405, 143 [arXiv:0912.4277]
  * Angulo et al. (2021) Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 5869 [arXiv:2004.06245]
  * Archibald et al. (2018) Archibald, A. M., Gusinskaia, N. V., Hessels, J. W. T., et al. 2018, Nature, 559, 73 [arXiv:1807.02059]
  * Aricò et al. (2021a) Aricò, G., Angulo, R. E., Contreras, S., et al. 2021a, Monthly Notices of the Royal Astronomical Society, 506, 4070 [arXiv:2011.15018]
  * Aricò et al. (2021b) Aricò, G., Angulo, R. E., & Zennaro, M. 2021b, arXiv e-prints, arXiv:2104.14568 [arXiv:2104.14568]
  * Avilez & Skordis (2014) Avilez, A. & Skordis, C. 2014, Phys. Rev. Lett., 113, 011101 [arXiv:1303.4330]
  * Ballardini et al. (2020) Ballardini, M., Braglia, M., Finelli, F., et al. 2020, J. Cosmology Astropart. Phys., 2020, 044 [arXiv:2004.14349]
  * Ballardini et al. (2022) Ballardini, M., Finelli, F., & Sapone, D. 2022, Journal of Cosmology and Astroparticle Physics, 2022, 004
  * Ballardini et al. (2016) Ballardini, M., Finelli, F., Umiltà, C., & Paoletti, D. 2016, J. Cosmology Astropart. Phys., 2016, 067 [arXiv:1601.03387]
  * Ballardini et al. (2019) Ballardini, M., Sapone, D., Umiltà, C., Finelli, F., & Paoletti, D. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 049–049
  * Bellini et al. (2018) Bellini, E., Barreira, A., Frusciante, N., et al. 2018, Phys. Rev. D, 97, 023520 [arXiv:1709.09135]
  * Bellini et al. (2020) Bellini, E., Sawicki, I., & Zumalacárregui, M. 2020, J. Cosmology Astropart. Phys., 2020, 008 [arXiv:1909.01828]
  * Bertotti et al. (2003) Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374
  * Blas et al. (2011) Blas, D., Lesgourgues, J., & Tram, T. 2011, J. Cosmology Astropart. Phys., 2011, 034 [arXiv:1104.2933]
  * Bose et al. (2020) Bose, B., Cataneo, M., Tröster, T., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 4650 [arXiv:2005.12184]
  * Brando et al. (2022) Brando, G., Fiorini, B., Koyama, K., & Winther, H. A. 2022, J. Cosmology Astropart. Phys., 2022, 051 [arXiv:2203.11120]
  * Brans & Dicke (1961) Brans, C. & Dicke, R. H. 1961, Physical Review, 124, 925
  * Casas et al. (2017) Casas, S., Kunz, M., Martinelli, M., & Pettorino, V. 2017, Physics of the Dark Universe, 18, 73 [arXiv:1703.01271]
  * Cataneo et al. (2019) Cataneo, M., Lombriser, L., Heymans, C., et al. 2019, MNRAS, 488, 2121 [arXiv:1812.05594]
  * Clifton et al. (2012) Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep, 513, 1 [arXiv:1106.2476]
  * DESI Collaboration (2016) DESI Collaboration. 2016, arXiv e-prints, arXiv:1611.00036 [arXiv:1611.00036]
  * Dicke (1962) Dicke, R. H. 1962, Physical Review, 125, 2163
  * Euclid Collaboration et al. (2021) Euclid Collaboration, Knabenhans, M., Stadel, J., et al. 2021, MNRAS, 505, 2840 [arXiv:2010.11288]
  * Ezquiaga & Zumalacárregui (2017) Ezquiaga, J. M. & Zumalacárregui, M. 2017, Phys. Rev. Lett., 119, 251304 [arXiv:1710.05901]
  * Fidler et al. (2017) Fidler, C., Tram, T., Rampf, C., et al. 2017, J. Cosmology Astropart. Phys., 2017, 043 [arXiv:1702.03221]
  * Fiorini et al. (2023) Fiorini, B., Koyama, K., & Baker, T. 2023, J. Cosmology Astropart. Phys., 2023, 045 [arXiv:2310.05786]
  * Freire et al. (2012) Freire, P. C. C., Wex, N., Esposito-Farèse, G., et al. 2012, MNRAS, 423, 3328 [arXiv:1205.1450]
  * Frusciante et al. (2023) Frusciante, N., Pace, F., Cardone, V. F., et al. 2023, arXiv e-prints, arXiv:2306.12368 [arXiv:2306.12368]
  * Giblin et al. (2019) Giblin, B., Cataneo, M., Moews, B., & Heymans, C. 2019, Monthly Notices of the Royal Astronomical Society, 490, 4826–4840
  * Heitmann et al. (2016) Heitmann, K., Bingham, D., Lawrence, E., et al. 2016, ApJ, 820, 108 [arXiv:1508.02654]
  * Horndeski (1974) Horndeski, G. W. 1974, International Journal of Theoretical Physics, 10, 363
  * Howlett et al. (2015) Howlett, C., Manera, M., & Percival, W. J. 2015, Astronomy and Computing, 12, 109 [arXiv:1506.03737]
  * Joudaki et al. (2022) Joudaki, S., Ferreira, P. G., Lima, N. A., & Winther, H. A. 2022, Phys. Rev. D, 105, 043522 [arXiv:2010.15278]
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193 [arXiv:1110.3193]
  * Lawrence et al. (2017) Lawrence, E., Heitmann, K., Kwan, J., et al. 2017, ApJ, 847, 50 [arXiv:1705.03388]
  * Li et al. (2013) Li, Y.-C., Wu, F.-Q., & Chen, X. 2013, Phys. Rev. D, 88, 084053 [arXiv:1305.0055]
  * LSST Science Collaboration (2009) LSST Science Collaboration. 2009, arXiv e-prints, arXiv:0912.0201 [arXiv:0912.0201]
  * Mauland et al. (2023) Mauland, R., Winther, H. A., & Ruan, C.-Z. 2023, arXiv e-prints, arXiv:2309.13295 [arXiv:2309.13295]
  * Mead et al. (2021) Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401 [arXiv:2009.01858]
  * Mead et al. (2016) Mead, A. J., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 [arXiv:1602.02154]
  * Mead et al. (2015) Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958 [arXiv:1505.07833]
  * Moran et al. (2023) Moran, K. R., Heitmann, K., Lawrence, E., et al. 2023, MNRAS, 520, 3443 [arXiv:2207.12345]
  * Nagata et al. (2004) Nagata, R., Chiba, T., & Sugiyama, N. 2004, Phys. Rev. D, 69, 083512 [arXiv:astro-ph/0311274]
  * Ooba et al. (2016) Ooba, J., Ichiki, K., Chiba, T., & Sugiyama, N. 2016, Phys. Rev. D, 93, 122002 [arXiv:1602.00809]
  * Ooba et al. (2017) Ooba, J., Ichiki, K., Chiba, T., & Sugiyama, N. 2017, Progress of Theoretical and Experimental Physics, 2017, 043E03 [arXiv:1702.00742]
  * Orjuela-Quintana & Nesseris (2023) Orjuela-Quintana, J. B. & Nesseris, S. 2023, J. Cosmology Astropart. Phys., 2023, 019 [arXiv:2303.14251]
  * Perivolaropoulos & Skara (2022) Perivolaropoulos, L. & Skara, F. 2022, New A Rev., 95, 101659 [arXiv:2105.05208]
  * Ruan et al. (2024) Ruan, C.-Z., Cuesta-Lazaro, C., Eggemeier, A., et al. 2024, MNRAS, 527, 2490 [arXiv:2301.02970]
  * Sáez-Casares et al. (2024) Sáez-Casares, I., Rasera, Y., & Li, B. 2024, MNRAS, 527, 7242 [arXiv:2303.08899]
  * Smith et al. (2003) Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [arXiv:astro-ph/0207664]
  * Sola et al. (2020) Sola, J., Gomez-Valent, A., de Cruz Perez, J., & Moreno-Pulido, C. 2020, arXiv e-prints, arXiv:2006.04273 [arXiv:2006.04273]
  * Solà Peracaula et al. (2019) Solà Peracaula, J., Gómez-Valent, A., de Cruz Pérez, J., & Moreno-Pulido, C. 2019, ApJ, 886, L6 [arXiv:1909.02554]
  * Takahashi et al. (2012) Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [arXiv:1208.2701]
  * Tassev et al. (2013) Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, J. Cosmology Astropart. Phys., 2013, 036 [arXiv:1301.0322]
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337 [arXiv:astro-ph/0111367]
  * Umiltà et al. (2015) Umiltà, C., Ballardini, M., Finelli, F., & Paoletti, D. 2015, J. Cosmology Astropart. Phys., 2015, 017 [arXiv:1507.00718]
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Naess, S., Genel, S., et al. 2018, The Astrophysical Journal, 867, 137
  * Voisin et al. (2020) Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020, A&A, 638, A24 [arXiv:2005.01388]
  * Will (2014) Will, C. M. 2014, Living Reviews in Relativity, 17, 4 [arXiv:1403.7377]
  * Winther et al. (2019) Winther, H. A., Casas, S., Baldi, M., et al. 2019, Phys. Rev. D, 100, 123540 [arXiv:1903.08798]
  * Winther et al. (2017) Winther, H. A., Koyama, K., Manera, M., Wright, B. S., & Zhao, G.-B. 2017, J. Cosmology Astropart. Phys., 2017, 006 [arXiv:1703.00879]
  * Wu & Chen (2010) Wu, F.-Q. & Chen, X. 2010, Phys. Rev. D, 82, 083003 [arXiv:0903.0385]
  * Wu et al. (2010) Wu, F.-Q., Qiang, L.-E., Wang, X., & Chen, X. 2010, Phys. Rev. D, 82, 083002 [arXiv:0903.0384]
  * Yahya et al. (2015) Yahya, S., Bull, P., Santos, M. G., et al. 2015, MNRAS, 450, 2251 [arXiv:1412.4700]
  * Zahn & Zaldarriaga (2003) Zahn, O. & Zaldarriaga, M. 2003, Phys. Rev. D, 67, 063002 [arXiv:astro-ph/0212360]
  * Zumalacárregui et al. (2017) Zumalacárregui, M., Bellini, E., Sawicki, I., Lesgourgues, J., & Ferreira, P. G. 2017, J. Cosmology Astropart. Phys., 2017, 019 [arXiv:1605.06102]


# Qualitative Lipschitz to bi-Lipschitz decomposition

David Bate<EMAIL_ADDRESS>Zeeman Building, University of Warwick,
Coventry CV4 7AL. ORCiD: 0000-0003-0808-2453

###### Abstract.

We prove that any Lipschitz map that satisfies a condition inspired by the
work of David may be decomposed into countably many bi-Lipschitz pieces.

This work was supported by the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. 948021)

Fix complete metric spaces $(X,d)$ and $(Y,\rho)$ and $s>0$. A map $f\colon
X\to Y$ is _Lipschitz_ if there exists $L\geq 0$ such that
$\rho(f(x),f(y))\leq Ld(x,y)$ for all $x,y\in X$. An injective Lipschitz map
is _bi-Lipschitz_ if $f^{-1}$ is Lipschitz. Let $\mathcal{H}^{s}$ denote the
$s$-dimensional Hausdorff (outer) measure on a metric space and, for $S$ a
subset of a metric space, let

$\Theta^{*,s}(S,x):=\limsup_{r\to 0}\frac{\mathcal{H}^{s}((B(x,r)\cap
S))}{(2r)^{s}}\quad\text{and}\quad\Theta^{s}_{*}(S,x):=\liminf_{r\to
0}\frac{\mathcal{H}^{s}((B(x,r)\cap S))}{(2r)^{s}}.$

######  Theorem.

Let $f\colon X\to Y$ be Lipschitz, $\mathcal{H}^{s}(X)<\infty$ and, for
$\mathcal{H}^{s}$-a.e. $x\in X$, suppose that

(1) $\limsup_{r\to 0}\frac{\mathcal{H}^{s}(B(f(x),\lambda_{x}r)\setminus
f(B(x,r)))}{(2\lambda_{x}r)^{s}}<\frac{1}{2}\Theta_{*}^{s}(Y,f(x))$

for some $0<\lambda_{x}\leq 1$. Then there exists a countable Borel
decomposition $X=N\cup\bigcup_{i}X_{i}$ with $\mathcal{H}^{s}(N)=0$ such that
each $f|_{X_{i}}$ is bi-Lipschitz.

This theorem appeared in Bate and Li [3, Theorem 1.2, iii) $\Rightarrow$ (R)]
as an intermediate step in order to prove characterisations of rectifiable
subsets of a metric space. These characterisations initiated a line of
research of geometric measure theory in metric spaces and form a cornerstone
of the results in [1, 2]. The Theorem is a qualitative analogue of results of
David [4] and Semmes [7] and the proof in [3] adapts the 85 page argument in
[7] to the qualitative setting. We give a new self contained proof.

A form of the following condition is present in [7, 4, 3].

######  Definition.

For $0<\kappa<\lambda\leq 1$ and $0<\xi<1$, a function $f\colon V\subset X\to
Y$ satisfies the condition $D(\lambda,\kappa,\xi)$ on a set $S\subset V$ if,
for all $r<\operatorname{diam}S$,

$\mathcal{H}^{s}(B(f(x),\lambda r)\setminus f(V\cap
B(x,r)))<\frac{1}{2}\xi\mathcal{H}^{s}(B(f(x),\kappa r)).$

A key idea of David [4] uses $D(\lambda,\kappa,\xi)$ to deduce bi-Lipschitz
bounds on $f$.

######  Lemma.

Let $0<\kappa<\lambda\leq 1$, $0<\xi<1$ and suppose $f\colon V\subset X\to Y$
satisfies $D(\lambda,\kappa,\xi)$ on $S\subset V$. Let $x,y\in S$, set
$r=d(x,y)/4$ and suppose $\mathcal{H}^{s}(B(f(x),\kappa
r))\geq\mathcal{H}^{s}(B(f(y),\kappa r))$. If

(2) $\mathcal{H}^{s}(f(V\cap B(x,r))\cap f(V\cap
B(y,r)))\leq(1-\xi)\mathcal{H}^{s}\left(B\left(f(x),\kappa r\right)\right)$

then

(3) $\rho(f(x),f(y))\geq(\lambda-\kappa)\frac{d(x,y)}{4}.$

_Proof._ Suppose that (3) does not hold. Then by the triangle inequality,

$B\left(f(x),\lambda r\right)\cap B\left(f(y),\lambda r\right)\supset
B\left(f(x),\kappa r\right).$

Combining this with $D(\lambda,\kappa,\xi)$ negates (2). Indeed, it gives

$\displaystyle\mathcal{H}^{s}(f(V\cap B(x,r))\cap f(V\cap B(y,r)))$
$\displaystyle>\mathcal{H}^{s}\left(B\left(f(x),\lambda r\right)\cap
B\left(f(y),\lambda r\right)\right)$
$\displaystyle\quad-\xi\mathcal{H}^{s}\left(B\left(f(x),\kappa
r\right)\right)$
$\displaystyle\geq(1-\xi)\mathcal{H}^{s}\left(B\left(f(x),\kappa
r\right)\right).$

###### Proof of the Theorem.

First let $V\subset X$ be Borel and note that (1) holds for the function
$f|_{V}$ and for $\mathcal{H}^{s}$-a.e. $x\in V$. Indeed, for
$\mathcal{H}^{s}$-a.e. $x\in V$, $\Theta^{*,s}(X\setminus V,x)=0$ (see [5,
§2.10.18]) and for such an $x$,

$\displaystyle\limsup_{r\to
0}\frac{\mathcal{H}^{s}(B(f(x),\lambda_{x}r)\setminus f(V\cap
B(x,r)))}{(2\lambda_{x}r)^{s}}$ $\displaystyle\leq\limsup_{r\to
0}\frac{\mathcal{H}^{s}(B(f(x),\lambda_{x}r)\setminus
f(B(x,r)))}{(2\lambda_{x}r)^{s}}$ $\displaystyle\quad+\limsup_{r\to
0}\frac{L^{s}\mathcal{H}^{s}(B(x,r)\setminus V)}{(2\lambda_{x}r)^{s}}$ (4)
$\displaystyle<\frac{1}{2}\Theta_{*}^{s}(Y,f(x))+\frac{L^{s}}{\lambda_{x}^{s}}\Theta^{*,s}(X\setminus
V,x).$

The first inequality uses $\mathcal{H}^{s}(f(A))\leq L^{s}\mathcal{H}^{s}(A)$
for any $A\subset X$, see [5, Corollary 2.10.11]. A consequence of (Proof of
the Theorem.) is

(5) $\mathcal{H}^{s}(V)>0\Rightarrow\mathcal{H}^{s}(f(V))>0.$

Now let $V_{1}\subset X$ be Borel with $\mathcal{H}^{s}(V_{1})>0$. By the
coarea formula (see [5, §2.10.25]), for $\mathcal{H}^{s}$-a.e. $y\in Y$,
$\operatorname{card}f^{-1}(y)<\infty$. Hence, by (5), for
$\mathcal{H}^{s}$-a.e. $x\in V_{1}$,

(6) $\operatorname{card}\\{x^{\prime}\in V_{1}:f(x^{\prime})=f(x)\\}<\infty.$

Let $V_{2}\subset V_{1}$ be a positive measure Borel set for which (6) holds
for all $x\in V_{2}$. By the Lusin-Novikov theorem (see [6, Exercise 18.14]
111In fact, since $f$ is continuous, (6) holds and $V_{2}$ may be chosen
compact, the result is elementary. ), there exists a Borel function $g\colon
f(V_{2})\to V_{2}$ such that $V_{3}:=g(f(V_{2}))$ is Borel and $f(g(y))=y$ for
all $y\in f(V_{2})$. By (5),
$\mathcal{H}^{s}(f(V_{3}))=\mathcal{H}^{s}(f(V_{2}))>0$ and hence, since $f$
is Lipschitz, $\mathcal{H}^{s}(V_{3})>0$.

Note that if $x\in X$ satisfies (1), then it also satisfies (1) for all
$0<\lambda\leq\lambda_{x}$. For $i\in\mathbb{N}$ let $1\geq\lambda_{i}\searrow
0$ and define $S_{i}$ to be the set of $x\in V_{3}$ for which

$\sup_{0<r<\lambda_{i}}\frac{\mathcal{H}^{s}(B(f(x),\lambda_{i}r)\setminus
f(V_{3}\cap
B(x,r)))}{(2\lambda_{i}r)^{s}}<\inf_{0<r^{\prime}<\lambda_{i}}\frac{1}{2}\left(1-\lambda_{i}\right)^{s}\frac{\mathcal{H}^{s}(B(f(x),r^{\prime}))}{(2r^{\prime})^{s}}.$

Then, by (Proof of the Theorem.), the $S_{i}$ monotonically increase to a full
measure subset of $V_{3}$. Therefore, there exist $i\in\mathbb{N}$ and
$S^{\prime}\subset S_{i}$ with $\mathcal{H}^{s}(S^{\prime})>0$ and
$\operatorname{diam}S^{\prime}\leq\lambda_{i}$. For any $0<r<\lambda_{i}$,
setting $r^{\prime}=(1-\lambda_{i}^{2})\lambda_{i}r$ shows that $f|_{V_{3}}$
satisfies $D(\lambda_{i},(1-\lambda_{i}^{2})\lambda_{i},(1+\lambda_{i})^{-s})$
on $S^{\prime}$. Since $f|_{V_{3}}$ is injective, (2) holds for all $x,y\in
S^{\prime}$ and hence the Lemma implies that $f|_{S^{\prime}}$ is bi-
Lipschitz.

The bi-Lipschitz condition extends to the closure of $S^{\prime}$. Hence
$\overline{S^{\prime}}\cap V_{1}$ is a Borel subset of $V_{1}$ of positive
measure on which $f$ is bi-Lipschitz. Since $V_{1}$ is an arbitrary Borel
subset of positive measure, the conclusion follows by exhaustion (e.g. as in
[5, §3.2.14]). ∎

## References

  * [1] David Bate  “Characterising rectifiable metric spaces using tangent spaces”  In _Invent. Math._ 230.3, 2022, pp. 995–1070  DOI: 10.1007/s00222-022-01136-7
  * [2] David Bate  “Purely unrectifiable metric spaces and perturbations of Lipschitz functions”  In _Acta Math._ 224.1, 2020, pp. 1–65  DOI: 10.4310/acta.2020.v224.n1.a1
  * [3] David Bate and Sean Li  “Characterizations of rectifiable metric measure spaces”  In _Ann. Sci. Éc. Norm. Supér. (4)_ 50.1, 2017, pp. 1–37  DOI: 10.24033/asens.2314
  * [4] Guy David  “Morceaux de graphes lipschitziens et intégrales singulières sur une surface”  In _Rev. Mat. Iberoamericana_ 4.1, 1988, pp. 73–114  DOI: 10.4171/RMI/64
  * [5] Herbert Federer  “Geometric measure theory”  In _Classics in Mathematics_ Springer Berlin Heidelberg, 1996, pp. xiv + 676  DOI: 10.1007/978-3-642-62010-2
  * [6] Alexander S. Kechris  “Classical descriptive set theory” 156, Graduate Texts in Mathematics  Springer-Verlag, New York, 1995, pp. xviii+402  DOI: 10.1007/978-1-4612-4190-4
  * [7] Stephen Semmes  “Measure-preserving quality within mappings”  In _Rev. Mat. Iberoamericana_ 16.2, 2000, pp. 363–458  DOI: 10.4171/RMI/279