# Non-linear density-velocity dynamics in $f(R)$ gravity from spherical
collapse.

Sharvari Nadkarni-Ghosh 1 and Sandip Chowdhury1
1Department of Physics, I.I.T. Kanpur, Kanpur, U.P. 208016 India
<EMAIL_ADDRESS><EMAIL_ADDRESS>(Corresponding
<EMAIL_ADDRESS>
###### Abstract

We investigate the joint density-velocity evolution in $f(R)$ gravity using
smooth, compensated spherical top-hats as a proxy for the non-linear regime.
Using the Hu-Sawicki model as a working example, we solve the coupled
continuity, Euler and Einstein equations using an iterative hybrid Lagrangian-
Eulerian scheme. The novel aspect of this scheme is that the metric potentials
are solved for analytically in the Eulerian frame. The evolution is assumed to
follow GR at very early epochs and switches to $f(R)$ at a pre-determined
epoch. Choosing the ‘switching epoch’ too early is computationally expensive
because of high frequency oscillations; choosing it too late potentially
destroys consistency with $\Lambda$CDM. To make an informed choice, we perform
an eigenvalue analysis of the background model which gives a ballpark estimate
of the magnitude of oscillations. There are two length scales in the problem:
the comoving Compton wavelength of the associated scalar field and the width
of the top-hat. The evolution is determined by their ratio. When the ratio is
large, the evolution is scale-independent and the density-velocity divergence
relation (DVDR) is unique. When the ratio is small, the evolution is very
close to GR, except for the formation of a spike near the top-hat edge, a
feature which has been noted in earlier literature. We are able to
qualitatively explain this feature in terms of the analytic solution for the
metric potentials, in the absence of the chameleon mechanism. In the
intermediate regime, the evolution is profile-dependent and no unique DVDR
exists.

###### keywords:

cosmology: theory - large-scale structure of Universe - dark energy

##  1 Introduction

Observations of type IA supernovae have established that our Universe is
undergoing a phase of accelerated expansion. The best explanation so far, is
that this acceleration is due to an additional constant energy density
component - called the cosmological constant ($\Lambda$). On galactic scales,
observations of galaxy rotation curves have established the need for a non-
baryonic matter component. In the standard $\Lambda$CDM model, our Universe
today is made up of about 68% of $\Lambda$, 28% of dark matter and 4% baryonic
or visible matter (Planck Collaboration et al., 2020). While this model agrees
remarkably well with a host of other observational constraints, there are
still some glaring issues. The physical origin of $\Lambda$ and dark matter
are still unknown. On the observational side, recent measurements indicate
that CMB estimates of the Hubble constant ($H_{0}$) and the amplitude of
fluctuations ($\sigma_{8}$) are in disagreement with other probes of the same
quantities (for e.g., Macaulay et al. 2013; Di Valentino et al. 2021). Thus,
the need arises to explore models beyond $\Lambda$CDM. There are two broad
avenues that have been considered. One path continues to assume that the
observed acceleration is due to an additional ‘dark energy’ component and a
plethora of phenomenological models have been explored. The other path assumes
that there is no additional energy density, but instead postulates that
Einstein’s General Relativity (GR) is incorrect and needs modification (see
reviews by Sotiriou & Faraoni 2010; De Felice & Tsujikawa 2010; Clifton et al.
2012; Joyce et al. 2015, 2016; Nojiri et al. 2017 for various aspects of this
subject). Current and future cosmological surveys aim to constrain modified
gravity (MG hereafter) models through a host of observables. The number counts
of clusters, their density and velocity profiles, the splashback radius and
weak lensing maps are some of the important cluster-scale probes that will
help constrain model parameters. Similarly, redshift-space distortions that
give estimates of galaxy peculiar velocities also provide complementary
constraints (see the recent review by Baker et al. 2021).

In standard GR, the perturbations on sub-horizon scales, in the single stream
limit, obey the coupled continuity, Euler and Poisson equations for the
fractional overdensity $\delta$, the peculiar velocity ${\bf v}$ and the
peculiar gravitational potential $\Psi$. In standard GR, in the linear regime,
the density and peculiar velocity are related through the continuity equation
as

$\Theta=-f(\Omega_{m})\delta$ (1)

where $\Theta=H^{-1}\nabla_{r}\cdot{\bf v}$ is the scaled divergence of the
peculiar velocity w.r.t. the physical coordinate ${\bf r}$. $f=d\ln\delta/d\ln
a$ is defined to be the linear growth rate and is usually expressed as
$f(\Omega_{m})=\Omega_{m}^{\gamma}$, where $\gamma$ is called the growth index
and is sensitive to the cosmological model. Indeed, one of the main aims of
surveys like DESI111https://www.desi.lbl.gov/ (e.g.,Alam et al. 2020),
SDSS222https://www.sdss.org (e.g., Alam et al. 2017) and
Euclid333http://sci.esa.int/euclid/. (e.g., Amendola et al. 2018) is to test
the above relation and constrain $f$. This relation has also been used
extensively to constrain the growth rate from local flow measurements (Nusser
et al. 2012; Nusser 2017; Lilow & Nusser 2021) or for measuring the velocity
of the local group (e.g., Nusser et al. 2014). However, this relation breaks
down on non-linear scales and this deviation can contribute to the error
budget in the determination of $f$ (Nusser et al., 2012). Velocities provide a
bias-independent measure of the underlying density field and hence the non-
linear velocity-gravity relation has been of interest since the early 1990s.
Analytic estimates were given based on perturbation theory both in Eulerian
and Lagrangian frames (for e.g., Bertschinger & Dekel 1989; Nusser et al.
1991; Bernardeau 1992; Gramann 1993b, a; Chodorowski 1997; Chodorowski & Lokas
1997; Susperregi & Buchert 1997; Chodorowski et al. 1998; Kitaura et al. 2012)
and were validated and extended by numerical simulations (for e.g, Mancinelli
& Yahil 1995; Mancinelli et al. 1993; Bernardeau & van de Weygaert 1996;
Bernardeau et al. 1999; Zaroubi et al. 1999; Kudlicki et al. 2000; Ciecielg et
al. 2003; Scoccimarro 2004; Colombi et al. 2007; Hahn et al. 2015). Fitting
forms for $\Lambda$CDM were provided based on spherical collapse (Bilicki &
Chodorowski 2008) and extended to other dark energy models and triaxial
collapse (Nadkarni-Ghosh 2013; Nadkarni-Ghosh & Singhal 2016; Mandal &
Nadkarni-Ghosh 2020).

All of the above mentioned investigations of the density-velocity divergence
relation (DVDR), also sometimes called the velocity-gravity relation have
assumed standard GR. One of the important features of standard GR is that the
‘Newtonian’ gravitational potential $\Psi$, which dictates particle dynamics
and appears in the Euler equation is the same as the potential $\Phi$ that
governs the spatial curvature and is related to the density through the
Poisson equation. This equality is violated in most MG models and an extra
equation is necessary to evolve the two potentials simultaneously. This
changes the fundamental structure of the evolution equations and the interplay
between the continuity and Euler equations which couples the density and
velocity is different for MG as compared to GR (see figures 1 of Uzan 2007 and
Song & Doré 2009 for nice schematics).

In this paper, we aim to understand the joint density-velocity dynamics using
the spherical collapse model as a proxy for the non-linear regime. In the
linear regime, the DVDR arises from imposing the ‘slaving condition’. It
implies that the velocity field is proportional to the acceleration field and
is such that there are no perturbations at the big bang epoch. In first order
Eulerian perturbation theory, this is equivalent to ignoring the decaying mode
in the linear solution for $\delta$ (Peebles, 1980) and in Lagrangian
perturbation theory, it is embedded in the Zeldovich approximation
(Zel’Dovich, 1970; Buchert, 1992). In Nadkarni-Ghosh (2013), it was shown that
it is possible to extend this same idea to the non-linear regime to obtain the
non-linear DVDR. In a spherically symmetric system, the condition of no
perturbations at the big bang, sets a specific relation between the non-linear
$\delta$ and $\Theta$ at any epoch. This relation traces out a special curve
in the two dimensional $\delta-\Theta$ phase space, which was termed as the
‘Zeldovich curve’ and it was shown that this curve is an invariant set of the
dynamical system given by the joint continuity and Euler equations, restricted
to spherical symmetry. This means that the curve depends only on the
cosmological parameters that govern evolution and not on initial conditions.
Perturbations starting anywhere in phasespace asymptote to the Zeldovich curve
as they evolve. Nadkarni-Ghosh (2013) also showed that this invariance can be
potentially exploited to break parameter degeneracy between parameters that
govern initial conditions such as $\sigma_{8}$ or $n_{s}$ and those that
govern evolution such as $\Omega_{m}$ or the equation of state $w$. In this
paper, we evolve smoothed, compensated spherical top-hat overdensities and
examine the resulting density and velocity profiles using a similar phase
space description. Is there a unique relation between the density and velocity
divergence akin to GR that is invariant of the initial conditions ? For
simplicity, we focus on the $f(R)$ model and in particular consider the form
given by Hu & Sawicki (2007).

Spherical collapse (SC) in the context of MG has been considered by various
authors. Schäfer & Koyama (2008) used the Birkoff’s theorem to write equations
of motions for the outer radius of a top-hat for a phenomenological model that
interpolated between $\Lambda$CDM and the DGP model. They argued that although
Birkoff’s theorem is known to be violated in most modified gravity models (Dai
et al., 2008), the error associated with this approximation was small. Schmidt
et al. (2010) coupled SC with the halo model to compute the halo mass
function, halo bias and the non-linear matter power spectrum for the DGP
theory. Chakrabarti & Banerjee (2016) investigated collapse in $f(R)$ theories
with a perfect fluid to understand the nature of the singularity. In standard
GR, SC is most commonly used to compute the critical overdensity for collapse
$\delta_{c}$, an important ingredient of excursion set based approaches such
as the Press-Schecter formalism and its extensions (Press & Schechter, 1974;
Sheth, Mo, & Tormen, 2001; Paranjape, Sheth, & Desjacques, 2013). The same is
true in MG models. Herrera et al. (2017) computed $\delta_{c}$ for $f(R)$,
whereas Li and collaborators gave excursion estimates of the mass function in
Galileon and chameleon models and also investigated the effect of the
environment (Li & Lam, 2012; Li & Efstathiou, 2012; Barreira et al., 2013).
Martino et al. (2009) computed cluster counts in MG for a Yukawa-like
potential. Very recently, SC has also been combined with Large Deviations
Theory to make analytic predictions of the matter density probability
distribution function (Cataneo et al., 2021). Many of the above investigations
model SC in modified gravity similar to that in standard GR, but with a
modified effective Newton’s constant.

The treatment of spherical collapse in this paper is akin to that of Borisov
et al. (2012) and Kopp et al. (2013) who solve the coupled continuity, Euler,
Poisson and scalar field equations iteratively for an $f(R)$ model. Borisov et
al. (2012) solve it in the Lagrangian frame by computing the acceleration of
spherical shells at each step, whereas, Kopp et al. (2013) solve the system
directly in the Eulerian frame as a coupled non-linear PDE using realistic
profiles obtained from peaks theory (Bardeen et al., 1986). However, our
treatment differs from these papers in some important aspects. Firstly, we
assume that the Compton wavelength is purely time-dependent and does not
depend on the density of the perturbation. With this approximation, we cannot
model the chameleon screening mechanism wherein the Compton wavelength changes
with the density of the perturbation (Khoury & Weltman, 2004). However, we
illustrate the effect of screening by considering different scales for the
width of the top-hat and demonstrate the scale-dependent nature of non-linear
evolution. Secondly, we assume that the perturbations obey GR at early epochs
and switch to the $f(R)$ evolution at a relatively late time denoted as
$a_{switch}$. This is necessary because the $f(R)$ model exhibits high
frequency oscillations at early epochs (Song et al., 2007; Hu & Sawicki, 2007)
making the equations numerically stiff and computationally expensive. However,
starting the evolution too late, destroys consistency with the background
$\Lambda$CDM evolution. Thus, the switching epoch, needs to be chosen
judiciously. In this paper, we perform an eigenvalue analysis of the
background equations to estimate the magnitude of oscillations and make an
informed choice of $a_{switch}$. Such an analysis has not been presented in
the literature thus far. Finally, we do not solve the equations for the scalar
field, but instead, recast the new degree of freedom in terms of the variable
$\chi$ which is proportional to the difference of the two perturbation
potentials $\Phi$ and $\Psi$. We present a closed form analytic solution for
$\chi$ in spherical symmetry. We solve the system using a hybrid Lagrangian-
Eulerian iterative scheme, wherein, the spatial equations for the two
potentials are solved analytically in the Eulerian frame and the continuity
and Euler equation are solved in the Lagrangian frame by recasting them as
second order evolution equations for the radii of a series of concentric
shells.

The rest of the paper is organized as follows. In §2 we define the Hu-Sawicki
system and perform the eigenvalue analysis of the same. §3 lays down the set
of equations and initial conditions governing the non-linear evolution. §4
outlines the iterative method of solution and gives analytic expressions for
the potential $\chi$. In order to illustrate screening, we plot the static
solutions to the potentials as a function of the radial coordinate for
different choices of the perturbation scale. §5 evolves the compensated top-
hat profile in the linear regime and §6 extends the evolution to the non-
linear regime. The late time density and velocity fields are analysed by
plotting them on the density-velocity divergence phase space and fitting
formulae are provided in certain limits. We summarize in §7. The appendix §A
gives the details of the initial top-hat setup and §B gives the details of the
error tests performed to validate the numerical code.

##  2 Background cosmology in the Hu-Sawicki model

###  2.1 Equations of motion

We consider a $f(R)$ model of modified gravity for which the action has the
form

$S=\frac{1}{2\kappa^{2}}\int d^{4}x\sqrt{-g}\\{R+f(R)\\}+S_{m},$ (2)

where $\kappa^{2}=8\pi G$ and $S_{m}$ is a minimally coupled matter action.
The field equations obtained by varying the action with respect to the metric
$g_{\mu\nu}$ are

$(1+f_{R})R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}(R+f)+(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu})f_{R}=\kappa^{2}T_{\mu\nu},$
(3)

where $T_{\mu\nu}$ is the stress-energy tensor, $R_{\mu\nu}$ and $R$ are the
metric-dependent Ricci tensor and scalar and $f_{R}=df/dR$. For the spatially
flat FRW background metric given by
$ds^{2}=-dt^{2}+a(t)^{2}(dx^{2}+dy^{2}+dz^{2})$,

$R=12H^{2}+6HH^{\prime},$ (4)

where $H$ is the Hubble parameter and prime denotes derivative w.r.t. $\ln a$.
Assuming that dark matter is well described by a perfect fluid with zero
pressure, the stress-energy tensor is

$T_{\mu\nu}=\rho U_{\mu}U_{\nu},$ (5)

where $U^{\mu}$ is the rest-frame four-velocity, $\rho$ is the energy density.
The Friedman equation for $f(R)$ models is then given by

$H^{2}-f_{R}(HH^{\prime}+H^{2})+\frac{1}{6}f+H^{2}f_{RR}R^{\prime}=\frac{\kappa^{2}{\bar{\rho}}}{3}.$
(6)

In standard GR, the Friedmann equation is a second order equation for the
scale factor $a$. In $f(R)$ gravity, replacing $R$ using eq. (4) recasts the
Friedmann equation as a fourth order equation in time for the two unknown
functions $f$ and $a$. There are two approaches to solve eq. (6). One can
assume a particular functional form for $f(R)$ and eq. (6) reduces to a second
order equation for the scale factor akin to standard GR. Examples include the
Starobinsky and Hu-Sawicki models (Starobinsky, 1980; Hu & Sawicki, 2007).
Alternatively, one can adopt the ‘designer approach’ wherein one demands that
the background evolution be identical to the $\Lambda$CDM evolution and solves
the resulting equations for $f$ as a function of time (Song et al., 2007;
Pogosian & Silvestri, 2008). In this paper, we consider the former approach
and in particular the Hu-Sawicki model as a specific example for which $f(R)$
is given by

$f(R)=-m^{2}\frac{c_{1}(R/m^{2})^{n}}{c_{2}(R/m^{2})^{n}+1}.$ (7)

$m^{2}=H_{0}^{2}\Omega_{m,0}$ denotes the mass scale, where $H_{0}$ and
$\Omega_{m,0}$ are the Hubble constant and matter density parameters today.
There are three parameters in this model: $c_{1},c_{2}$ and $n$. Requiring
that at high redshifts (large $R$) the $f(R)$ model should reduce to standard
GR gives,

$\frac{c_{1}}{c_{2}}\approx 6\frac{\Omega_{\Lambda,0}}{\Omega_{m,0}}.$ (8)

The parameters $n$ and $c_{1}/c_{2}^{2}$ are related to the value of the
derivative $f_{R0}$ through

$f_{R0}\approx-n\frac{c_{1}}{c_{2}^{2}}\left(\frac{12}{\Omega_{m,0}}-9\right)^{-n-1}.$
(9)

Thus, choosing $n$ and $f_{R0}$ completely specifies the model. Throughout
this paper, we choose $n=1$. Any theory of modified gravity needs to satisfy
certain viability conditions (Pogosian & Silvestri, 2008). They are:
$f_{RR}>0$ to ensure a stable high-curvature regime, $1+f_{R}>0$ to have a
positive effective Newton’s constant, $f_{R}<0$ and must asymptote to zero for
large $R$ so as to ensure GR at early epochs and finally, $f_{R0}$ should be
small enough to satisfy solar system and galactic constraints. The last
condition imposes the bound: $|f_{R0}|\lesssim 10^{-6}$ (Hu & Sawicki, 2007).

We use the value $f_{R0}=-10^{-6}$ for most of this paper. However, in §2.2
and §2.3 we use $f_{R0}=-0.01$ because the new aspects of the analysis are
best illustrated with a larger $f_{R0}$. For the $\Lambda$CDM model,

$\displaystyle H_{\Lambda CDM}^{2}(a)$ $\displaystyle=$ $\displaystyle
H_{0}^{2}\left(\frac{\Omega_{m,0}}{a^{3}}+\Omega_{\Lambda,0}\right)$ (10)
$\displaystyle{\rm and}\;\;R_{\Lambda CDM}(a)$ $\displaystyle=$ $\displaystyle
3H_{0}^{2}\left(\frac{\Omega_{m,0}}{a^{3}}+4\Omega_{\Lambda,0}\right).$ (11)

###  2.2 Consistency with $\Lambda$CDM

The explicit time dependence of $f(a)$ can be obtained by substituting eq.
(11) in eq. (7). With this substitution, eqs. (4) and (6) become a two
dimensional coupled system for the variables $R$ and $H$. An important, but
not obvious, assumption in $f(R)$ models is that the solutions for $R$ and $H$
obtained from solving eqs. (4) and (6) are ‘consistent’ with the $\Lambda$CDM
values, namely, the evolution given by eqs. (10) and (11). We explicitly
perform this ‘consistency check’ to validate this assumption. Following Hu &
Sawicki (2007), we define the variables

$y_{H}=\frac{H^{2}}{m^{2}}-a^{-3}\;\;\;{\rm
and}\;\;\;y_{R}=\frac{R}{m^{2}}-3a^{-3}.$ (12)

These variables have the advantage that their values are constant for the
$\Lambda$CDM model given by $y_{H,\Lambda
CDM}=\Omega_{\Lambda,0}/\Omega_{m,0}$ and $y_{R,\Lambda
CDM}=12\Omega_{\Lambda,0}/\Omega_{m,0}$. In terms of these variables, eqs. (4)
and (6) become

$\displaystyle y_{H}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{3}y_{R}-4y_{H},$ (13) $\displaystyle y_{R}^{\prime}$
$\displaystyle=$
$\displaystyle\frac{9}{a^{3}}-\frac{1}{{\tilde{f}}_{{\tilde{R}}{\tilde{R}}}(y_{H}+a^{-3})}\left[y_{H}-{\tilde{f}}_{\tilde{R}}\left(\frac{y_{R}}{6}-y_{H}-\frac{1}{2a^{3}}\right)+\frac{\tilde{f}}{6}\right],$
(14)

where

${\tilde{f}}=\frac{f}{m^{2}}\;\;\;{\rm and}\;\;\;{\tilde{R}}=\frac{R}{m^{2}}.$
(15)

The choice of the initial epoch is arbitrary and GR is assumed to be valid at
that time. The effective equation of state 444$w_{eff}$ is defined through the
relation
$H^{2}=H_{0}^{2}\left(\frac{\Omega_{m,0}}{a^{3}}+\Omega_{de,0}\exp\left[{-\int_{a}^{1}\left\\{1+w_{eff}(u)\right\\}\frac{du}{u}}\right]\right)$,
$u$ being a dummy variable.  in terms of these variables is

$1+w_{eff}(z)=-\frac{1}{3}\frac{y_{H}^{\prime}}{y_{H}}.$ (16)

Figure 1:  The deviation of the effective equation of state from its
$\Lambda$CDM value as defined in eq. (16) for $f_{R0}=-0.01$ (left panel) and
$f_{R0}=-10^{-6}$ (right panel) plotted for three different starting epochs:
$z_{init}=4$ (blue), 9 (red) and 19 (black). The solutions are oscillatory.
The amplitude of deviation is high at earlier epochs and decreases as the
evolution proceeds. Thus, consistency with $\Lambda$CDM at low redshifts i.e.,
in the ‘dark energy era’ is better if the evolution is started earlier.
However, the oscillation frequency is higher when the evolution is started
earlier making the system numerically ‘stiff’. The right panel shows that this
deviation is significantly smaller for smaller values of $|f_{R0}|$, however
the oscillation frequency is higher (see figure 2).

Consistency with $\Lambda$CDM demands that $w_{eff}$ be close to -1. We find
that $w_{eff}$ is oscillatory and in particular the deviation from the
$\Lambda$CDM value depends sensitively on the starting epoch $a_{init}$.
Figure 1 shows $w_{eff}$ for two values of $f_{R0}$ and three values of
$a_{init}:0.2,0.1$ and $0.05$ corresponding to $z=4$ (blue), 9 (red) and $19$
(black) respectively. In order to compare with Hu & Sawicki (2007), we choose
$\Omega_{m,0}=0.24$ and $\Omega_{\Lambda,0}=0.76$ for this analysis but the
discussion is not sensitive to this choice. From the left panel, it is clear
that the deviation from $\Lambda$CDM at lower redshifts is less when the
system has a higher starting redshift. However, the oscillation frequency is
higher at higher redshifts. Computationally, this means that the equations are
more ‘stiff’ requiring a higher number of steps when starting at a higher
redshift. This behaviour poses a challenge: starting evolution too early is
computationally expensive whereas starting too late implies greater deviation
from the $\Lambda$CDM equation of state. This challenge is somewhat mitigated
when smaller values of $|f_{R0}|$ are considered (right panel). The equations
remain computationally stiff at high redshifts, however, the oscillation
amplitude is significantly reduced implying an evolution very close to
$\Lambda$CDM. For $a_{init}=0.2$, the deviation of the equation of state, from
its $\Lambda$CDM value, in the dark energy dominated era can be more than 5 %
when $f_{R0}=-0.01$ and is less than 0.05 % when $f_{R0}=-10^{-6}$. Thus, for
parameter values of this model that satisfy solar system constraints, the
oscillatory behaviour can be circumvented by starting the evolution later. In
order to understand this behaviour better and make an informed choice about
the starting redshift, we examine the above system by analysing the
eigenvalues which give insights into the nature of oscillations.

###  2.3 Eigenvalue analysis of the Hu-Sawicki background cosmology

Figure 2: Instantaneous eigenvalues of the linearized matrix ${\mathcal{A}}$
given by eq. (18). The red dashed line indicates $f_{R0}=-10^{-6}$ and the
black solid line indicates $f_{R0}=-0.01$; both models have $n=1$. The
eigenvalues are complex implying oscillations. The real part (left panel) is
negative implying that the oscillation amplitude decreases as $a$ increases.
The imaginary part, which determines the oscillation frequency is large at
early epochs and decreases as a power law in $a$ and levels off as $a$
approaches 1. $a_{trans}$ denotes the epoch where the two asymptotes cross. At
the same redshift, the instantaneous oscillation frequency is larger for
smaller values of $|f_{R0}|$ but the transition value is independent of
$f_{R0}$. The initial amplitude is lower for smaller values of $|f_{R0}|$ (not
shown here), but the rate of decrease is independent of $|f_{R0}|$ as
indicated by the real part.

We recast the system given by eqs. (13) and (14) in a more general form as

$\left\\{y^{\prime}_{H},y^{\prime}_{R}\right\\}=\\{{\mathcal{A}}_{1}(y_{H},y_{R}),{\mathcal{A}}_{2}(y_{H},y_{R})\\},$
(17)

where ${\mathcal{A}}_{1}$ and ${\mathcal{A}}_{2}$ correspond to the right hand
side of eqs. (13) and (14). If ${\mathcal{A}}_{1}$ and ${\mathcal{A}}_{2}$
were linear combinations of $y_{H}$ and $y_{R}$ with constant (in time)
coefficients, then the system would be a linear, autonomous system. For such a
system the eigenvalues of the operator matrix determine the nature of the
solution. Real, negative eigenvalues indicate a stable solution, real positive
eigenvalues indicate an unstable solution that grows with time and complex
eigenvalues indicate oscillatory solutions. The rate of change of amplitude
and frequency of the oscillations are determined by the real and complex part
of the eigenvalues respectively. For a non-linear autonomous system, the
eigenvalues of the linearized operator give information about the local
behaviour of the solution. This interpretation of eigenvalues breaks down when
the system is non-autonomous i.e., when the coefficients of the dependent
variables are time varying (Slotine & Weiping, 1991; Strogatz, 1994). The
system given by eqs. (13) and (14) and is both non-linear and non-autonomous.
Thus, while the eigenvalues of the linearized system cannot give information
about the global behaviour, they can still give qualitative information about
the instantaneous solution. In particular they can indicate the presence or
absence of oscillations and provide an estimate of the instantaneous
oscillation frequency. The linearized differential operator is given by

${\mathcal{A}}_{lin}=\left(\begin{array}[]{cc}\frac{\partial{\mathcal{A}}_{1}}{\partial
y_{H}}&\frac{\partial{\mathcal{A}}_{1}}{\partial y_{R}}\\\ \\\
\frac{\partial{\mathcal{A}}_{2}}{\partial
y_{H}}&\frac{\partial{\mathcal{A}}_{2}}{\partial y_{R}},\end{array}\right)$
(18)

where

$\displaystyle\frac{\partial{\mathcal{A}}_{1}}{\partial y_{H}}$
$\displaystyle=$
$\displaystyle-4,\;\;\;\frac{\partial{\mathcal{A}}_{1}}{\partial
y_{R}}=\frac{1}{3},$ (19)
$\displaystyle\frac{\partial{\mathcal{A}}_{2}}{\partial y_{H}}$
$\displaystyle=$
$\displaystyle-\frac{1+{\tilde{f}}_{\tilde{R}}}{{\tilde{f}}_{{\tilde{R}}{\tilde{R}}}\left(a^{-3}+y_{H}\right)}+\frac{\frac{{\tilde{f}}}{6}+y_{H}-{\tilde{f}}_{\tilde{R}}\left(-\frac{1}{2a^{3}}-y_{H}+\frac{y_{R}}{6}\right)}{{\tilde{f}}_{{\tilde{R}}{\tilde{R}}}\left(a^{-3}+y_{H}\right)^{2}}\;\;\;\;\;{\rm
and}\;\;\;\;\;\frac{\partial{\mathcal{A}}_{2}}{\partial
y_{R}}=\frac{{\tilde{f}}_{\tilde{R}}}{6{\tilde{f}}_{{\tilde{R}}{\tilde{R}}}\left(a^{-3}+y_{H}\right)}.$
(20)

${\mathcal{A}}_{lin}$ can be evaluated at any time by using the values of
${\tilde{f}},{\tilde{f}}_{\tilde{R}},{\tilde{f}}_{{\tilde{R}}{\tilde{R}}},y_{H}$
and $y_{R}$ at that time. The values of $y_{H}$ and $y_{R}$ can be read off
from the numerical solution to eqs. (13) and (14). We found that using the GR
values for $y_{H}$ and $y_{R}$ did not give significantly different answers to
the eigenvalues.

Figure 2 shows the eigenvalues of ${\mathcal{A}}_{lin}$. The black solid line
and red dashed lines correspond to $f_{R0}=-0.01$ and $f_{R0}=-10^{-6}$
respectively. The eigenvalues are complex. The left panel shows the real part
of the eigenvalue, which is negative, indicating that the amplitude of the
oscillations decreases as evolution proceeds. The right panel shows the
imaginary part which indicates the instantaneous frequency of oscillations.
The figure shows two distinct regimes. At very high redshifts, the oscillation
frequency is large and decreases as a power-law in $a$: $|\rm{Im}[\nu]|\sim
a^{\alpha}$. The best fit value was found to be $\alpha\sim-3$. At low
redshifts, the instantaneous frequency is almost a constant. The epoch where
these two asymptotes cross is marked as the transition epoch $a_{trans}\sim
0.34$ corresponding to a redshift of around $z\sim 2$. For a smaller value of
$|f_{R0}|$, the magnitude of the imaginary part is higher implying higher
frequency oscillations at the same redshift. However, we find that the
transition epoch is independent of $|f_{R0}|$ and also independent of the
value of $\Omega_{m,0}$ for flat models. For $a_{init}<<a_{trans}$, the high
frequency oscillations make the system “numerically stiff” and evolving the
system is computationally expensive.

The perturbed system described in the next section also exhibits oscillations
akin to the background. In this paper, we have expressed the perturbation
equations in terms of two temporal evolution equations for the density and
velocity and two spatial equations for the metric potentials. This form of the
equations is not amenable to the eigenvalue analysis discussed above. However,
in order to overcome the issue of numerical stiffness we adopt the following
strategy. We evolve the perturbations in standard GR starting from an early
epoch $a_{init}=0.001$ to a late epoch $a_{switch}$, after which the evolution
‘switches’ over to $f(R)$ gravity. Following the insight gained from the
eigenanalysis of the background we choose $a_{switch}$ to be in the vicinity
of $a_{trans}$. In particular, we choose $a_{switch}=0.1$. It is possible to
choose $a_{switch}$ in a more quantitative manner by setting tolerances for
deviation of equation of state or numerical stiffness, but we do not follow
such a strategy here. We are interested in the qualitative differences in non-
linear growth between the GR and $f(R)$ models. Furthermore, the DVDR, if such
a unique relation exists, should depend only the evolution equations and not
on the choice of $a_{switch}$. Thus, it suffices to only consider only a
single choice for $a_{switch}$ and we do not investigate the dependence of
non-linear growth on this parameter. Clearly, choosing $a_{switch}$ earlier
(later) will increase (decrease) the effect of the modification.

High frequency oscillations in the Hu-Sawicki model have been referred to in
the original paper (Hu & Sawicki, 2007) as well as in follow up studies
(Elizalde et al., 2012). They have also been encountered in the context of
evolution of designer models (Song et al., 2007) as well as in linear
perturbations (see for example Lima & Liddle 2013; Pogosian & Silvestri 2008;
Silvestri et al. 2013). The issue of $\Lambda$CDM-consistency has been
explicitly addressed in the literature (for example Starobinsky 2007; Appleby
& Battye 2007; Ceron-Hurtado et al. 2016; Chakraborty et al. 2021). Some of
these studies have been in the context of finding stable $f(R)$ models. We
note that stability in the Hu-Sawicki model is guaranteed by the choice of
parameters. The system is oscillatory but stable since the amplitude of
oscillations decays as evolution proceeds. To the best of our knowledge, an
understanding of the oscillatory behaviour based eigenvalue analysis of this
system has not been presented so far in the literature.

###  2.4 Compton wavelength

In the Einstein frame, the $f(R)$ action can be recast as a standard Einstein-
Hilbert action with non-minimally coupled scalar field corresponding to the
additional degree of freedom represented by $f_{R}$ (Magnano & Sokołowski,
1994; Chiba, 2003; Sotiriou & Faraoni, 2010). The mass of this field is

$m^{2}_{f_{R}}=\frac{m^{2}}{3}\left(\frac{1+{\tilde{f}}_{\tilde{R}}}{\tilde{f}_{{\tilde{R}}{\tilde{R}}}}-{\tilde{R}}\right).$
(21)

The associated physical Compton wavelength $\lambda_{C}$ 555Putting in the
dimensions this gives
$\lambda_{C}=\frac{c\sqrt{3{\tilde{f}_{{\tilde{R}}{\tilde{R}}}}}}{H_{0}\sqrt{\Omega_{m,0}}}$.
is given by

$\lambda_{C}=\frac{2\pi}{m_{f_{R}}}\approx\frac{2\pi\sqrt{3{\tilde{f}_{{\tilde{R}}{\tilde{R}}}}}}{m}.$
(22)

We have used the fact that ${\tilde{f}}_{{\tilde{R}}{\tilde{R}}}<<1$ and
${\tilde{f}}_{{\tilde{R}}{\tilde{R}}}^{-1}>>{\tilde{R}}$ throughout the
evolution. The comoving Compton wavelength $x_{C}$ and the reduced comoving
Compton wavelength ${\bar{x}}_{C}$ are defined as

$x_{C}=\frac{\lambda_{C}}{a}\;\;\;{\rm
and}\;\;\;{\bar{x}}_{C}=\frac{x_{C}}{2\pi}\approx\frac{\sqrt{3{\tilde{f}_{{\tilde{R}}{\tilde{R}}}}}}{am},$
(23)

where $m^{2}=H_{0}^{2}\Omega_{m0}$. This wavelength defines the ‘range’ of the
‘fifth-force’ corresponding to the additional degree of freedom represented by
the scalar field. The effects of modified gravity manifest only on scales less
than this range. The $f(R)$ model, obeys solar system constraints through the
‘Chameleon screening’ mechanism, wherein the Compton wavelength decreases in
high density regions reducing the range of the force (Khoury & Weltman, 2004).

Figure 3 shows the relevant scales and epochs on a plot in the $k-a$ plane.
$a_{eq}$ corresponds to the epoch when the dark energy density is equal to the
matter density in the background evolution. $a_{trans}$ denotes the epoch when
frequency of oscillations in the Hu-Sawicki background starts to level off.
Starting evolution for $a<<a_{trans}$ is computationally expensive where as
$a>>a_{trans}$ implies significant deviations from the $\Lambda$CDM
background, particularly for larger values of $|f_{R0}|$. $a_{trans}$ provides
a estimate of the optimal epoch to start evolution. At early epochs, the
cosmologically relevant scales ($\sim 1h^{-1}$Mpc) are well above the Compton
wavelength and one cannot expect to see departures from GR on those scales.
Coincidentally, for $f_{R0}$ values which are allowed by solar system and
galactic constraints the Compton wavelength becomes greater than $\sim
1h^{-1}$Mpc at around $a\sim a_{trans}$.

Figure 3: Important scales in the model marked on a $k-a$ plane (left panel).
The red dotted and orange dashed lines denote the comoving Compton wavelength
$x_{C}$ of the scalaron field implied by the Hu-Sawicki model for
$f_{R0}=-0.01$ and $f_{R0}=-10^{-6}$ respectively. At any epoch, the effects
of modified gravity begin to manifest when $x\lesssim x_{C}$. The blue solid
line indicates the scale when the given mode crosses the horizon. $a_{eq}$
denotes the epoch of dark energy-matter equality and $a_{trans}$ is the epoch
when the frequency of oscillations in the Hu-Sawicki background start to level
off to a constant.

##  3 Equations for non-linear perturbations

###  3.1 Equations

Assuming only scalar perturbations of the background, the metric in conformal
Newtonian gauge is (Ma & Bertschinger, 1995)

$ds^{2}=-a^{2}(1+2\Psi)d\tau^{2}+a^{2}(1-2\Phi)(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}).$
(24)

$\Psi$ corresponds to the gravitational potential in the Newtonian limit which
determines how particles move and $\Phi$ denotes the ‘curvature’ fluctuation
which is seeded by the matter density. ${\bf x}=(x_{1},x_{2},x_{3})$ is the
comoving coordinate and $d\tau=dt/a$ is the conformal time. The Universe is
assumed to consist only of cold dark matter modelled to be a pressureless
fluid with no anisotropic stresses. Let $\delta$ denote the fractional density
and $v_{i}$ the coordinate velocity, which is equivalent to the peculiar
velocity $adx_{i}/dt$. We employ the Quasi Static Approximation (QSA) which
comprises of two independent assumptions: (1) the time derivatives of the
potentials are small compared to the spatial derivatives and (2) the length
scales are much smaller than the horizon scale, which is the usual ‘sub-
horizon’ approximation. In standard $\Lambda$CDM, these two assumptions are
identical. The potentials vary on a time scale of $1/H$ and both assumptions
correspond to $ck>>aH$, where $c$ is the speed of light. However, in modified
theories, the potentials can vary on a time scale much smaller than the Hubble
time. In particular, the potentials can be oscillatory, similar to the
behaviour seen in the background. Suppose $\Psi\sim e^{i\omega\ln a}\Psi_{0}$,
where $\omega$ is the frequency of oscillation in the variable $\ln a$, then
the first assumption reads $ck>>a\omega H$, which need not be satisfied even
if the scales are sub-horizon since $\omega$ can be large. Thus, the two
conditions need to be assumed separately. For the Hu-Sawicki model, with the
parameters that satisfy solar system constraints, these oscillations are
undetectable, justifying this assumption (Hojjati et al., 2012; Silvestri et
al., 2013). In this limit, the Einstein equations are (Oyaizu, 2008; Schmidt
et al., 2009; Borisov et al., 2012)

$\displaystyle\nabla_{x}^{2}\delta f_{R}$ $\displaystyle=$
$\displaystyle\frac{a^{2}\delta
R}{3c^{2}}-\frac{H^{2}a^{2}\Omega_{m}\delta}{c^{2}}$ (25)
$\displaystyle\nabla_{x}^{2}\Psi$ $\displaystyle=$ $\displaystyle
2H^{2}a^{2}\Omega_{m}\delta-\frac{a^{2}\delta R}{6},$ (26)

where $\delta f_{R}$ and $\delta R$ are perturbations to $f_{R}$ and the Ricci
scalar $R$ respectively. Note that $\delta R$ is an implicit non-linear
function of $\delta f(R)$. However, assuming that the change in $R$ is small,
Taylor expanding around the background gives

$\delta f_{R}\approx f_{RR}\delta R,$ (27)

where $f_{RR}$ is evaluated at the background $R$ given by eq. (11). In the
limit of $|f_{R}|<<1$ and $|f/R|<<1$, the difference between the two
potentials $\Phi$ and $\Psi$ is related to $\delta f_{R}$ through (Hu &
Sawicki (2007))

$\Phi-\Psi=c^{2}\delta f_{R}.$ (28)

This relation is also follows from the anisotropy part of the Einstein
equations in the QSA. It is further convenient to define the variables
(Pogosian & Silvestri 2008)

$\displaystyle\chi$ $\displaystyle=$ $\displaystyle\Phi-\Psi$ (29)
$\displaystyle\Phi_{+}$ $\displaystyle=$ $\displaystyle\frac{\Phi+\Psi}{2}.$
(30)

Substituting eqs. (27), (28), (29) and (30) in eqs. (25) and (26) and
combining with the conservation equations, which stay unchanged in $f(R)$,
gives the system

$\displaystyle\frac{\partial\delta}{\partial t}+\left(\frac{{\bf
v}}{a}\cdot\nabla_{x}\right)\delta$ $\displaystyle=$
$\displaystyle-\frac{(1+\delta)}{a}(\nabla_{x}\cdot{\bf v})$ (31)
$\displaystyle\frac{\partial{\bf v}}{\partial t}+\left(\frac{{\bf
v}}{a}\cdot\nabla_{x}\right){\bf v}+H{\bf v}$ $\displaystyle=$
$\displaystyle-\frac{1}{a}\nabla_{x}\Psi$ (32)
$\displaystyle\nabla_{x}^{2}\Phi_{+}$ $\displaystyle=$
$\displaystyle\frac{3}{2}H^{2}a^{2}\Omega_{m}\delta$ (33)
$\displaystyle\nabla_{x}^{2}\chi-\frac{a^{2}}{3c^{2}f_{RR}}\chi$
$\displaystyle=$ $\displaystyle-H^{2}a^{2}\Omega_{m}\delta.$ (34)

Note that $\Psi=\left(\Phi_{+}-\frac{\chi}{2}\right)$. We have retained $\Psi$
in the right hand side of eq. (32) to illustrate that the form of eqs. (31)
and (32) is unchanged from the standard GR case. $\Omega_{m}$ refers to the
time-dependent matter density parameter, which is related to its value today
as

$\Omega_{m}(a)=\frac{H_{0}^{2}\Omega_{m,0}a_{0}^{3}}{H^{2}a^{3}}.$ (35)

This is a coupled system of equations for the variables $\delta$, ${\bf v}$,
$\Phi_{+}$ and $\chi$ which are functions of space as well as time. These have
to be solved given the initial profiles for $\delta$ and ${\bf v}$ and
boundary conditions for $\Phi_{+}$ and $\chi$. Assuming no sources at infinity
(homogeneous on large enough scales) and no forces at the origin (solution is
not singular at $x=0$), gives the boundary conditions for eqs. (33) and (34)

$\Phi_{+}(x\rightarrow\infty)=0,\;\;\;\chi(x\rightarrow\infty)=0,\;\;\;\left.\frac{\partial\Phi_{+}}{\partial
x}\right|_{x=0}=0\;\;\;{\rm and}\;\;\;\left.\frac{\partial\chi}{\partial
x}\right|_{x=0}=0.$ (36)

In what follows, we set $\Omega_{m,0}=0.32$ and $\Omega_{\Lambda,0}=0.68$ in
accordance with Planck Collaboration et al.
(2020)666https://www.cosmos.esa.int/web/planck and we use the terms ‘standard
GR’ and $\Lambda$CDM interchangeably when we refer to this model. The density
and velocity profiles are known as a function of the radial coordinate at some
initial time $a_{init}$. Throughout this work, we consider compensated top-
hats as initial conditions. The details of the set-up are given in appendix A.
The perturbation is characterised by a single scale, namely, the width of the
top-hat $x_{top}$. We note that the coefficient of $\chi$ in eq. (34) can be
re-written as ${\bar{x}}_{C}^{-2}$ using eq. (23) and eq. (34) becomes

$\nabla_{x}^{2}\chi-\frac{1}{{\bar{x}_{C}}^{2}}\chi=-H^{2}a^{2}\Omega_{m}\delta.$
(37)

The solution to $\chi$ depends on the ratio777The $Q$ defined here is
different but related to the $Q$ defined in Fourier space in Pogosian &
Silvestri (2008).

$Q=\frac{{\bar{x}}_{C}}{x_{top}}.$ (38)

When $Q>>1$, eq. (37) reduces to
$\nabla^{2}\chi\approx-H^{2}a^{2}\Omega_{m}\delta$ and $\chi$ is proportional
to $\Phi_{+}$. In this limit, the gravitational potential $\Psi$ is enhanced
by a factor of $4/3$ as compared to its value in standard GR (denoted by
$\Psi_{GR}$) i.e., $\Psi=4/3\Psi_{GR}$. When $Q<<1$, the solution for $\chi$
is given by $\chi={\bar{x}}_{C}^{2}H^{2}a^{2}\Omega_{m}\delta$. The
dimensionless value of the extra field $\chi/(H^{2}x_{top}^{2})$ is suppressed
by a factor of $Q^{2}$ and the evolution is close to GR. $Q<<1$ has been
referred to as the ‘Compton condition’ in Hu & Sawicki (2007).

In this work, we assume that ${\bar{x}_{C}}$ depends only on the value of
$f_{RR}$ in the background and does not change with the perturbation. As a
result, we do not capture the ‘chameleon screening mechanism’ wherein the
Compton wavelength is small in high density regions; this is responsible for
satisfying solar system constraints (Khoury & Weltman, 2004; Hu & Sawicki,
2007). Nevertheless, we are able to illustrate ‘screening’ by considering
perturbations of differing length scales $x_{top}$ for a fixed $f_{R0}$. We
consider three values of $x_{top}=0.002,2$ and 200 $h^{-1}$Mpc, which for
$f_{R0}=-10^{-6}$ gives $Q=1217,1.21,0.012$ at $a=1$ respectively. We refer to
$Q>>1$ and $Q<<1$ as the ‘strong’ and ‘weak’ field regimes respectively. The
evolution of $Q$ for these three scales is shown in figure 4. Note that, for
$x_{top}\sim 0.002h^{-1}{\rm Mpc}$, $Q>1$ and for $x_{top}\sim 200h^{-1}$Mpc,
$Q<1$ from $a=0.1$ to today and they can be considered as being ‘strong’ and
‘weak’ perturbations over the entire domain of evolution. In what follows, we
will refer to the three perturbations with their value of $Q$ today and ignore
the temporal evolution of $Q$.

Figure 4: The parameter $Q={\bar{x}}_{C}/(x_{top})$ plotted as a function of
$a$ for the three values of top-hat widths considered in this paper. The
dotted line denotes $a_{switch}$, which is the epoch at which evolution is
assumed to switch from GR to $f(R)$. In the text, we characterize these three
fluctuations by their value of $Q$ today and ignore the time dependence.

##  4 Method of solution

###  4.1 Outline

To proceed we further restrict ourselves to a spherically symmetric
perturbation which is modelled as a series of concentric shells. Equations
(31) and (32) involve time-derivative operators for the variables $\delta$ and
${\bf v}$, whereas eqs. (33) and (34) contain purely spatial derivatives. The
gravitational potential depends on the instantaneous density making eqs. (31)
to (34) an integro-differential equation. To solve this system we employ an
iterative scheme as follows. The system is evolved from the initial epoch,
$a_{init}$, to the intermediate epoch, $a_{switch}$, in a single step using
the equations of GR. The time domain from $a_{switch}$ to $a_{final}$ is
divided into $N_{t}$ intervals. The density and velocity profiles are known at
the beginning of each interval. We first solve the spatial differential
equations eqs. (33) and (34) in Eulerian space (fixed radial coordinates) to
compute the gravitational potential corresponding to the density at the start
of the step. Assuming that this potential is constant over the small time
interval, we propagate the system using eqs. (31) and (32) to obtain the
density and velocity at the start of the next step. This temporal evolution is
carried out in Lagrangian coordinates. The advantage of using Lagrangian
coordinates is that the Lagrangian time derivative is the total time
derivative

$\frac{d}{dt}=\frac{\partial}{\partial t}+\frac{1}{a}({\bf
v}\cdot\nabla_{x}).$ (39)

This converts the non-linear coupled p.d.e. given by eqs. (31) and (32) to a
second order ordinary differential equation. We perform error tests to check
that this iterative hybrid Lagrangian-Eulerian scheme is convergent.

###  4.2 Lagrangian Coordinates

Let origin be the centre of the sphere. Define the Lagrangian coordinate of a
shell to be its initial comoving coordinate, denoted by $q$ (we henceforth
drop the vectors since the dynamics is only radial). The (Eulerian) comoving
coordinate at any later epoch is

$x\equiv x(q,a)=A(q,a)q,$ (40)

where $A(q,a)$ can be thought of as a ‘scale factor’ of the shell at $q$. The
$q$ dependence of the scale factor arises because, in general, the density is
radially dependent. By definition, for all shells $A(q,a_{init})=1$. The
peculiar velocity at any later time $t$ is

$v(q,a)=a{\dot{A}}(q,a)q=aH(a)A^{\prime}(q,a)q,$ (41)

where ‘dot’ denotes the total derivative w.r.t. time $t$ and ‘prime’ denotes
derivatives w.r.t. $\ln a$. If $\delta(q,a_{init})$ is the density at the
initial epoch, mass conservation for each shell implies that the density at
any later epoch is

$1+\delta(q,a)=\left(1+\delta(q,a_{init})\right)\frac{q^{2}dq}{x^{2}dx}=\frac{(1+\delta_{init})}{A^{3}\left|1+\frac{q}{A}\frac{dA}{dq}\right|}.$
(42)

Using the definitions in eqs. (40) and (41) and changing the derivatives to
Lagrangian coordinates, it can be easily checked that the density given by the
above expression satisfies the continuity equation given by eq. (31). The
spherically averaged density is defined as

$\Delta(x)=\frac{3}{x^{3}}\int_{0}^{x}\delta(x,a)x^{2}dx.$ (43)

It satisfies the condition

$1+\Delta(q,a)=\frac{(1+\Delta_{init})A^{3}_{init}}{A^{3}}.$ (44)

###  4.3 Equations in the hybrid Eulerian-Lagrangian scheme

Substituting eq. (39) and eq. (41) in eq. (32), imposing spherical symmetry in
eqs. (33) and (34), changing the time variable to $\ln a$ and defining
${\tilde{\Phi}}=\Phi/H^{2},{\tilde{\chi}}=\chi/H^{2}$ and
${\tilde{\Psi}}=\Psi/H^{2}$ recasts the system as

$\displaystyle
A^{\prime\prime}+\left(2-\frac{3}{2}\Omega_{m}(a)\right)A^{\prime}$
$\displaystyle=$ $\displaystyle-\frac{1}{a^{2}q}\nabla_{q}{\tilde{\Psi}}(q,a)$
(45) $\displaystyle\frac{\partial^{2}{\tilde{\Phi}_{+}}}{\partial
x^{2}}+\frac{2}{x}\frac{\partial{\tilde{\Phi}}_{+}}{\partial x}$
$\displaystyle=$ $\displaystyle\frac{3}{2}a^{2}\Omega_{m}(a)\delta(x,a).$ (46)
$\displaystyle\frac{\partial^{2}{\tilde{\chi}}}{\partial
x^{2}}+\frac{2}{x}\frac{\partial{\tilde{\chi}}}{\partial
x}-\frac{1}{{\bar{x}}_{C}^{2}}{\tilde{\chi}}$ $\displaystyle=$
$\displaystyle-a^{2}\Omega_{m}(a)\delta(x,a),$ (47)

where the ‘prime’ denotes derivative with respect to $\ln a$. The boundary
conditions for each step become

${\tilde{\Phi}_{+}}(x\rightarrow\infty)=0,\;\;\;{\tilde{\chi}}(x\rightarrow\infty)=0,\;\;\;\left.\frac{\partial{\tilde{\Phi}_{+}}}{\partial
x}\right|_{x=0}=0,\;\;\;{\rm
and}\;\;\;\left.\frac{\partial{\tilde{\chi}}}{\partial x}\right|_{x=0}=0.$
(48)

From eq. (41) and the definition of the Lagrangian coordinate, initial
conditions for $A$ are:

$A(q,a_{init})=1\;\;\;{\rm
and}\;\;\;A^{\prime}(q,a_{init})=\frac{v_{init}(q,a_{init})}{a_{init}H(a_{init})q}.$
(49)

In the GR regime, ${\tilde{\Phi}_{+}}={\tilde{\Psi}}$. Using the solution to
Poisson’s equation in terms of $\Delta$, and the equations reduce to a second
order ODE for $A$:

$A^{\prime\prime}+\left(2-\frac{3}{2}\Omega_{m}(a)\right)A^{\prime}=-\frac{\Omega_{m}(a)A}{2}\left(\frac{1+\Delta_{init}}{A^{3}}-1\right),$
(50)

where $\Delta_{init}\equiv\Delta(q,a_{init})$ is the average density defined
in eq. (43), evaluated at the initial epoch.

###  4.4 Initial conditions

We evolve the system assuming GR from $a_{init}=0.001$ to the switching epoch
$a_{switch}=0.1$. The initial profile for $\delta(x,a=0.001)$ is a smooth
compensated top-hat. Its amplitude is chosen small enough so that the profile
does not collapse in either $f(R)$ or GR but large enough to become non-
linear. The details of the profiles used in various sections of the paper are
given in appendix A. The velocity profile at $a=0.001$ is always set by
assuming the Zeldovich condition (Zel’Dovich, 1970; Brenier et al., 2003). For
a radially varying profile, this gives

$v(x,a_{ainit})=-\frac{f(\Omega_{m})Ha\Delta(x,a_{init})x}{3},$ (51)

where $\Delta(x,a)$ is defined in eq. (43) and
$f(\Omega_{m})\sim\Omega_{m}^{0.55}$ is the linear growth factor in the GR
regime. For $a_{init}=0.001$, $\Omega_{m}\approx 1$. In Lagrangian
coordinates, at the initial time, by definition, $x=q$ and the initial density
and velocity are $\Delta(x,a_{init})=\Delta(q,a_{init})$ and
$v(x,a_{init})=v(q,a_{init})$. Using eqs. (40) and (41) gives the initial
conditions

$A(q,a_{init})=1\;\;\;{\rm for}\;{\rm all}\;\;q\;\;\;\;\;{\rm
and}\;\;\;A^{\prime}(q,a_{init})=-\frac{\Delta(q,a_{init})}{3}.$ (52)

###  4.5 Algorithm

Notation: Let $a_{init}$ and $a_{final}$ denote the initial and final epochs
of interest and $a_{switch}$ denote the epoch where the evolution switches
from GR to $f(R)$. The time between $a_{switch}$, to $a_{final}$ is divided
into $N_{t}$ intervals equi-spaced in $\ln a$. Let $a_{0}=a_{init}$,
$a_{1}=a_{switch}$, $a_{N_{t}}=a_{final}$ and $a_{2},a_{3}\dots$ denote
intermediate epochs. At any epoch $a_{n}$, the data at the start of the time
interval is defined on a 1-d grid with $N_{s}$ points uniformly spaced in $\ln
x$ space. These correspond to radial shells. $x_{i}(a_{n}),i=1,2,\dots N_{s}$
denotes the Eulerian position of the $i$-th shell at the start of the $n$-th
step. At $a_{n}$, the Lagrangian coordinate of the shell $i$-th shell is
defined as $q_{n,i}=x_{i}(a_{n})$. The density and velocity are known at the
time $a_{n}$ on the uniform grid and are given by $\delta_{i}(q_{n},a_{n})$
and $v_{i}(q_{n},a_{n})$. Here we drop the subscript ‘i’ on $q_{n}$ but it is
understood that $q_{n}$ is different for each ‘i’.

At the end of the step, i.e. at $a_{n+1}$, the Eulerian position, density and
velocity of the $i$-th shell are denoted as $x_{i}(q_{n},a_{n+1})$,
$\delta_{i}(q_{n},a_{n+1})$ and $v_{i}(q_{n},a_{n+1})$. However, the Eulerian
positions $x_{i}(q_{n},a_{n+1}),i=1,2,\dots N_{s}$ are not equi-spaced. At
$a_{n+1}$, a new grid is defined in Eulerian space which is equispaced in $\ln
x$ over the range $x_{min}$ to $x_{max}$, where $x_{min}$ and $x_{max}$
correspond to the Eulerian location of the innermost and outermost shell at
$a_{n+1}$. This uniform grid is denoted as
$x_{i^{\prime}}(a_{n+1}),i^{\prime}=1,2\dots N_{s}$. The $i^{\prime}$
indicates that the shell locations are re-defined after each step. Thus,
$x_{i}(q_{n},a_{n+1})\neq x_{i^{\prime}}(a_{n+1})$. The number of shells are
the same at each step since no shells have collapsed or crossed. The
Lagrangian coordinate at $a_{n+1}$ is given as
$q_{n+1,i^{\prime}}=x_{i^{\prime}}(a_{n+1})$. The density and peculiar
velocity are obtained by interpolating $\delta_{i}(q_{n},a_{n+1})$ and
$v_{i}(q_{n},a_{n+1})$ onto the uniform grid. They are denoted as
$\delta_{i^{\prime}}(q_{n+1},a_{n+1})$ and $v_{i^{\prime}}(q_{n+1},a_{n+1})$
respectively.
Step 1. GR evolution:

  1. 1.

Evolve the system from $a_{0}=a_{init}$ to $a_{1}=a_{switch}$ using eq. (50)
and initial conditions eq. (52).

  2. 2.

Construct the Eulerian position, density and velocity at $a_{1}$ using eqs.
(41) and (42). This gives $x_{i}(q_{0},a_{1})$, $\delta_{i}(q_{0},a_{1})$ and
$v_{i}(q_{0},a_{1})$ for each shell $i$, $i=1\dots N_{s}$.

  3. 3.

Re-initialize the system by defining a uniform grid in $\ln x$ space, define
the new Lagrangian coordinate $q_{1}$ and interpolate the density and peculiar
velocity to get $\delta_{i^{\prime}}(q_{1},a_{1})$ and
$v_{i^{\prime}}(q_{1},a_{1})$.

Step 2. $f(R)$ evolution:

  1. 1.

At $a_{1}$, and at any $a_{n}$ thereafter, first solve the spatial equations,
eqs. (46) and (47), setting $a=a_{n}$ and $x=q_{n}$. The solutions for
$\Phi_{+}$ and $\chi$ can be obtained analytically and are given in the next
section. These expressions are evaluated numerically using the values of
$\delta_{i}(q_{n},a_{n})$ known at the start of the step $a=a_{n}$.

  2. 2.

Using the solutions for $\Phi_{+}$ and $\chi$ compute the r.h.s of eq. (45)
i.e., $\nabla_{q}\Psi(q,a)=\nabla_{q}\Psi(q,a_{n})$ at each shell $i$. We
assume that throughout the interval $a_{n}$ to $a_{n+1}$,
$\nabla_{q}\Psi_{i}(q_{n},a)=\nabla_{q}\Psi_{i}(q_{n},a_{n})$ i.e., the force
is constant.

  3. 3.

Compute the temporal evolution between $a_{n}$ and $a_{n+1}$ using eq. (45).
There are $N_{s}$ such equations, one for each shell. Let $A_{i}$ denote the
scale factor of the $i$-th shell. The initial conditions are

$A_{i}(q_{n},a_{n})=1\;\;\;{\rm
and}\;\;\;A_{i}^{\prime}(q_{n},a_{n})=\frac{v(q_{n},a_{n})}{a_{n}H(a_{n})q_{n}}$
(53)

  4. 4.

Construct the density and velocity at $a_{n+1}$ using eqs. (41) and (42). This
gives $x_{i}(q_{n},a_{n+1})$ and $\delta_{i}(q_{n},a_{n+1})$ and
$v_{i}(q_{n},a_{n+1})$ for the $i$-th shell. Compute the end-points $x_{min}$
and $x_{max}$ at $a_{n+1}$ to be the Eulerian positions of the innermost and
outermost shells.

  5. 5.

Re-initialize: define a new uniform grid over the range $x_{min}$ and
$x_{max}$ denoted by $x_{i^{\prime}}(a_{n+1}),i=1,2\dots N_{s}$. Define the
new Lagrangian coordinate $q_{n+1}=x_{i^{\prime}}(a_{n+1})$ and interpolate
the data given by $\delta_{i}(q_{n},a_{n+1})$ and $v_{i}(q_{n},a_{n+1})$ onto
this uniform grid to get the initial data for the next step
$\delta_{i^{\prime}}(q_{n+1},a_{n+1})$ and $v_{i^{\prime}}(q_{n+1},a_{n+1})$
at equi-spaced shell positions.

  6. 6.

Go to step 2(i) until $a_{n+1}=a_{final}$.

###  4.6 Analytic Solutions for the potentials

The solutions to eqs. (46) and (47) subject to the boundary conditions eq.
(48) are

$\displaystyle{\tilde{\Phi}_{+}}(x)$ $\displaystyle=$
$\displaystyle-\frac{\mathcal{P}_{1}}{3}\int_{x}^{\infty}y\Delta(y)dy,$ (54)
$\displaystyle{\tilde{\chi}}(x)$ $\displaystyle=$
$\displaystyle\frac{\mathcal{P}_{3}}{\sqrt{\mathcal{P}_{2}}}\frac{1}{x}\left[\exp(-\sqrt{\mathcal{P}_{2}}x)I_{1}(x)+\sinh(\sqrt{\mathcal{P}_{2}}x)I_{2}(x)\right]$
(55)

where

$\displaystyle I_{1}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{x}\left\\{\sinh(\sqrt{\mathcal{P}_{2}}y)\delta(y)y\right\\}dy,$
(56) $\displaystyle I_{2}(x)$ $\displaystyle=$
$\displaystyle\int_{x}^{\infty}\left\\{\exp(-\sqrt{\mathcal{P}_{2}}y)\delta(y)y\right\\}dy.$
(57)

$y$ is a dummy variable in the integrals and $\Delta$ represents the
spherically averaged fractional density contrast defined in eq. (43).
${\mathcal{P}_{1}}$, ${\mathcal{P}_{2}}$ and ${\mathcal{P}_{3}}$ are purely
time-dependent functions defined as

${\mathcal{P}}_{1}(a)=\frac{3\Omega_{m}(a)a^{2}}{2},\;\;\;{\mathcal{P}_{2}}(a)=\frac{1}{{\bar{x}}_{C}^{2}}\;\;\;{\rm
and}\;\;\;{\mathcal{P}_{3}}(a)=a^{2}\Omega_{m}(a).$ (58)

The ‘force’ that moves the shells is related to the derivatives of the
potentials:

$\displaystyle\nabla_{x}{\tilde{\Phi}_{+}}$ $\displaystyle=$
$\displaystyle\frac{\mathcal{P}_{1}}{3}x\Delta(x),$ (59)
$\displaystyle\nabla_{x}{\tilde{\chi}}$ $\displaystyle=$
$\displaystyle-\frac{{\tilde{\chi}}}{x}+\frac{{\mathcal{P}_{3}}}{x}\left[-\exp(-\sqrt{\mathcal{P}_{2}}x)I_{1}(x)+\cosh(\sqrt{\mathcal{P}_{2}}x)I_{2}(x)\right].$
(60)

Figure 5: The spatial solutions for dimensionless $\chi$, $\nabla\chi$, $\Phi$
and $\Psi$ at $a=1$ corresponding to a smoothed top-hat source $\delta(x)$
(top-panel). In the strong field limit $Q>>1$, the effect of the modification
is the strongest corresponding to the largest difference between $\Phi$ and
$\Psi$. In the weak field limit, $Q<<1$, the solution for $\chi$ is
proportional to $\delta$ and the ‘extra force’ which is proportional to
$\nabla\chi$ is largest at the edge of the top-hat. The amplitude of this
additional force is small compared to the potentials and they are effectively
equal as is expected in standard GR.

Figure 5 shows the solution for the dimensionless $\chi$ (second panel),
$\nabla\chi$ (third panel) and the potentials $\Phi$ and $\Psi$ (fourth panel)
as a function of $x$ for the density profile $\delta$ corresponding to a
smooth compensated top-hat (top-panel). The solution for $\chi$ is evaluated
at $a=1$ for $f_{R0}=-10^{-6}$. The ${\bar{x}}_{C}$ for this $f_{R0}$ at $a=1$
is 2.43 $h^{-1}{\rm Mpc}$. To illustrate how $\chi$ changes with the choice of
scale we pick three values for the width of the top-hat denoted by
$x_{top}=0.002,2$ and 200 $h^{-1}{\rm Mpc}$. This corresponds to
$Q=1217,1.21,0.012$ respectively. Changing the choice of scale is also
equivalent to changing the value of $f_{R0}$. The effect of the fifth force
depends on $\nabla\chi$. It is clear from the third panel, that when $Q<<1$,
the extra force $\nabla\chi$ is non-zero only at the edge of the top-hat. In
such a system, the evolved density develops a ‘spike’ at the top-hat boundary
(see §6.4). This phenomenon has been reported in the literature before
(Borisov et al., 2012; Kopp et al., 2013; Lombriser, 2016), but in the context
of the chameleon mechanism. Here, we provide a mathematical explanation of
this feature based on the analytic solution, but in the absence of chameleon
screening. A more detailed discussion can be found in §6.4. The lowest panel
shows the two potentials $\Phi$ and $\Psi$. When $Q>>1$, the two potentials
clearly differ, when $Q<<1$, the effect of the ‘extra force’ is effectively
‘screened’ and the two potentials almost coincide as is expected in standard
GR. The difference in the two potentials is related to $\chi$, which is
suppressed by a factor of $Q^{2}$ in this regime.

Equation (46) implies that the potential $\Phi_{+}$ has the same equation as
the Newtonian potential in standard GR. We denote $\Phi_{+}=\Psi_{GR}$. In the
strong field limit, $\chi\approx-2/3\Phi_{+}$ which gives $\Psi\approx
4/3\Psi_{GR}$, whereas in the weak field limit $\Psi\approx\Psi_{GR}$. Figure
6 shows the ratio $\Psi/\Psi_{GR}$, evaluated at the innermost shell, as a
function of the parameter $Q$ of the profile. The strong and weak field
asymptotic values are clearly in agreement with expectations. Since, the three
values of $Q$ above were not sufficient to characterize this plot, the static
potentials were evaluated for additional profiles with different values of
$Q$. The details can be found in appendix A.

Figure 6: The ratio $\Psi/\Psi_{GR}$ evaluated at the innermost shell for
profiles with different values of $Q$. In $f(R)$, the ‘Newtonian potential’
$\Psi$ i.e., the potential $\Psi$, that appears in Euler’s equation is
enhanced. In the strong field limit, $Q>>1$, the ratio $\Psi/\Psi_{GR}\sim
4/3$. In the weak field limit, $Q<<1$, the ratio asymptotes to unity.

##  5 Linear evolution

Figure 7: Evolution in the linear regime. The left panel shows the initial
density profile (orange, dashed), the expected final density in GR (black,
solid), the theoretically expected (cyan) and the numerically (red, dotted)
evolved final density in the $f(R)$ model at two final epochs. The cyan and
red dotted curves coincide implying the correct linear evolution for $Q>>1$
given by the scale-independent eq. (62). The right panel shows the same data
for a profile with $Q\sim 1$. The cyan and red-dotted curves deviate since the
scale-independent linear growth equation is no longer valid. The deviation
from GR is higher in the strong field regime ($Q>>1$) as compared to the
intermediate field regime ($Q\sim 1$) as expected.

In GR, a top-hat remains a top-hat until shells cross. This follows as a
consequence of the inverse-square law of Newtonian gravity. The evolution of
each shell depends only on the mass enclosed inside the shell. The position
dependence of the equations is implicit for a radially varying density
profile. This is not the case for $f(R)$ models. The force $\nabla_{x}\Psi$ in
the Euler equation eq. (45) has a more complicated dependence on $x$ and this
introduces a location-dependent evolution for each shell even for a top-hat
density distribution. In linear theory, this manifests itself as a scale-
dependent linear growth factor. However, as seen above, on scales much smaller
than the Compton wavelength, $\Psi\approx 4/3\Psi_{GR}$ and the evolution is
scale-independent. On these scales, the linear growth equation for $\delta$
which in standard GR is

${\rm
Std.GR:}\;\;\;\;\;\;\;\;\;\delta^{\prime\prime}+\left(2+\frac{H^{\prime}}{H}\right)\delta^{\prime}-\frac{3}{2}\Omega_{m}(a)\delta=0.$
(61)

gets changed in modified gravity (MG) models to

${\rm
MG(Q>>1):}\;\;\;\delta_{MG}^{\prime\prime}+\left(2+\frac{H^{\prime}}{H}\right)\delta_{MG}^{\prime}-2\Omega_{m}(a)\delta_{MG}=0.$
(62)

The solution to this system can be written as

$\delta_{GR}(a)=\delta(a_{init})\frac{D_{GR}(a)}{D_{GR}(a_{init})}=\delta(a_{init}){\tilde{D}}_{GR}(a)\;\;\;\;{\rm
and}\;\;\;\;\delta_{MG}(a)=\delta(a_{init})\frac{D_{MG}(a)}{D_{MG}(a_{init})}=\delta_{MG}(a_{init}){\tilde{D}}_{MG}(a),$
(63)

where $D_{GR}(a)$ corresponds to the growth factor given by the growing mode
solution and we define the scaled growth factor ${\tilde{D}}(a)$ through the
second equality. An analytic expression for $D(a)$ can be obtained for dark
energy models (Dodelson, 2003). To compute $D_{MG}$ for the $f(R)$ model, we
solve eqs. (61) and (62) with $\delta(a_{init})=0.001$ from $a_{init}=0.1$ to
$a_{final}=1$ with Zeldovich initial conditions
($\delta^{\prime}(a_{init})=\delta(a_{init})$). The ratio
$\delta_{MG}(a)/\delta(a_{init})$ gives ${\tilde{D}}_{MG}(a)$.

To confirm the expectations in the linear regime, we evolve a compensated top-
hat density from $a_{init}=0.1$ up to two final epochs $a_{final}=0.2$ and
$a_{final}=0.4$ and for two values of $Q$: $Q\sim 1217$ and $Q\sim 1.21$. The
$Q<<1$ is not expected to show significant deviation from GR. Since the
evolution is over a small time interval, it suffices to use a small number of
steps: $N_{t}=40$. The initial velocity is set by assuming Zeldovich initial
conditions. The initial and final densities are plotted in figure 7. The
dashed orange corresponds to the initial profile and the black line is the
theoretically expected, linearly evolved profile in GR, the red-dotted line is
the numerically evolved profile using the non-linear code for the Hu-Sawicki
model with $f_{R0}=-10^{-6}$ and the cyan line is its theoretical expectation
given by multiplying the profile with the scaled linear growth factor defined
in eq. (63). For the $Q>>1$, the cyan and the red-dotted lines coincide at
both $a=0.2$ and $a=0.4$. Their deviation from the corresponding GR value is
larger for larger final epochs. This serves as a code check of the non-linear
code in the linear regime. For $Q\sim 1$, the red-dotted lines and the cyan
lines do not coincide because eq. (62) is no longer valid. The red-dotted
lines are now closer to the GR values because the strength of the modification
reduces as $Q$ decreases. It is important to note that the linear growth given
by eq. (62) is in real space. In Fourier space, it is possible to write a
linear growth equation for each mode $\delta_{k}$ which is valid for all
scales, but is scale-dependent (see for e.g., Bean et al. 2007; Pogosian &
Silvestri 2008). Most current and upcoming large scale structure surveys aim
to constrain modified gravity parameter by measuring the linear growth rate of
perturbations. Most of these analyses are scale-independent i.e., they only
look for a deviation from the growth rate in GR (for e.g. Alam et al. 2020).
Although, a systematic scale-dependent measurement of the growth rate is
ideal, scale-independent measurements can be used to place upper limits on the
Compton wavelength of the scalaron field.

##  6 Non-linear regime

Figure 8: Initial $\delta$ at $a=0.001$ and the non-linearly evolved $\delta$,
$\Delta$, infall velocity and relative Hubble velocity ($\delta_{v}$) at $a=1$
for three values of the top-hat scale $x_{top}=0.002,2,200$ $h^{-1}{\rm Mpc}$
corresponding to $Q>>1$, $Q\sim 1$ and $Q<<1$ respectively. The red and blue
dotted curves correspond to $\sigma=0.0025$ (profile I) and $\sigma=0.1$
(profile II) evolved with the Hu-Sawicki model with $f_{R0}=-10^{-6}$. The
solid curves denote the evolution using standard GR from $a=0.001$ to $a=1$
for the same. The evolution in GR is scale independent. Thus, the final
profiles are the same for all values of $Q$. For the $f(R)$ model, When
$Q>>1$, the evolution equations are structurally similar to GR with an
enhanced potential. So the evolved profiles I and II coincide in the interior
similar to the GR case, but with an increased final density. When $Q<<1$, the
scale is outside the Compton scale making the system seemingly
indistinguishable from GR. When $Q\sim 1$, the evolution is scale-dependent
and hence profile dependent.

###  6.1 Non-linear density and velocity

We evolve the compensated top-hat profile from $a_{init}=0.001$ to
$a_{final}=1$ using the algorithm outlined in §4.5 with $a_{switch}=0.1$. We
consider three values of the top-hat scale $x_{top}=0.002,2,200$ $h^{-1}{\rm
Mpc}$ and two values of the smoothing parameter $\sigma=0.0025$ (red, dotted,
profile I) and $\sigma=0.1$ (blue, dotted, profile II), for each scale. The
details of the initial profile are given in appendix A. Figure 8 shows the
relevant quantities. The top-most panel shows the initial density profile
which is the same for both models. The density and its spherical average are
defined in eq. (42) and (43). The infall velocity and the relative Hubble
velocity are defined as

$\displaystyle v_{infall}$ $\displaystyle=$ $\displaystyle-v(q,a)$ (64)
$\displaystyle\delta_{v}$ $\displaystyle=$
$\displaystyle\frac{1}{H}\frac{\dot{r}}{r}-1=\frac{v(q,a)}{aH(a)x(q,a)}=\frac{A^{\prime}}{A}.$
(65)

The black curves denote the evolution using standard GR from $a=0.001$ to
$a=1$ and the dotted curves denote the evolution using the $f(R)$ model.

In the $\Lambda$CDM model, the radial evolution of a shell depends only on the
initial average density at the shell location and is independent of the size
of the shell (see eq. (50)). Consequently, in the GR case, the final density
profiles for the two $\sigma$ values, coincide in the interior of the top-hat
since the initial density in the interior is same for both profiles. Because
of the scale-independent evolution, the final density profiles in GR are the
same for all values of $Q$. The magnitude of the velocity depends on the scale
but the dimensionless quantity $\delta_{v}$ defined in eq. (65) is also
independent of $Q$.

In the $f(R)$ model, the evolution is, in general, scale-dependent. The
evolution of the extra potential $\chi$ depends sensitively on the parameter
$Q$. In the strong field regime when $Q>>1$, the structure of the equations is
similar to that of GR. The final densities in the interior of the top-hat are
independent of $\sigma$ as is the case in GR, but because the potential is
enhanced by a factor of $4/3$, the final density of the top-hat is higher as
compared to GR. In the intermediate regime, when $Q\sim 1$, the evolution is
profile dependent. An un-smooth top-hat is described uniquely by a single
scale $x_{top}$; smoothing broadens the top-hat edge increasing the effective
width of the top-hat and decreasing the effective $Q$. Thus, the profile I
(smaller $\sigma$) shows a greater density enhancement than profile II (larger
$\sigma$). In the weak field regime, when $Q<<1$, the extra force is suppresed
and the evolution is close to GR. There is expected to be some density
enhancement at the edge of the top-hat as was discussed in §4.6 but it is too
small to be seen on the scale of the plot. This effect is demonstrated in §6.

###  6.2 Phasespace evolution

In standard gravity, the fractional overdensity and the divergence of the
peculiar velocity are related through the coupled continuity and Euler
equations. The temporal evolution of this system requires two initial
conditions: the initial density and initial velocity field. However, it is
physically reasonable to assume that there are ‘no perturbations at the big
bang’ and this condition relates the initial density and velocity uniquely. At
first order, in Eulerian perturbation theory, this condition is implemented by
ignoring the decaying mode in the linear solution for $\delta$ and in
Lagrangian perturbation theory it is embedded in the ‘Zeldovich
approximation’. This is usually referred to as the slaving condition since the
velocity field is ‘slaved’ to the acceleration field (Brenier et al., 2003).
Thus, the linear density-velocity divergence (DVDR) relation is

$\Theta=-f\delta,$ (66)

where, $\Theta=\nabla_{r}\cdot{\bf v}/H$ is the scaled velocity divergence.
The linear growth rate, $f$, is primarily dependent on the matter density
$\Omega_{m}$ and is usually expressed as $f\equiv\Omega_{m}^{\gamma}$, where
$\gamma$, sometimes called the growth index, is a sensitive probe of
cosmology888$\gamma$ defined here differs from the post-Newtonian parameter
sometimes used in literature (Bertschinger & Zukin, 2008; Joyce et al.,
2016).. For a pure matter universe $\gamma=0.6$ (Peebles, 1980), for a
$\Lambda$CDM model $\gamma=0.55$ (Linder, 2005) and for the DGP model of
gravity, $\gamma=0.68$ (Dvali et al., 2000; Linder & Cahn, 2007).

In Nadkarni-Ghosh (2013), N13, hereafter, this relation was extended to the
non-linear regime using the spherical top-hat by imposing the condition ‘no
perturbations at the big bang’ in the full solution. At any epoch, given an
initial $\delta$, it is possible to compute the corresponding $\Theta$ that
satisfies the condition using the exact evolution equation for the outer edge
of the top-hat. This unique relationship traces out a curve in the two
dimensional $\delta-\delta_{v}$ phasespace which was called the ‘Zeldovich
curve’ because it encapsulates the spirit of the Zeldovich approximation. It
is not obvious that this is the desired late time non-linear DVDR. In order to
check this, N13, examined the dynamics of perturbations using the coupled
continuity and Euler equations, which for a top-hat in standard GR reduce to

$\displaystyle\delta^{\prime}$ $\displaystyle=$
$\displaystyle-(1+\delta)\Theta$ (67) $\displaystyle\Theta^{\prime}$
$\displaystyle=$
$\displaystyle-\frac{3}{2}\Omega_{m}(a)\delta-\Theta\left(2+\frac{H^{\prime}}{H}\right)-\frac{\Theta^{2}}{3}.$
(68)

By evolving many initial $\delta-\Theta$ pairs using the above equations, it
was shown that the ‘Zeldovich curve’ is indeed an invariant of the dynamical
system described by eqs. (67) and (68). This means that perturbations that
start anywhere in phasespace, not necessarily on the curve, eventually
converge to it. Those that start on the curve, stay on the curve. Thus, the
curve traced out by the ‘no-perturbations at the big bang’ condition indeed is
the late time, non-linear DVDR.

Finding such invariants of the dynamical system can be potentially useful in
breaking parameter degeneracies as was illustrated in N13. This is because
invariant sets depend only on cosmological parameters that appear as
coefficients in the dynamical equations and are insensitive to the parameters
that describe initial conditions, such as, the amplitude $\sigma_{8}$ or
spectral index $n_{s}$. One of the primary aims of this paper is to understand
if such invariants exist for the $f(R)$ system. Because of the complexity of
the equations, it is not practically feasible to implement the algorithm used
in N13 to select the $\delta-\Theta$ pairs that satisfy the ‘slaving
condition’. Instead, we evolve the compensated top-hat profiles which start
with Zeldovich initial conditions, compute the resulting density-velocity
divergence pairs along the evolved curves and examine where they lie in phase
space. Is the resulting curve in phase-space universal ? Or does it depend on
the details of the initial profile ?

Since the density profile is radially varying, it is easier to use the
spherically averaged $\Delta$ and the relative Hubble velocity defined in eq.
(65) 999N13, used the variables $\delta-\delta_{v}$ to describe the phase-
space. For constant density perturbations considered in that paper, the two
are equivalent $\Theta=3\delta_{v}$.. In standard GR, using relations eqs.
(44), (50) and (65) the dynamical equations for the $\Delta-\delta_{v}$ system
read:

$\displaystyle\Delta^{\prime}$ $\displaystyle=$
$\displaystyle-3(1+\Delta)\delta_{v}$ (69) $\displaystyle\delta_{v}^{\prime}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\Omega_{m}(a)\Delta-\delta_{v}\left(2+\frac{H^{\prime}}{H}\right)-\delta_{v}^{2}.$
(70)

The form of the coupled system given by eqs. (67)-(68) or eqs. (69)-(70)
hinges on the fact that the potential in Euler’s equation is related to the
density through the Poisson’s equation. For $f(R)$ theories, in general, this
closure is not possible because of the extra degree of freedom encapsulated in
$\chi$. However, in the limit $Q>>1$, $\Psi=4/3\Psi_{GR}$ and the system can
be written as

$\displaystyle\delta_{MG}^{\prime}$ $\displaystyle=$
$\displaystyle-(1+\delta_{MG})\Theta_{MG},$ (71)
$\displaystyle\Theta_{MG}^{\prime}$ $\displaystyle=$
$\displaystyle-2\Omega_{m}(a)\delta_{MG}-\Theta_{MG}\left(2+\frac{H^{\prime}}{H}\right)-\frac{\Theta_{MG}^{2}}{3},$
(72)

or alternately

$\displaystyle\Delta_{MG}^{\prime}$ $\displaystyle=$
$\displaystyle-3(1+\Delta_{MG})\delta_{v,MG},$ (73)
$\displaystyle\delta_{v,MG}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}\Omega_{m}(a)\Delta_{MG}-\delta_{v,MG}\left(2+\frac{H^{\prime}}{H}\right)-\delta_{v,MG}^{2}.$
(74)

Figure 9 shows the $\Delta-\delta_{v}$ pairs plotted in the two dimensional
phasespace for the profiles plotted in figure 8. These pairs are directly
calculated from the evolved profiles at $a=1$ plotted in figure 8 for all
three values of $Q$. The red (blue) lines correspond to profiles I and II with
$\sigma=0.0025$ and $0.1$ respectively. Note that, when $Q>>1$, although cases
I and II have different profiles for $\Delta$ and $\delta_{v}$, the pairs,
when plotted in phasespace, trace out the same curve. Thus, this curve is an
invariant of the dynamics described by eqs. (73) and (74). Similarly, in the
weak field limit, when $Q<<1$, the $\Delta-\delta_{v}$ pairs trace out the
same curve for profiles I and II. In this case, this curve coincides with the
GR limit and is an invariant of the dynamics given by eqs. (69) and (70). In
the intermediate regime, when $Q\sim 1$, the $\Delta-\delta_{v}$ curve is
different for profiles I and II. This is related to the fact that the dynamics
in this regime cannot be described by a closed form dynamical system of the
kind discussed above.

Figure 9: The spherically averaged density $\Delta$ and fractional Hubble
parameter $\delta_{v}$ plotted on the $\Delta-\delta_{v}$ phase space for
profiles I (red, dotted) and II (blue, dot-dashed). In standard GR, the
‘slaving condition’ imposes a specific relation between the non-linear density
and velocity which traces out a special curve in the phasespace. This curve is
an invariant of the dynamical system described by eqs. (69) and (70). In the
$f(R)$ model, in the strong field regime ($Q>>1$), the evolution is scale-
independent and the dynamics can be described by eqs. (73) and (74) akin to
standard GR and the non-linear relation between $\Delta$ and $\delta_{v}$ is
profile independent. In the weak field limit ($Q<<1$), the curves coincide
with the GR ($\Lambda$CDM) case. However, in the intermediate regime ($Q\sim
1$), the Newtonian force is scale dependent and in the $\Delta-\delta_{v}$
plane, the dynamics is profile dependent.

There are various parameters that have been used in the literature to
characterize the departure of modified gravity models from GR (see for example
Silvestri et al. 2013; Lombriser 2016) and it is useful to correlate the
departure of the DVDR relation from GR with such parameters. We define the
parameter (in real space)

$\eta=\frac{\Psi}{\Phi}.$ (75)

This is distinguished from $\eta_{k}=\frac{\Phi_{k}}{\Psi_{k}}$ defined by
Pogosian & Silvestri (2008) which is ratio of the Fourier components of $\Phi$
and $\Psi$. $\eta_{k}$ varies smoothly from $0.5$ to 1 while traversing the
limit from modified gravity to GR. We found that the real space $\eta$ was
constant along the radial direction for the entire range of the perturbation
rising sharply to $\infty$ after the compensated region where
$\Psi\sim\Phi\sim 0$, since it is ill-defined on that domain. In comparing the
deviations, we consider only the domain of the perturbation where $\eta<1$. In
figure 10, we plot the deviation of the curves in phasespace defined as
follows. At a given radial point $x$, the pair $\\{\Delta,\delta_{v,fR}\\}$ is
well-defined from the numerical density and velocity profiles evolved in the
$f(R)$ theory and plotted in figure 8. The pair $\\{\Delta,\delta_{v,GR}\\}$
corresponds to the expected value of $\delta_{v}$ if the evolution assumed
standard GR. The latter was computed by evolving the same initial profile in
standard GR and interpolating the resulting $\Delta-\delta_{v}$ curve; it can
also be computed using the fit given in N13. The deviation between the two
curves is given by

${\rm deviation}(x)=|\delta_{v,GR}(\Delta(x))-\delta_{v,fR}(\Delta(x))|.$ (76)

The red and blue lines correspond to profiles I and II respectively. It is
clear that there is a direct correspondence between $\eta$ and the deviation
in phasespace. Note that $\eta$ and the deviation both remain approximately
constant throughout the radius of the perturbation, although the density,
velocity and potentials are radially varying. As expected, $\eta$ closer to
unity implies as smaller deviation.

Figure 10: Deviation of the DVDR relation in modified gravity from that in GR
at $a=1$ corresponding to the curves in figure 9. The response of the velocity
field to the matter is different in GR and modified gravity. The difference is
directly correlated to the difference between the two potentials $\Psi$ (which
governs dynamics) and $\Phi$ (which governs the curvature). It is interesting
to note that while $\Phi$ and $\Psi$ vary along the radial direction (see
figure 5), their ratio is a constant and directly correlated with deviation of
the DVDR curve from its GR counterpart in phasespace.

###  6.3 Fitting form for the DVDR in the strong field regime.

The non-linear density velocity divergence relation has been extensively
investigated in the past using perturbation theory both in the Eulerian
(Bernardeau, 1992; Chodorowski, 1997; Chodorowski & Lokas, 1997) and
Lagrangian frames (Susperregi & Buchert, 1997) as well as numerical
simulations (Bernardeau & van de Weygaert, 1996; Chodorowski et al., 1998;
Zaroubi et al., 1999; Kudlicki et al., 2000; Ciecielg et al., 2003; Kitaura et
al., 2012). The formula given by Bernardeau (1992) (B92 hereafter) in standard
GR reads

$\Theta=\frac{3\Omega_{m}^{0.6}}{2}\left\\{1-(1+\delta)^{2/3}\right\\}$ (77)

and is valid in the regime $-1\leq\delta\leq 2$. Bilicki & Chodorowski (2008),
BC08 hereafter, gave a fitting form valid for a larger range of $\delta$ but
is based on the exact analytic solution of the spherical collapse and hence
was derived for pure matter cosmologies. N13 gave a formula based on the B92
and BC08, but with the growth index modified to account for a general dark
energy component with a constant equation of state $w$. Mandal & Nadkarni-
Ghosh (2020) checked that this formula also holds for early dark energy
scenarios where $w$ may vary with time.

In this paper we find a fitting form for DVDR in the strong field regime of
$f(R)$ where the dynamics is scale independent. We evaluated the
$\Delta-\delta_{v}$ pairs for final density and velocity profiles in the range
from $0.5\leq a\leq 1$ using a function of the form
$\delta_{v}=A(\Omega_{m})[1-(1+\Delta)^{B(\Omega_{m}})]$ where $A$ and $B$
were coefficients which depend on the instantaneous $\Omega_{m}$. $B$ was
found to be independent of $\Omega_{m}$, whereas $A$ varies as a power law.
The form

$\delta_{v}=0.64\times\Omega_{m}^{0.54}[1-(1+\Delta)^{0.61}]$ (78)

gave a good fit over the range $0\leq\Delta\leq 2$ and $0.5\leq a\leq 1$. In
terms of $\delta$ and $\Theta$ this relation reads

$\Theta=1.92\times\Omega_{m}^{0.54}[1-(1+\delta)^{0.61}].$ (79)

The r.m.s. relative error between the fit and the data was about 17% over the
entire range of $\Delta$ and $a$. In this regime, we do not expect the non-
linear DVDR to depend on the exact choice of $a_{switch}$. We did not attempt
to fit the void region since the top-hats considered here are overdense giving
a positive $\Delta$.

Figure 11 shows the fitting form for the DVDR relation for the $\Lambda$CDM
case (black) and the $f(R)$ model (red) in the strong field regime as a
function of redshift. The red points represent the numerically evolved
$(\Delta,\delta_{v})$ points corresponding to the selected epoch. The blue
arrows show the streamlines of the flow given by eqs. (73) and (74). They show
the instantaneous direction of the vector defined by
$(\Delta^{\prime},\delta_{v}^{\prime})$. The system defined by eqs. (73) and
(74) is a non-autonomous system since the coefficients of $\Delta$ and
$\delta_{v}$ in the r.h.s. are time-dependent. Hence, the streamlines are not
constant at all epochs and neither is the DVDR curve. This can be contrasted
with the case of a pure $\Omega_{m}=1$ flat universe, where the dynamical
system for $\Delta-\delta_{v}$ is an autonomous system and the DVDR relation
is time-independent (see N13 for a detailed discussion).

Figure 11: The DVDR in the strong field limit ($Q>>1$). The red line is the
fitting form given in eq. (78) and the red dots show the $\Delta-\delta_{v}$
pairs computed from the radial profile at the selected epoch. The solid black
line shows the DVDR for the reference $\Lambda$CDM model considered here and
the black dotted line shows the DVDR for the flat EdS model with
$\Omega_{m}=1$. The dynamical system of equations governing $\Delta$ and
$\delta_{v}$ is non-autonomous in the case of $\Lambda$CDM and $f(R)$ and
hence the curves are time-dependent. On the other hand, the EdS curve is the
same for all epochs.

The DVDR relation based on spherical collapse is one-to-one. However, for
realistic density and velocity fields, there is a scatter around a mean
relation. Numerical simulations in the case of $\Lambda$CDM have shown that
the mean relation is close to that given by the spherical collapse model
(Bernardeau et al., 1999; Kudlicki et al., 2000; Bilicki & Chodorowski, 2008).
This scatter is usually attributed to the non-local nature of gravity. The
scatter can be related to the shear component of the velocity field as was
shown by Chodorowski (1997) using perturbation theory and using triaxial
collapse by Nadkarni-Ghosh & Singhal (2016). It remains to be checked whether
the mean relation from the joint density-velocity PDF obtained from
simulations agrees with the estimate based on spherical collapse given in this
paper.

###  6.4 Density enhancement in non-linear evolution

In this paper, we do not undertake a systematic study of the highly non-linear
regime, in particular the process of collapse. However, we illustrate,
qualitatively, the difference between highly non-linear evolution in the
strong and weak limits. We consider two profiles: both with $\sigma=0.0025$,
but with different $x_{top}$. One case corresponds to $x_{top}=0.002h^{-1}$Mpc
and the other with $x_{top}=200h^{-1}$Mpc. $Q>>1$ for the former and $Q<<1$
for the latter over the entire range of evolution from $a=0.1$ to $a=1$. In
each case, the initial amplitude was chosen so that the final evolved density
in the $f(R)$ model has $\delta\sim 100$. This implies that in the strong
field limit, the initial amplitude of the perturbation is less than that for
the weak field limit. The details are in appendix A. The resulting non-linear
$\Delta$ is shown in the top-panel of figure 12 and the difference between the
$f(R)$ and GR profiles is shown in the bottom panel. For $Q>>1$, the
amplification is large since the effective force is $4/3$ times stronger. For
$Q<<1$, the evolved profiles for $f(R)$ and GR apparently coincide. However,
the lower panel shows that the $f(R)$ model has a higher density at the edge
of the top-hat. This can be explained in terms of the solution for $\chi$ as
was discussed in §4.6. $\chi$ is proportional to $\delta$ in this regime and
the force which is proportional to the gradient is non-zero only at the edge
of the top-hat.

The density enhancement at the edge has been reported earlier (Borisov et al.,
2012; Kopp et al., 2013; Lombriser, 2016) but in the presence of the chameleon
screening. The appearance of a ‘spike’ at the edge of the top-hat is akin to
the presence of a ‘thin-shell’ in the chameleon mechanism, although in this
paper we have not explicitly modelled the chameleon screening. Hence, it is
worth understanding the relation between these two features. Following the
seminal paper on ‘chameleon cosmology’ (Khoury & Weltman, 2004), we consider a
scalar field $\phi$, evolving according to

$\nabla_{r}^{2}\phi=V_{eff}(\phi)_{,\phi},$ (80)

where $V_{eff}=V(\phi)+\rho e^{\sqrt{8\pi G}\beta\phi}$ and $\beta$ is a
coupling constant of order unity. This equation is solved for an overdense
sphere of density $\rho$ embedded in a background of density
$\rho_{\infty}<\rho$, with boundary conditions
$\phi(r\rightarrow\infty)=\phi_{\infty}$, where $\phi_{\infty}$ is the value
of the field that minimises the effective potential at $\infty$ and
$d\phi/dr(r\rightarrow 0)=0$. The ‘chameleon solution’ refers to a special
solution of this equation where the effective potential is at its minimum
almost everywhere inside the sphere, except for a thin shell at the surface of
the sphere. This solution corresponds to $\nabla_{r}^{2}\phi\approx
V_{eff,\phi}\approx 0$ and the ‘extra fifth force’, which is given by
$\nabla\phi$ is zero inside the sphere except near the surface. When $\rho$ is
large, as is the case inside the sphere, the second derivative of the
effective potential at the minimum is larger i.e., the potential is steeper.
The mass of the scalar field for this solution is large and correspondingly,
the Compton wavelength is small inside the sphere. Outside the sphere, the
density is low, the effective potential is shallower, the Compton wavelength
is large and the fifth force is not screened.

In the $f(R)$ model, the scalar field corresponds to $\chi=\delta f_{R}$. The
$Q<<1$ solution corresponds to $\nabla_{x}^{2}\chi\approx 0$, which gives
$\chi={\bar{x}}_{c}^{2}H^{2}\Omega_{m}a^{2}\delta$. Since the ‘fifth force’
depends on $\nabla_{x}\chi$, it is zero everywhere except at the edge of the
top-hat, where the derivative of the density is non-zero. However, we do not
refer to this as ‘chameleon screening’ since the Compton wavelength is
constant both inside and outside the top-hat and is independent of $\delta$.
This can also be seen by defining an ‘effective potential’
$V_{eff}(\chi)=\frac{\chi^{2}}{2{\bar{x}}_{C}^{2}}+H^{2}\Omega_{m}a^{2}\delta\chi$,
recasting the equation for $\chi$ as
$\nabla_{x}^{2}\chi=V_{eff}(\chi)_{,\chi}$ and noting that the second
derivative, $V_{eff,\chi\chi}$, which is inversely related to the Compton
wavelength, is independent of $\delta$. In this paper, we have used the first
order of the Taylor expansion $\delta f_{R}\sim f_{RR}\delta R+\mathcal{O}(R)$
in obtaining a linear equation for $\chi$. Ideally, as has been argued by Hu-
Sawicki, this linearization is not valid when $\delta R$ is large (high
curvature limit). Alternately, it is clear that when the Compton wavelength is
small i.e., $f_{RR}$ is small, higher order corrections to this expansion may
become important making the equation for $\chi$ non-linear. Higher order terms
can be modelled by allowing the Compton wavelength to depend on $\chi$. Thus,
the chameleon solution is essentially a solution of the non-linear equation
for $\chi$.

Physically, the density spike at the edge is related to the enhanced force on
the shells near the edge. For a constant density top-hat, all shells up to the
edge experience the same acceleration and the top-hat maintains its shape
until the innermost shell collapses first. On the other hand, in $f(R)$
gravity, even though at the initial time, all shells start with the same
acceleration, as evolution proceeds, the shells near the edge accelerate
faster than they would in GR. The inner shells accelerate at the same rate as
GR since the fifth force is screened inside. Mass starts to accumulate at the
edge giving rise to an enhanced density, until, eventually, the faster moving
shells cross the inner slower shells giving rise to a caustic at the edge.
This was observed in Borisov et al. (2012). It is important to note that this
density enhancement is pronounced for a top-hat due to the presence of a steep
gradient. The ratio of the additional force to the GR force is

$\frac{\nabla_{x}\chi}{\nabla_{x}\Phi_{GR}}\sim{\bar{x}}_{C}\frac{d\ln\delta}{dx}.$
(81)

Enhancement will be prominent only when the density changes over a scale much
smaller than the Compton wavelength of the scalar field. For a smooth radial
density, enhancement will be minimal in the weak field regime. Thus, the edge
effect is observed when the scale of the density perturbation is larger than
the Compton wavelength, but the scale of the density gradients are larger than
the Compton wavelength. In the strong field regime, when $Q>>1$, the the
Compton wavelength is large compared to the scale of the perturbation. There
are no edge effects since the ‘fifth force’ and the associated density
enhancement is uniform; the final density is higher as compared to GR. When
$Q\sim 1$, there is a slight edge effect (see middle panel of $\nabla\chi$ in
figure 5), but is suppressed compared to the $Q<<1$ case because the Compton
wavelength is larger. Understanding this behaviour for the full non-linear
equation is left for future work.

Figure 12: Demonstration of non-linear growth in strong vs. weak regimes. In
the upper panel, the red dashed denotes $f(R)$ and black solid is the GR value
for $\Delta$ at $a=1$. In the strong regime $Q>>1$, the profile is amplified
as compared to GR, whereas for the weak field regime, where $Q<<1$, the two
profiles coincide by eye. The lower panel profile develops a cusp at the edge
of the top-hat. This is related to the solution for $\chi$ which is
proportional to the density in the weak field limit. Since the force is
proportional to the $\nabla\chi$, this manifests itself as a cusp at the edge
of the top-hat.

##  7 Summary and discussion

In this work, we have investigated the non-linear evolution of cosmological
perturbations in $f(R)$ models using the compensated spherical top-hat as a
proxy for the non-linear regime. Our algorithm is an iterative hybrid
Lagrangian-Eulerian scheme; the continuity equation and Euler equation are
solved in Lagrangian coordinates and the equations for the metric potentials
are purely spatial and are solved analytically in Eulerian space.
Specifically, we chose the Hu & Sawicki (2007) model for simplicity, with
majority of the results obtained by setting $f_{R0}=-10^{-6}$ and $n=1$. The
evolution equations depend upon the ratio of $Q={\bar{x}_{C}}/{x_{top}}$,
where ${\bar{x}_{C}}$ is the reduced comoving Compton wavelength of the model
given by eq. (23). The Compton wavelength is assumed to depend only on the
background cosmology; thus we do not model the chameleon screening mechanism
(Khoury & Weltman 2004). However, we demonstrate the role played by the
Compton wavelength by changing the scale of the top-hat ($x_{top}$). We
consider three values of $x_{top}=0.002,2$ and 200 $h^{-1}{\rm Mpc}$,
corresponding to $Q=1217,1.21$ and 0.012 at $a=1$. There are several new
features of this work.

  * •

Analysis of oscillations in the background evolution
It is well known that the $f(R)$ model exhibits high frequency oscillations at
high redshifts, both in the background and perturbation evolution making the
system numerically stiff. To overcome this, it is favourable to assume GR at
early epochs and ‘switch’ to $f(R)$ evolution at a late epoch denoted by
$a_{switch}$. By performing an eigenvalue analysis of the linearized
background equations, we were able to assess the nature of oscillations and
make an informed choice of $a_{switch}$. Such an eigenanalysis has been
performed for the Boltzmann hierarchy in order to gain insights for optimizing
solvers (Nadkarni-Ghosh & Refregier 2017) and a similar analysis can also be
performed for the linear perturbation system of $f(R)$ or other models of
modified gravity. In particular, it could help in putting the Quasi-Static-
Approximation on firmer grounds (Hojjati et al. 2012; Sawicki & Bellini 2015;
Pace et al. 2021).

  * •

Analytic solutions for the metric potentials and density-enhancement
Growth of perturbations in standard GR is governed by the joint continuity,
Euler and Poisson equations. In $f(R)$ models, there is an extra equation
corresponding to the ‘extra’ degree of freedom, usually encapsulated by the
variable $f_{R}$. In the limit $|f_{R}|<<1$ and $|f/R|<<1$, the perturbations
in $f_{R}$ can be expressed as the difference of the two potentials
$\Phi-\Psi$, where $\Phi$ and $\Psi$ are the metric perturbations
corresponding to the ‘curvature’ potential and ‘Newtonian’ potential
respectively. We derive an analytic solution for $\chi=\Phi-\Psi$ in real
space in spherical symmetry valid for all values of $Q$. In the strong field
limit, when the scale of the perturbation is much less than its Compton
wavelength ($Q>>1$), $\chi$ obeys the Poisson equation and $\Psi$ is 4/3 times
its value in standard GR. In the weak field limit, when the scale of the
perturbation is much greater than its Compton wavelength ($Q<<1$),
$\chi\propto\delta$. This means that the force, which is proportional to the
gradient of $\chi$ is non-zero only at the edge of the top-hat. This results
in a density enhancement at the top-hat edge, a phenomenon that has been
previously reported in the literature (Borisov et al., 2012; Kopp et al.,
2013) in the context of ‘chameleon screening’. Here we demonstrate the same
effect in the absence of a chameleon mechanism. If the Compton wavelength is
allowed to change within the perturbation, then the resulting equation for
$\chi$ is non-linear, possibly enhancing this effect. A more quantitative
exploration is left for future investigations. Nevertheless, the analytic
solutions presented here can also serve as a ‘test case’ to check other
numerical schemes, similar in spirit to the use of spherical collapse as a
‘control case’ for numerical simulations (see Joyce & Labini 2012).

  * •

Phasespace analysis of evolution of perturbations and the density-velocity
divergence DVDR relation in the strong field limit
We evolve the compensated top-hats and examine the evolution of the
perturbations in the 2-D density-velocity divergence phase space. For the
radially varying profiles, it is easier to compute the $\Delta-\delta_{v}$
paris instead of the $\delta-\Theta$ pairs, where $\Delta$ is the spherically
averaged density and $\delta_{v}$ is the fractional Hubble velocity of the
shell. In the strong field limit, where the evolution of perturbations is
scale-invariant, it is possible to define a unique DVDR relation and we give a
fitting form for the same. In the intermediate field limit, $Q\sim 1$, the
evolution of perturbations is profile dependent and an invariant relation
between $\delta$ and $\Theta$ does not exist. In the weak field limit, the
dynamics is approximately that of GR.

There are several caveats in the present work. The spherical collapse is a
simplistic model: the system does not take into account any effect of the
environment or rotations and ignores mode-mode coupling. The first step
forward is to consider triaxial collapse. Triaxial collapse, which is an
important ingredient in mock catalog generators like PINOCCHIO (Monaco, 2016;
Rizzo et al., 2017; Moretti et al., 2020; Song et al., 2021) has not been
studied in great detail in the context of modified gravity, albeit, with the
exception of some recent efforts (Burrage et al. 2015, 2018; Ruan et al.
2020). To capture the effects of mode coupling in the analytic framework, it
is necessary to consider perturbation theory, either in the Eulerian or
Lagrangian frame and related schemes (see e.g., Aviles & Cervantes-Cota 2017;
Aviles et al. 2019; Valogiannis et al. 2020; Aviles et al. 2021; Baojiu 2018).
The hybrid Eulerian-Lagrangian scheme outlined here does not rely on spherical
geometry. Equations (31)-(34) are general and the hybrid scheme can also be
incorporated with a multi-step Lagrangian perturbation theory (Nadkarni-Ghosh
& Chernoff, 2011, 2013), which guarantees convergence in $\Lambda$CDM models
up until shell crossing. However, convergence and shell-crossing related
issues in modified gravity scenarios remain to be investigated (see for
example Rampf & Frisch 2017; Rampf & Hahn 2020 for a discussion of these
matters in standard GR). A code, called SELCIE, based on Finite Element
Methods to solve the chameleon equations of motion in arbitrary mass
distributions has recently been developed by Briddon et al. (2021). However,
it is currently not set up to perform temporal evolution of the matter
distribution in the presence of chameleon screening. It may be possible to
couple iterative methods used in this paper with such tools. In addition, for
the $f(R)$ model considered here, the relevant scales for which the system is
in the strong field limit, are not cosmological scales, but correspond to
smaller bodies governed by additional astrophysical processes which have also
been ignored in this treatment.

The $f(R)$ model, in the metric formalism can be considered as a special case
of Brans-Dickie theories, which in turn belongs to the more general class of
Horndeski theories (De Felice & Tsujikawa, 2010; Kobayashi, 2019). More
recently, the Effective Field Theory of dark energy has become a popular
formalism which is aimed at capturing a wide class of modified theory models
(see review by Frusciante & Perenon 2020). Other kinds of $f(R)$ modifications
have also been introduced in the context of understanding galaxy rotation
curves (for e.g., Dey, Bhattacharya, & Sarkar 2013, 2015) or in the context of
unified models of dark energy and inflation (for e.g., Nojiri & Odintsov 2007;
Cognola et al. 2008; Nojiri & Odintsov 2011). Spherically symmetric solutions
play an important role in understanding the collapse processes in such systems
(Astashenok et al., 2019; Chowdhury et al., 2020; Jaryal & Chatterjee, 2021).
Furthermore, many future tests of modified gravity will involve measurements
on a very wide range of astrophysical and cosmological scales depending upon
the magnitude of modification to the ‘curvature’ and ‘potential’ fields (Baker
et al., 2015). Many of these systems are often modelled by spherical symmetry.
We hope that some of the techniques presented in this work may be extended to
gain insights into these generalized models of modified gravity in such
systems.

##  8 Acknowledgements

S.N. would like to acknowledge the Department of Science and Technology, Govt.
of India, for the grant no. WOS-A/PM-21/2018. S.C. would like to acknowledge
the Council of Scientific & Industrial Research (CSIR) for the grant no.
09/092(0930)/2015-EMR-I. Both authors would like to thank Tapobrata Sarkar for
useful discussions through the course of this project and for detailed
comments on the manuscript. S. N. would like to thank Alessandra Silvestri for
pointing out some useful references.

##  9 Data Availability

Data generated by the numerical runs available on request from the
corresponding author.

## References

  * Alam et al. (2017) Alam S. et al., 2017, Monthly Notices of the Royal Astronomical Society, 470, 2617
  * Alam et al. (2020) Alam S. et al., 2020, arXiv e-prints, 2011, arXiv:2011.05771
  * Amendola et al. (2018) Amendola L. et al., 2018, Living Reviews in Relativity, 21, 2
  * Appleby & Battye (2007) Appleby S., Battye R., 2007, Physics Letters B, 654, 7
  * Astashenok et al. (2019) Astashenok A. V., Mosani K., Odintsov S. D., Samanta G. C., 2019, International Journal of Geometric Methods in Modern Physics, 16, 1950035
  * Aviles & Cervantes-Cota (2017) Aviles A., Cervantes-Cota J. L., 2017, Physical Review D, 96
  * Aviles et al. (2019) Aviles A., Cervantes-Cota J. L., Mota D. F., 2019, Astronomy & Astrophysics, 622, A62
  * Aviles et al. (2021) Aviles A., Valogiannis G., Rodriguez-Meza M. A., Cervantes-Cota J. L., Li B., Bean R., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 039
  * Baker et al. (2021) Baker T. et al., 2021, Reviews of Modern Physics, 93, 015003
  * Baker et al. (2015) Baker T., Psaltis D., Skordis C., 2015, The Astrophysical Journal, 802, 63
  * Baojiu (2018) Baojiu L., 2018, International Journal of Modern Physics, 27
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, The Astrophysical Journal, 304, 15
  * Barreira et al. (2013) Barreira A., Li B., Baugh C. M., Pascoli S., 2013, Journal of Cosmology and Astro-Particle Physics, 11, 056
  * Bean et al. (2007) Bean R., Bernat D., Pogosian L., Silvestri A., Trodden M., 2007, Physical Review D, 75, 064020
  * Bernardeau (1992) Bernardeau F., 1992, The Astrophysical Journal, 390, L61
  * Bernardeau et al. (1999) Bernardeau F., Chodorowski M. J., Lokas E. L., Stompor R., Kudlicki A., 1999, Monthly Notices of the Royal Astronomical Society, 309, 543
  * Bernardeau & van de Weygaert (1996) Bernardeau F., van de Weygaert R., 1996, Monthly Notices of the Royal Astronomical Society, 279, 693
  * Bertschinger & Dekel (1989) Bertschinger E., Dekel A., 1989, The Astrophysical Journal, 336, L5
  * Bertschinger & Zukin (2008) Bertschinger E., Zukin P., 2008, Physical Review D, 78, 024015
  * Bilicki & Chodorowski (2008) Bilicki M., Chodorowski M. J., 2008, Monthly Notices of the Royal Astronomical Society, 391, 1796
  * Borisov et al. (2012) Borisov A., Jain B., Zhang P., 2012, Physical Review D, 85, 63518
  * Brenier et al. (2003) Brenier Y., Frisch U., Hénon M., Loeper G., Matarrese S., Mohayaee R., Sobolevskĩi A., 2003, Monthly Notices of the Royal Astronomical Society, 346, 501
  * Briddon et al. (2021) Briddon C., Burrage C., Moss A., Tamosiunas A., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 043
  * Buchert (1992) Buchert T., 1992, Monthly Notices of the Royal Astronomical Society, 254, 729
  * Burrage et al. (2018) Burrage C., Copeland E. J., Moss A., Stevenson J. A., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 056?056
  * Burrage et al. (2015) Burrage C., Copeland E. J., Stevenson J. A., 2015, Physical Review D, 91
  * Cataneo et al. (2021) Cataneo M., Uhlemann C., Arnold C., Gough A., Li B., Heymans C., 2021, The matter density pdf for modified gravity and dark energy with large deviations theory
  * Ceron-Hurtado et al. (2016) Ceron-Hurtado J. J., He J.-h., Li B., 2016, Physical Review D, 94, 064052
  * Chakrabarti & Banerjee (2016) Chakrabarti S., Banerjee N., 2016, General Relativity and Gravitation, 48, 57
  * Chakraborty et al. (2021) Chakraborty S., MacDevette K., Dunsby P., 2021, arXiv e-prints, 2103, arXiv:2103.02274
  * Chiba (2003) Chiba T., 2003, Physics Letters B, 575, 1
  * Chodorowski (1997) Chodorowski M. J., 1997, Monthly Notices of the Royal Astronomical Society, 292, 695
  * Chodorowski & Lokas (1997) Chodorowski M. J., Lokas E. L., 1997, Monthly Notices of the Royal Astronomical Society, 287, 591
  * Chodorowski et al. (1998) Chodorowski M. J., Lokas E. L., Pollo A., Nusser A., 1998, Monthly Notices of the Royal Astronomical Society, 300, 1027
  * Chowdhury et al. (2020) Chowdhury S., Pal K., Pal K., Sarkar T., 2020, European Physical Journal C, 80, 902
  * Ciecielg et al. (2003) Ciecielg P., Chodorowski M. J., Kiraga M., Strauss M. A., Kudlicki A., Bouchet F. R., 2003, Monthly Notices of the Royal Astronomical Society, 339, 641
  * Clifton et al. (2012) Clifton T., Ferreira P. G., Padilla A., Skordis C., 2012, Physics Reports, 513, 1
  * Cognola et al. (2008) Cognola G., Elizalde E., Nojiri S., Odintsov S. D., Sebastiani L., Zerbini S., 2008, Physical Review D, 77
  * Colombi et al. (2007) Colombi S., Chodorowski M. J., Teyssier R., 2007, Monthly Notices of the Royal Astronomical Society, 375, 348
  * Dai et al. (2008) Dai D.-C., Maor I., Starkman G., 2008, Physical Review D, 77, 64016
  * De Felice & Tsujikawa (2010) De Felice A., Tsujikawa S., 2010, Living Reviews in Relativity, 13, 3
  * Dey et al. (2013) Dey D., Bhattacharya K., Sarkar T., 2013, Physical Review D, 87
  * Dey et al. (2015) Dey D., Bhattacharya K., Sarkar T., 2015, General Relativity and Gravitation, 47
  * Di Valentino et al. (2021) Di Valentino E. et al., 2021, Classical and Quantum Gravity, 38, 153001
  * Dodelson (2003) Dodelson S., 2003, Modern Cosmology. Academic Press
  * Dvali et al. (2000) Dvali G., Gabadadze G., Porrati M., 2000, Physics Letters B, 485, 208
  * Elizalde et al. (2012) Elizalde E., Odintsov S. D., Sebastiani L., Zerbini S., 2012, European Physical Journal C, 72, 1843
  * Frusciante & Perenon (2020) Frusciante N., Perenon L., 2020, Physics Reports, 857, 1?63
  * Gramann (1993a) Gramann M., 1993a, The Astrophysical Journal, 405, 449
  * Gramann (1993b) Gramann M., 1993b, The Astrophysical Journal Letters, 405, L47
  * Hahn et al. (2015) Hahn O., Angulo R. E., Abel T., 2015, Monthly Notices of the Royal Astronomical Society, 454, 3920
  * Herrera et al. (2017) Herrera D., Waga I., Jorás S., 2017, Physical Review D, 95
  * Hojjati et al. (2012) Hojjati A., Pogosian L., Silvestri A., Talbot S., 2012, Physical Review D, 86
  * Hu & Sawicki (2007) Hu W., Sawicki I., 2007, Physical Review D, 76, 104043
  * Jaryal & Chatterjee (2021) Jaryal S. C., Chatterjee A., 2021, The European Physical Journal C, 81
  * Joyce et al. (2015) Joyce A., Jain B., Khoury J., Trodden M., 2015, Physics Reports, 568, 1
  * Joyce et al. (2016) Joyce A., Lombriser L., Schmidt F., 2016, Annual Review of Nuclear and Particle Science, 66, 95
  * Joyce & Labini (2012) Joyce M., Labini F. S., 2012, Monthly Notices of the Royal Astronomical Society, 429, 1088
  * Khoury & Weltman (2004) Khoury J., Weltman A., 2004, Physical Review D, 69, 044026
  * Kitaura et al. (2012) Kitaura F.-S., Angulo R. E., Hoffman Y., Gottlöber S., 2012, Monthly Notices of the Royal Astronomical Society, 425, 2422
  * Kobayashi (2019) Kobayashi T., 2019, Reports on Progress in Physics, 82, 086901
  * Kopp et al. (2013) Kopp M., Appleby S. A., Achitouv I., Weller J., 2013, Physical Review D, 88, 084015
  * Kudlicki et al. (2000) Kudlicki A., Chodorowski M., Plewa T., Różyczka M., 2000, Monthly Notices of the Royal Astronomical Society, 316, 464
  * Li & Efstathiou (2012) Li B., Efstathiou G., 2012, Monthly Notices of the Royal Astronomical Society, 421, 1431
  * Li & Lam (2012) Li B., Lam T. Y., 2012, Monthly Notices of the Royal Astronomical Society, 425, 730
  * Lilje & Lahav (1991) Lilje P. B., Lahav O., 1991, The Astrophysical Journal, 374, 29
  * Lilow & Nusser (2021) Lilow R., Nusser A., 2021, Monthly Notices of the Royal Astronomical Society
  * Lima & Liddle (2013) Lima N. A., Liddle A. R., 2013, Physical Review D, 88, 043521
  * Linder (2005) Linder E. V., 2005, Physical Review D, 72, 43529
  * Linder & Cahn (2007) Linder E. V., Cahn R. N., 2007, Astroparticle Physics, 28, 481
  * Lombriser (2016) Lombriser L., 2016, Journal of Cosmology and Astro-Particle Physics, 11, 039
  * Ma & Bertschinger (1995) Ma C.-P., Bertschinger E., 1995, The Astrophysical Journal, 455, 7
  * Macaulay et al. (2013) Macaulay E., Wehus I. K., Eriksen H. K., 2013, Physical Review Letters, 111
  * Magnano & Sokołowski (1994) Magnano G., Sokołowski L. M., 1994, Physical Review D, 50, 5039
  * Mancinelli & Yahil (1995) Mancinelli P. J., Yahil A., 1995, The Astrophysical Journal, 452, 75
  * Mancinelli et al. (1993) Mancinelli P. J., Yahil A., Canon G., Dekel A., 1993, in On nonlinear approximations to cosmic problems with mixed boundary conditions, p. 215
  * Mandal & Nadkarni-Ghosh (2020) Mandal A., Nadkarni-Ghosh S., 2020, Monthly Notices of the Royal Astronomical Society, 498, 355
  * Martino et al. (2009) Martino M. C., Stabenau H. F., Sheth R. K., 2009, Phys. Rev. D, 79, 084013
  * Monaco (2016) Monaco P., 2016, Galaxies, 4, 53
  * Moretti et al. (2020) Moretti C., Mozzon S., Monaco P., Munari E., Baldi M., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1153
  * Nadkarni-Ghosh (2013) Nadkarni-Ghosh S., 2013, Monthly Notices of the Royal Astronomical Society, 428, 1166
  * Nadkarni-Ghosh & Chernoff (2011) Nadkarni-Ghosh S., Chernoff D. F., 2011, Monthly Notices of the Royal Astronomical Society, 410, 1454
  * Nadkarni-Ghosh & Chernoff (2013) Nadkarni-Ghosh S., Chernoff D. F., 2013, Monthly Notices of the Royal Astronomical Society, 431, 799
  * Nadkarni-Ghosh & Refregier (2017) Nadkarni-Ghosh S., Refregier A., 2017, Monthly Notices of the Royal Astronomical Society, 471, 2391
  * Nadkarni-Ghosh & Singhal (2016) Nadkarni-Ghosh S., Singhal A., 2016, Monthly Notices of the Royal Astronomical Society, 457, 2773
  * Nojiri et al. (2017) Nojiri S., Odintsov S., Oikonomou V., 2017, Physics Reports, 692, 1–104
  * Nojiri & Odintsov (2007) Nojiri S., Odintsov S. D., 2007, Physics Letters B, 657, 238–245
  * Nojiri & Odintsov (2011) Nojiri S., Odintsov S. D., 2011, Physics Reports, 505, 59–144
  * Nusser (2017) Nusser A., 2017, Monthly Notices of the Royal Astronomical Society, 470, 445
  * Nusser et al. (2012) Nusser A., Branchini E., Davis M., 2012, The Astrophysical Journal, 744, 193
  * Nusser et al. (2014) Nusser A., Davis M., Branchini E., 2014, The Astrophysical Journal, 788, 157
  * Nusser et al. (1991) Nusser A., Dekel A., Bertschinger E., Blumenthal G. R., 1991, The Astrophysical Journal, 379, 6
  * Oyaizu (2008) Oyaizu H., 2008, Physical Review D, 78, 123523
  * Pace et al. (2021) Pace F., Battye R. A., Bellini E., Lombriser L., Vernizzi F., Bolliet B., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 017
  * Paranjape et al. (2013) Paranjape A., Sheth R. K., Desjacques V., 2013, Monthly Notices of the Royal Astronomical Society, 431, 1503–1512
  * Peebles (1980) Peebles P., 1980, The Large-Scale Structure of the Universe. Princeton University Press
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, Astronomy and Astrophysics, 641, A6
  * Pogosian & Silvestri (2008) Pogosian L., Silvestri A., 2008, Physical Review D, 77, 023503
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, Astrophysical Journal, 187, 425
  * Rampf & Frisch (2017) Rampf C., Frisch U., 2017, Monthly Notices of the Royal Astronomical Society, 471, 671
  * Rampf & Hahn (2020) Rampf C., Hahn O., 2020, Monthly Notices of the Royal Astronomical Society: Letters, 501, L71
  * Rizzo et al. (2017) Rizzo L. A., Villaescusa-Navarro F., Monaco P., Munari E., Borgani S., Castorina E., Sefusatti E., 2017, Journal of Cosmology and Astroparticle Physics, 2017, 008
  * Ruan et al. (2020) Ruan C.-Z., Zhang T.-J., Hu B., 2020, Monthly Notices of the Royal Astronomical Society, 492, 4235
  * Sawicki & Bellini (2015) Sawicki I., Bellini E., 2015, Physical Review D, 92, 084061
  * Schäfer & Koyama (2008) Schäfer B. M., Koyama K., 2008, Monthly Notices of the Royal Astronomical Society, 385, 411
  * Schmidt et al. (2010) Schmidt F., Hu W., Lima M., 2010, Physical Review D, 81
  * Schmidt et al. (2009) Schmidt F., Lima M., Oyaizu H., Hu W., 2009, Physical Review D, 79, 83518
  * Scoccimarro (2004) Scoccimarro R., 2004, Physical Review D, 70, 83007
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, Monthly Notices of the Royal Astronomical Society, 323, 1–12
  * Silvestri et al. (2013) Silvestri A., Pogosian L., Buniy R. V., 2013, Physical Review D, 87, 104015
  * Slotine & Weiping (1991) Slotine J. E., Weiping L., 1991, Applied non-linear control. Prentice-Hall Inc., Englewood Cliffs, New Jersey 07632 USA
  * Song et al. (2021) Song Y., Moretti C., Monaco P., Hu B., 2021
  * Song & Doré (2009) Song Y.-S., Doré O., 2009, Journal of Cosmology and Astroparticle Physics, 2009, 025
  * Song et al. (2007) Song Y.-S., Hu W., Sawicki I., 2007, Physical Review D, 75, 044004
  * Sotiriou & Faraoni (2010) Sotiriou T. P., Faraoni V., 2010, Reviews of Modern Physics, 82, 451
  * Starobinsky (1980) Starobinsky A., 1980, Physics Letters B, 91, 99
  * Starobinsky (2007) Starobinsky A. A., 2007, Soviet Journal of Experimental and Theoretical Physics Letters, 86, 157
  * Strogatz (1994) Strogatz S. H., 1994, Non-linear dynamics and chaos. Perseus Publishing, Cambridge, Massachusets
  * Susperregi & Buchert (1997) Susperregi M., Buchert T., 1997, Astronomy and Astrophysics, 323, 295
  * Uzan (2007) Uzan J.-P., 2007, General Relativity and Gravitation, 39, 307, arXiv: astro-ph/0605313
  * Valogiannis et al. (2020) Valogiannis G., Bean R., Aviles A., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 055?055
  * Wolfram Research (2018) Wolfram Research I., 2018, Mathematica, version 11.3.0 edn. Wolfram Research, Inc., Champaign, Illinois
  * Zaroubi et al. (1999) Zaroubi S., Hoffman Y., Dekel A., 1999, The Astrophysical Journal, 520, 413
  * Zel’Dovich (1970) Zel’Dovich Y. B., 1970, Astronomy and Astrophysics, 5, 84

##  Appendix A Initial Profiles

The compensated top-hat is a one-dimensional function represented as

$\displaystyle\delta_{top}(x)$ $\displaystyle=$ $\displaystyle
A\;\;\;\;\;\;\;\;\;\;0<x\leq x_{top}$ (82) $\displaystyle=$
$\displaystyle-1\;\;\;x_{top}<x\leq x_{u}$ (83) $\displaystyle=$
$\displaystyle 0\;\;\;{\rm otherwise,}$ (84)

where $A$ is the amplitude of the top-hat. $x_{top}$ and $x_{u}$ are the
boundaries of the overdense region and the underdense compensating region and
are related to through conservation of mass as

$(1+A)x_{top}^{3}=x_{u}^{3}.$ (85)

A smooth top-hat is obtained by the following function

$\delta_{smooth}(x)=\frac{1}{\sqrt{2\pi\sigma^{2}x^{2}}}\left\\{\int_{0}^{x_{top}}A\left(e^{-\frac{(x-y)}{2\sigma^{2}}}-e^{-\frac{(x+y)}{2\sigma^{2}}}\right)ydy+\int_{x_{top}}^{x_{u}}(-1)\left(e^{-\frac{(x-y)}{2\sigma^{2}}}-e^{-\frac{(x+y)}{2\sigma^{2}}}\right)ydy\right\\},$
(86)

where $\sigma$ is the smoothing parameter. In this paper we have used two
different profiles with $\sigma=0.1$ and $\sigma=0.0025$. In the numerical
implementation, $L_{min}$ and $L_{max}$ denote the minimum and maximum values
of the grid points along the $x$-direction. The value of $A$ and $x_{top}$ are
chosen according to the calculation performed. $N_{t}$ gives the number of
time steps used to divide the time interval from the initial epoch
$a_{switch}$ to the final epoch $a_{final}$. $N_{s}$ denotes the number of
spatial grid points, equi-spaced in $\ln x$ along the radial direction. Tables
1 and 2 give the list of parameters used in different sections of the paper.

Profile Name | $\sigma$ | $x_{top}$ |  $Q$ (for $f_{R0}=-10^{-6}$) |  $L_{min}$ ($h^{-1}$Mpc) |  $L_{max}$ ($h^{-1}$Mpc)
---|---|---|---|---|---
1a | 0.0025 | 0.002 | 1217 | 0.0001 | 0.1
1b | 0.0025 | 2 | 1.21 | 0.1 | 100
1c | 0.0025 | 200 | 0.012 | 10 | 10000
1d | 0.0025 | 0\. 02 | 121 | 0.001 | 1
1e | 0.0025 | 0.1 | 24.34 | 0.005 | 5
1f | 0.0025 | 0.2 | 12.1 | 0.01 | 10
1g | 0.0025 | 1 | 2.43 | 0.05 | 50
1h | 0.0025 | 10 | 0.24 | 0.5 | 500
1i | 0.0025 | 20 | 0.12 | 1 | 1000
1j | 0.0025 | 100 | 0.024 | 5 | 5000
2a | 0.1 | 0.002 | 1217 | 0.0001 | 0.1
2b | 0.1 | 2 | 1.21 | 0.1 | 100
2c | 0.1 | 200 | 0.012 | 10 | 10000
2d | 0.1 | 0.02 | 121 | 0.001 | 1
Table 1: Table denoting the parameters of the two profiles used. The profiles
1e-1j are used only in figure 6 .

Figure | Profile used | $A$ | $N_{s}$ | $N_{t}$ | $a_{switch}$ | $a_{final}$
---|---|---|---|---|---|---
figure 5 | 1a,b,c | 7 | 3000 | - | - | -
figure 6 | 1b-1j | 7 | 3000 | - | - | -
figure 7 | 1a and 1b | 0.0007 | 500 | 40 | 0.1 | 0.2, 0.4
figures 8, 9, 10 | 1a,b,c and 2a,b,c | 0.0008 | 1600 | 320 | 0.1 | 1
figure 11 | 1a | 0.0008 | 1600 | 320 | 0.1 | 1
figure 12 | 1a and 1c | 0.0015 (for 1c) 0.002 (for 1d) | 1600 | 320 | 0.1 | 1
figure 13 | 1d | 0.0008 | $5\times N_{t}$ | 20 to 640 | 0.1 | 1

Table 2: Table characterizing the initial data used in various figures in the
text. Figures 5 and 6 correspond to static solutions and there are no
parameters for the temporal evolution

.

We found that it was numerically more stable to use compensated top-hats
because the potentials $\Phi$ and $\Psi$ analytically vanish after a finite
extent. For a pure top-hat or a more realistic profile such as the one based
on peaks theory (Bardeen et al., 1986; Lilje & Lahav, 1991) and used in
Nadkarni-Ghosh (2013) or Kopp et al. (2013) is not as stable because the
fields are non-zero at any finite radius (which is inevitable in the
discretization).

##  Appendix B Convergence tests

Figure 13: Convergence test of the algorithm. Comparison of successive
approximations for evolution from $a=0.1$ to $a=1$ for profile 1a for the
variable $\Delta$. The algorithm was applied to a $\Lambda$CDM model with
$\Phi=\Psi$ and an $f(R)$ model with $f_{R0}=-10^{-6}$.

The algorithm outlined in §4.5 was used to evolve an initial profile from
$a=0.001$ to $a=1$ with the transition from $f(R)$ to GR at $a_{switch}=0.1$.
The evolution from $a=0.1$ to $a=1$ was carried out using successive runs with
$N_{t}=20,40,80,160,320$ and 640 steps. For each run, the spatial domain was
divided into $N_{s}$ steps equi-spaced in the $\ln x$ direction. In order to
achieve convergence, the spatial domain needs to be refined along with the
refinement in the temporal domain. We choose $N_{s}=5\times N_{t}$. All
solutions to the second order ordinary differential equations (O.D.E.s) we
computed using the inbuilt function NDSolve in the software package
MATHEMATICA (Wolfram Research, 2018).

In the GR case, the exact answer can be computed in a single step. As a code
test, we computed approximate answers for the GR case using the same algorithm
but assuming $\Phi=\Psi$ and compared to the exact answer as well as compared
successive approximations to give the Cauchy error. For the $f(R)$ evolution
we compared only successive approximations. We define the errors as follows

$\displaystyle{\mathcal{E}}(N)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{N}\sum_{i=1}^{N_{s}}|f_{N}(x_{i})-f_{exact}(x_{i})|^{2}}$
(87) $\displaystyle{\mathcal{E}}_{Cauchy}(N)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{N}\sum_{i=1}^{N_{s}}|f_{N}(x_{i})-f_{2N}(x_{i})|^{2}},$
(88)

where $f$ denotes the function to be compared. The first error cannot be
computed in the $f(R)$ case for lack of knowledge of the exact solution. The
spatial domain has a different grid spacing in each approximation. Thus to
compare, we need to interpolate the approximation which has the finer grid and
compute it at the domain points corresponding to the coarser grid. Figure 13
plots the errors for both GR and $f(R)$ in the evolved average density
$\Delta$ for the profile 1d with the parameters as shown in table 2. We found
the trends to be the same for the variables $\delta$ and $v$ for profile 1e
and the same trends for all three variables for profile 2d.

##  Appendix C Derivation of the analytical solution

The equation for $\Phi_{+}$ is the usual Poisson equation with the solution
given by

$\Phi_{+}(x)=C_{1}+\int_{k_{1}}^{x}\left\\{\frac{C_{2}}{z^{2}}+\frac{1}{z^{2}}\int_{k_{2}}^{z}{\mathcal{P}}_{1}(a)\delta(y)y^{2}dy\right\\}dz,$
(89)

where $y$ and $z$ are dummy variables for the integrals and

${\mathcal{P}}_{1}(a)=\frac{3\Omega_{m}(a)a^{2}}{2}.$ (90)

$k_{1}$ and $k_{2}$ are arbitrary constants which are related to the boundary
points along the $x$ axis. Imposing the boundary condition at $x=\infty$,
gives

$C_{1}=-\int_{k_{1}}^{\infty}\left\\{\frac{C_{2}}{z^{2}}+\frac{1}{z^{2}}\int_{k_{2}}^{z}{\mathcal{P}}_{1}(a)\delta(y)y^{2}dy\right\\}dz,$
(91)

The derivative

$\left.\nabla_{x}{\tilde{\Phi}_{+}}\right|_{x=0}=\frac{C_{2}}{x^{2}}+\frac{1}{x^{2}}\int_{k_{2}}^{y}{\mathcal{P}}_{1}(a)\delta(y)y^{2}dy=0\implies
C_{2}=\int_{0}^{k_{2}}{\mathcal{P}}(a)\delta(y).y^{2}dy$ (92)

Thus, ${\tilde{\Phi}_{+}}$ and its derivative (force) are

$\displaystyle{\tilde{\Phi}_{+}}(x)$ $\displaystyle=$
$\displaystyle{\mathcal{P}}_{1}(a)\int_{\infty}^{x}\frac{1}{z^{2}}\int_{0}^{z}\delta(y)y^{2}dy=\frac{{\mathcal{P}}_{1}(a)}{3}\int_{\infty}^{x}z\Delta(z)dz$
(93) $\displaystyle\nabla_{x}{\tilde{\Phi}_{+}}$ $\displaystyle=$
$\displaystyle{\mathcal{P}}_{1}(a)\frac{1}{x^{2}}\int_{0}^{x}\delta(y)y^{2}dy=\frac{{\mathcal{P}}_{1}(a)}{3}x\Delta(x),$
(94)

where the spherically averaged fractional density contrast $\Delta$ is defined
as:

$\Delta(x)=\frac{3}{x^{3}}\int_{0}^{x}\delta(y)y^{2}dy.$ (95)

The general solution of eq. (47) for ${\tilde{\chi}}$

${\tilde{\chi}}(x)=C_{1}\frac{e^{-\sqrt{{\mathcal{P}}_{2}}x}}{x}+\frac{C_{2}}{2\sqrt{{\mathcal{P}}_{2}}}\frac{e^{\sqrt{{\mathcal{P}}_{2}}x}}{x}+\frac{{{\mathcal{P}}_{3}}}{2\sqrt{{\mathcal{P}}_{2}}}\frac{e^{-\sqrt{{\mathcal{P}}_{2}}x}}{x}\int_{k_{1}}^{x}e^{\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy-\frac{{{\mathcal{P}}_{3}}}{2\sqrt{{\mathcal{P}}_{2}}}\frac{e^{\sqrt{{\mathcal{P}}_{2}}x}}{x}\int_{k_{2}}^{x}e^{-\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy,$
(96)

where

${{\mathcal{P}}_{2}}=\frac{1}{{\bar{x}}_{C}^{2}}\;\;\;{\rm
and}\;\;\;{{\mathcal{P}}_{3}}=a^{2}\Omega_{m}(a).$ (97)

At $x=\infty$, the first and third terms vanish since $\delta(x)$ has a finite
extent. Thus, imposing the boundary condition at $x=\infty$ gives

$C_{2}={{\mathcal{P}}_{3}}\int_{k_{2}}^{\infty}e^{-\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy.$
(98)

This gives

${\tilde{\chi}}(x)=C_{1}\frac{e^{-\sqrt{{\mathcal{P}}_{2}}x}}{x}+\frac{{\mathcal{P}}_{3}}{2\sqrt{{\mathcal{P}}_{2}}x}\left[e^{\sqrt{{\mathcal{P}}_{2}}x}L_{2}(k_{2},\infty)+e^{-\sqrt{{\mathcal{P}}_{2}}x}L_{1}(k_{1},x)-e^{\sqrt{{\mathcal{P}}_{2}}x}L_{2}(k_{2},x)\right]$
(99)

where

$\displaystyle L_{1}(k_{1},x)$ $\displaystyle=$
$\displaystyle\int_{k1}^{x}e^{\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy$ (100)
$\displaystyle L_{2}(k_{2},x)$ $\displaystyle=$
$\displaystyle\int_{k_{2}}^{x}e^{-\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy$
(101)

Taking the derivative

$\displaystyle\frac{d{\tilde{\chi}}}{dx}=$ $\displaystyle-$
$\displaystyle\frac{C_{1}e^{-\sqrt{{\mathcal{P}}_{2}}x}}{x^{2}}(1+\sqrt{{\mathcal{P}}_{2}}x)$
(102) $\displaystyle+$
$\displaystyle\frac{{\mathcal{P}}_{3}}{2\sqrt{{\mathcal{P}}_{2}}x^{2}}\left[L_{2}(k_{2},\infty)(\sqrt{\mathcal{P}}_{2}x-1)e^{\sqrt{\mathcal{P}}_{2}x}-L_{1}(k,x)(1+\sqrt{\mathcal{P}}_{2}x)e^{-\sqrt{\mathcal{P}}_{2}x}-L_{2}(k,x)(\sqrt{\mathcal{P}}_{2}x-1)e^{\sqrt{\mathcal{P}}_{2}x}\right]$
(103)

Setting $\left.\frac{d{\tilde{\chi}}}{dx}\right|_{x\approx 0}=0$ gives

$C_{1}=-\frac{{\mathcal{P}}_{3}}{2\sqrt{{\mathcal{P}}_{2}}}\left[L_{1}(k_{1},0)-L_{2}(k_{2},0)+L_{2}(k_{2},\infty)\right]$
(104)

Substituting for $C_{1}$ and $C_{2}$, writing

$-L_{2}(k_{2},0)+L_{2}(k_{2},\infty)=\int_{0}^{\infty}e^{-\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy=\int_{0}^{x}e^{-\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy+\int_{x}^{\infty}e^{-\sqrt{{\mathcal{P}}_{2}}y}y\delta(y)dy$
(105)

and rearranging the terms gives

$\displaystyle{\tilde{\chi}}(x)$ $\displaystyle=$
$\displaystyle\frac{\mathcal{P}_{3}}{\sqrt{\mathcal{P}_{2}}}\frac{1}{x}\left[\exp(-\sqrt{\mathcal{P}_{2}}x)I_{1}(x)+\sinh(\sqrt{\mathcal{P}_{2}}x)I_{2}(x)\right],$
(106) $\displaystyle\nabla_{x}{\tilde{\chi}}$ $\displaystyle=$
$\displaystyle-\frac{{\tilde{\chi}}}{x}+\frac{{\mathcal{P}_{3}}}{x}\left[-\exp(-\sqrt{\mathcal{P}_{2}}x)I_{1}(x)+\cosh(\sqrt{\mathcal{P}_{2}}x)I_{2}(x)\right].$
(107)

where

$\displaystyle I_{1}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{x}\left\\{\sinh(\sqrt{\mathcal{P}_{2}}y)\delta(y)y\right\\}dy,$
(108) $\displaystyle I_{2}(x)$ $\displaystyle=$
$\displaystyle\int_{x}^{\infty}\left\\{\exp(-\sqrt{\mathcal{P}_{2}}y)\delta(y)y\right\\}dy.$
(109)