The splashback radius in symmetron gravity

TL;DR

This paper investigates how the splashback radius in dark matter haloes is affected by symmetron gravity, revealing potential observable deviations from General Relativity that could constrain modified gravity theories.

Contribution

It introduces a semi-analytical model of spherical collapse in symmetron gravity, showing how screening impacts the splashback radius and predicting measurable deviations.

Findings

Splashback radius can differ by up to 10% from GR predictions in symmetron models.

Screening effects become significant as overdensities grow, affecting the collapse dynamics.

Current experiments can potentially constrain symmetron parameters based on splashback measurements.

Abstract

The splashback radius $r_\mathrm{sp}$ has been identified in cosmological $N$-body simulations as an important scale associated with gravitational collapse and the phase-space distribution of recently accreted material. We employ a semi-analytical approach to study the spherical collapse of dark matter haloes in symmetron gravity and provide insights into how the phenomenology of splashback is affected. The symmetron is a scalar-tensor theory of gravity which exhibits a screening mechanism whereby higher-density regions are screened from the effects of a fifth force. In this model, we find that, as over-densities grow over cosmic time, the inner region becomes heavily screened. In particular, we identify a sector of the parameter space for which material currently sitting at $r_\mathrm{sp}$ has followed, during the collapse, the formation of this screened region. As a result, we find that for this part of the parameter space the splashback radius is maximally affected by the symmetron force and we predict changes in $r_\mathrm{sp}$ up to around $10\%$ compared to its General Relativity value. Because this margin is within the precision of present splashback experiments, we expect this feature to soon provide constraints for symmetron gravity on previously unexplored scales.

Figures

Equations

$$M(R,t)\propto a(t{)}^s,$$

$$\frac{4\pi }{3}R(t{)}^3{\rho }_b(t)={\left(\frac{4}{3\pi }\right)}^{2}M(R,t),$$

$$R(t)\propto a(t{)}^{1+s/3}.$$

$$\frac{d^{2}r}{d{t}^2}=-\frac{GM(r,t)}{r^2},$$

$$\xi =\frac{r}{r_{\ast }},\tau =\frac{t}{t_{\ast }};$$

$$M(r,t)=M(R,t)\mathcal{M}(r/R).$$

$$\Xi (\tau )={\tau }^{2/3+2s/9}.$$

$$\frac{d^2\xi }{d{\tau }^2}=-\frac{{\pi }^2}{8}\frac{{\tau }^{2s/3}}{{\xi }^2}\mathcal{M}\left(\frac{\xi }{\Xi (\tau )}\right),$$

$$\mathcal{M}(y)=\frac{2s}{3}{\int }_1^{\infty }\frac{d\tau }{{\tau }^{1+2s/3}}H\left(y-\frac{\xi (\tau )}{\Xi (\tau )}\right),$$

$$\rho (r,t)=\frac{1}{4\pi r^2}\frac{dM}{dr}(r,t),$$

$$S=S_{\phi }+S_M\left({\tilde{g}}_{\mu \nu },\Psi \right),$$

$$S_{\phi }=\int \sqrt{-g}{d}^{4}x\left[\frac{M_p^2}{2}R-\frac{1}{2}{\nabla }^{\mu }\phi {\nabla }_{\mu }\phi -V(\phi )\right],$$

$${\tilde{g}}_{\mu \nu }=A^2(\phi )g_{\mu \nu }.$$

$$\square \phi =V_{,\phi }-A^3(\phi )A_{,\phi }(\phi )\rho \equiv {\tilde{V}}_{,\phi }(\phi ),$$

$$F_{\phi }=-\nabla \log A(\phi ).$$

$$V(\phi )=-\frac{1}{2}\mu {\phi }^2+\frac{1}{4}\lambda {\phi }^4,$$

$$A(\phi )=1+\frac{1}{2}\frac{{\phi }^2}{M^2},$$

$$\tilde{V}(\phi )=\frac{1}{2}\left(\frac{\rho }{M^2}-{\mu }^2\right){\phi }^2+\frac{1}{4}\lambda {\phi }^4.$$

$$\rho >{\rho }_{ssb}\equiv M^2{\mu }^2$$

$${\rho }_b=\frac{1}{6\pi G{t}^2}\propto (1+z{)}^3.$$

$${\lambda }_0=\frac{1}{\sqrt{2}\mu },$$

$$\beta =\frac{{\phi }_0{M}_p}{M^2}.$$

$$F_{\phi }=-16\pi G{\beta }^2{\lambda }_0^2{\rho }_{ssb}\chi \nabla \chi .$$

$$\frac{d^2\chi }{d{r}^2}+\frac{2}{r}\frac{d\chi }{dr}=\frac{1}{2{\lambda }_0^2}\left[\left(\frac{\rho (r,t)}{{\rho }_{ssb}}-1\right)\chi +{\chi }^3\right].$$

$$\mathcal{L}[{\chi }_{i+1},{\chi }_{i-1};{\chi }_i]=0,$$

$$\mathcal{L}[{\chi }_{i+1},{\chi }_{i-1};{\chi }_i]\equiv {\mathcal{D}}_{\mathrm{K}}[{\chi }_{i+1},{\chi }_{i-1};{\chi }_i]-{\mathcal{D}}_{\mathrm{P}}[{\chi }_i,{\rho }_i]$$

$${\mathcal{D}}_{\mathrm{K}}\equiv \frac{{\chi }_{i+1}+{\chi }_{i-1}-2{\chi }_i}{h^2}+\frac{2}{r_i}\frac{{\chi }_{i+1}-{\chi }_{i-1}}{2h}$$

$${\mathcal{D}}_{\mathrm{P}}=\frac{1}{{\lambda }_0^2}\left(\left(\frac{{\rho }_i}{{\rho }_{\text{ssb}}}-1\right){\chi }_i+{\chi }_i^3\right).$$

$${\chi }^{\text{new}}(i)={\chi }_n(i)-{\frac{\mathcal{L}(\chi (i))}{\partial \mathcal{L}(\chi (i))/\partial \chi (i)}|}_{\chi (i)={\chi }_n(i)}.$$

$${\chi }_{n+1}(i)=\omega {\chi }^{\text{new}}+(1-\omega ){\chi }_n,$$

Full text

I Introduction

Gravity, one of the fundamental forces of nature, plays a crucial role in inferring our model of the cosmos as well as all the precision constraints placed on fundamental physics through cosmology. The theory of General Relativity (GR) introduced by Einstein a century ago Einstein (1916), provided a coherent theoretical framework within which to study all gravitational phenomena. While it is arguably one of the most successful theories of modern physics, having passed a host of empirical phenomena, there remain regimes of curvature and scale where GR has yet to be accurately tested. Its theoretical and phenomenological limitations are being fully explored, with an endeavour which is carried out at virtually all energy scales, ranging from the ultraviolet properties of the theory, down to energy scale of $H_{0}$, associated to the present-day expansion rate of the Universe Riess et al. (1998).

Upcoming large scale structure (LSS) surveys will provide unprecedented constraints on gravity on cosmological scales, allowing to discriminate among many theories alternative to GR. The phenomenology of theories of modified gravity (MG) on linear cosmological scales is fairly well understood, and it is commonly characterized in terms of modifications in the relation between matter density contrast and, respectively, the lensing and Newtonian potential Zhao et al. (2010); Ade et al. (2016); Pogosian and Silvestri (2016). On the other hand, it is well known that non-linear mechanisms in MG theories “screen away” the effects of additional degrees of freedom in high density regions. This ensures that any fifth force is suppressed and MG reduces to GR in regions where it has been tested with remarkable accuracy Will (1993).

A natural regime of interest is the intermediate range, between the screened and unscreened regimes, e.g. the regions of space at the boundaries of dark matter haloes. To this extent, a feature that is gaining prominence is the so-called splashback, which corresponds to an observable steepening of dark matter halo density profile close to the boundary Diemer and Kravtsov (2014). Locally, the position of this steepening contains interesting information about the clustering of dark matter shells and it can be understood as the dividing radius of single-stream and multistream sectors of the dark matter phase space. This feature has already been noticed in the self-similar spherical collapse framework developed and studied in Fillmore and Goldreich (1984); Bertschinger (1985), and generalized to 3D collapse in Lithwick and Dalal (2011). Self-similarity, however, is fully operational in a universe without a characteristic scale, such as the Einstein-de Sitter (EdS) universe with $\Omega_{m}=1$. Even though realistic applications of the same principle to $\Lambda$CDM universe are possible Shi (2016), in this paper we will focus on the collapse in EdS scenario and will leave more realistic scenarios for a future work.

The profiles of the largest dark matter haloes in the Universe, where galaxy clusters reside, can be mapped by measuring the deformation of background sources Kaiser and Squires (1993); Umetsu et al. (2011). This technique, known as lensing, has been used to measure the splashback feature around clusters Umetsu and Diemer (2017); Contigiani et al. (2018). It should be noted however that the most stringent constraints are obtained using the distribution of subhaloes traced by the cluster galaxy members More et al. (2016); Baxter et al. (2017); Chang et al. (2018); Shin et al. (2018). In this case the interpretation is nevertheless not straightforward and an accurate comparison with $N$-body $\Lambda$CDM simulations is required.

In this paper we consider the splashback radius in MG scenarios, investigating the microscopic effects of alternative theories of gravity on the dark matter shells accreting into the halo. Since we aim at gaining insight on the physical details, we do not resort to numerical simulations, but rather employ a semi-analytical method based on the framework of self-similar spherical collapse of Fillmore and Goldreich (1984). We focus on the class of theories of gravity that display the symmetron screening mechanism Hinterbichler and Khoury (2010). While we present an overview of the symmetron gravity in the main text, let us mention here that our analysis can be easily extended to other types of screening mechanisms, e.g. to Chameleon screening exhibited by $f(R)$ models Capozziello et al. (2003); Carroll et al. (2004), where the essence is the density dependence in the scalar field mass, rather than the field couplings.

We have organized our presentation as follows. In Section II we introduce the self-similar density profile and present the relevant equations of motion for the collapsing shells. Section III presents the basics of symmetron gravity. Section IV discusses our numerical methods and demonstrates the effect of the symmetron force on the phase space of the dark matter halo and the shift in the splashback radius.111In the interest of reproducibility we make our numerical codes available at http://github.com/contigiani/sym-splash. Finally, we discuss the implications of our findings and suggest potential further studies in Section V.

II Density profile

In order to study the motion of accreting material onto an overdensity, we first need to specify a matter density profile. In this work, we employ the so-called self-similar approximation to the problem of spherical collapse. In this context, the idea of self-similarity was introduced for the first time by Fillmore and Goldreich (1984), where it was shown that around EdS backgrounds, where the scale factor scales as a power-law of cosmic time, $a(t)\propto t^{2/3}$, the spherical collapse equations admit self-similar and self-consistent solutions.

Material surrounding a scale-free perturbation initially coupled to the Hubble flow eventually reaches turn-around and collapses onto a central overdensity. We denote by $R(t)$ and $M(r,t)$ the position of the turn-around radius at a time $t$ and the mass contained within the radius $r$, respectively. The mass within the turn-around radius can be written as a function of scale radius as:

[TABLE]

where the parameter $s$ is referred to as the accretion rate. In this model, $M(R,t)$ and $R(t)$ are related to each other through

[TABLE]

where $\rho_{b}(t)$ is the EdS background density at time $t$. This additionally implies that the position $R$ as a function of time also depends on $s$:

[TABLE]

Notice that $s$ and the mass of the present-day perturbation are the only free parameters of this model. In this work, we choose a fixed value of $s=1.5$ for the accretion rate, known to be representative for the low-redshift Universe in numerical simulations Correa et al. (2015); Diemer and Kravtsov (2014).

During spherical collapse, Gauss’s law ensures that the trajectory for each shell of material is influenced only by the mass contained within it. The equation of motion for each shell can be written as

[TABLE]

where the left-hand side is the Newtonian force $F_{N}(r)$ proportional to Newton’s gravitational constant $G$.

While before turn-around the mass within a shell is manifestly constant, afterwards this is not true: as multiple shells start orbiting the halo, their trajectories start intersecting. This phenomenon is known as shell-crossing and it is the principal reason why integrating Eq. 4 is not straightforward.

If we label each shell of material by its turn-around time $t_{\ast}$ and radius $r_{\ast}$, such that $R(t_{\ast})=r_{\ast}$, the trajectory for each shell is found to be independent from these quantities when self-similarity is satisfied. This can be verified by rewriting the equation of motion for the given shell in terms of the rescaled variables

[TABLE]

and by enforcing the mass profile $M(r)$ to be of the form:

[TABLE]

Notice that, from Eq. 3 it follows that the rescaling of the local turn-around radius $\Xi=\frac{R(t)}{r_{\ast}}$ can be also written as a function of $\tau$ alone:

[TABLE]

The system is then evolved through the following self-similarity equations for $\xi(\tau)$ and $\mathcal{M}\left(\xi/\Xi\right)$:

[TABLE]

[TABLE]

where $H(\dots)$ is the Heaviside step function, and the turn-around initial conditions for $\xi(\tau)$ are $\xi(\tau=1)=1$, $d\xi/d\tau(\tau=1)=0$. Notice that because these two equations are coupled to each other, they should be solved jointly to obtain self-consistent solutions for the orbits and the mass profile. This is done by starting from an initial guess for $\mathcal{M}(y)$ and then evaluating numerically the trajectories $\xi(\tau)$ using Eq. 8. The corresponding $\mathcal{M}(y)$, evaluated using Eq. 9, is then taken as an initial guess for the next iteration. This is repeated until convergence is reached and a final result for $M(r,t)$ is obtained. The corresponding density profile is then simply

[TABLE]

and it is shown in Figure 1. Notice in particular that its time-dependence is completely described by $\rho_{b}(t)$ and $R(t)$.

III Symmetron gravity

In this section we provide a brief overview of symmetron gravity and introduce the framework needed to study its effects on spherical collapse.

We consider a scalar-tensor theory of the form

[TABLE]

with

[TABLE]

$M_{p}$ being the Planck mass, and $S_{\mathrm{M}}$ the action for matter fields. The scalar field $\varphi$ couples to the Einstein frame metric $g_{\mu\nu}$ with Ricci scalar $R$, while matter fields (collectively represented by $\Psi$) couple to the Jordan frame metric $\tilde{g}_{\mu\nu}$. The two metrics are assumed to be related by the transformation

[TABLE]

Notice that such model is fully specified by the functions $A(\varphi)$ and $V(\varphi)$. Varying the action with respect to $\varphi$ gives us the equation of motion:

[TABLE]

where $\rho$ is the trace of the matter stress-energy tensor, equal to the local matter density, and $\tilde{V}(\varphi)$ is an effective potential. The fifth force per unit mass exerted by the field $\varphi$ and experienced by a matter test particle can then be written as:

[TABLE]

In this paper we will focus on a realization of such a theory, namely the symmetron model specified by the functions:

[TABLE]

[TABLE]

and effective potential:

[TABLE]

In this parametrization, the symmetron naturally assumes the form of an Effective Field Theory with ${\varphi\to-\varphi}$ symmetry.

In high density regions, where the condition

[TABLE]

is satisfied, the effective potential $\tilde{V}(\varphi)$ has only one minimum in $\varphi=0$ and the field is driven towards it, resulting in a null fifth force. In other words, high density regions are screened. In low density environments, on the other hand, the minimum is not located at zero. For example, for $\rho=0$ the vacuum expectation value is $\varphi_{0}=\mu/\sqrt{\lambda}$.

The fifth force can be constrained by local tests of gravity; to see in detail how local limits translate into bounds on the mass scale $M$ and the Mexican hat parameters $\mu,\lambda$ we refer the reader to Hinterbichler and Khoury (2010), for a general overview, and to the introduction of O’Hare and Burrage (2018), for a more recent analysis.

In an EdS background, the average matter density as a function of redshift $z$ is

[TABLE]

As the Universe expands, the symmetron can undergo spontaneous symmetry breaking (SSB) when ${\rho_{b}(z_{\mathrm{ssb}})=\rho_{\mathrm{ssb}}}$. For more details about the cosmological evolution of the symmetron field and the allowed expansion histories we refer the reader to Hinterbichler et al. (2011); Bamba et al. (2013). Let us stress however that we are not interested in the possibility of using the field $\varphi$ to drive the late-time expansion of the Universe, but we are only interested in the additional fifth force and its effects on spherical collapse.

In this paper we will work in terms of the dimensionless field $\chi=\varphi/\varphi_{0}$ and symmetron parameters composed by the average matter density at symmetry breaking $\rho_{ssb}$, the vacuum Compton wavelength

[TABLE]

and the dimensionless coupling

[TABLE]

Using these parameters, the fifth force sourced by the symmetron field can be written as:

[TABLE]

IV Spherical collapse with the symmetron

Having introduced the symmetron, let us now go back to the original goal of this paper, i.e. study spherical collapse in symmetron gravity with particular focus on splashback.

The splashback radius is commonly defined as the point where the density profile $\rho(r)$ is at its steepest. While this steepening is noteworthy because it can be detected as a departure from an equilibrium profile, this definition is clearly not suited for our study, where we assume a predefined density profile. Fortunately, the splashback radius is also known to be connected to the apocenter of recently accreted material and the location of the latest caustic visible in the density profile. Here we study the effects of the symmetron force on splashback by using this latter definition.

Our simulation is based on a system of equations that includes the spherical collapse equations, as discussed in Sec. II, coupled to the equation for the field profile of the symmetron field, discussed in Sec. III. We start by presenting our numerical method to compute both the symmetron field profile and the additional fifth force for the assumed density profile. We then proceed to integrate the shell equation to predict the fractional change in the splashback position in the presence of the symmetron force.

IV.1 Field profile

Assuming the temporal evolution of the field to be very fast compared to the other time-scales of the problem, i.e. the Hubble timescale and that of the clustering of matter, the dimensionless field profile $\chi(r)$ sourced by a density profile $\rho(r,t)$ satisfies the following equation:

[TABLE]

This quasi-static approximation is common in the literature Davis et al. (2012); Clampitt et al. (2012); Brax et al. (2012) and has been tested in the context of $N$-body simulations Llinares and Mota (2014); Noller et al. (2014). In order to provide a rough, order of magnitude justification for this assumption, let us just mention that the timescale associated to the field dynamics in vacuum is given by $\sim\lambda_{0}/c$. It is clear that in order for the symmetron field to be relevant for the dynamics of the spherical collapse, this $\lambda_{0}$ should be of the same order of magnitude as the scale of the cluster itself. The latter, of course, is several orders of magnitude smaller than $c/H_{0}$.

The static symmetron equation of motion (24) is a non-linear elliptical boundary value problem, for which we set the standard boundary conditions of vanishing spatial gradient of the field at $r=0$ and $r\rightarrow\infty$. We use a one-dimensional version of the Newton-Gauss-Seidel relaxation method for the numerical integration of the equation. This is a standard method used for obtaining the scalar field profiles in $N$-body simulations with modifications of gravity mentioned above.

In practice, we discretize our $1D$ static symmetron equation of motion on a regular grid of size $h$ and use a second order discretization scheme for all the derivatives.222We have tested some outputs of our integrator against the results of version where higher order discretization schemes are employed. For our particular problem we did not encounter significant differences in performance of the integrator and performed the main analysis with the version which employs the second order scheme. The resulting equation takes the form

[TABLE]

where

[TABLE]

contains the discretization of the Laplace operator

[TABLE]

and effective potential:

[TABLE]

The basic idea of the relaxation methods is to find a field profile from this equation which is closer to the real solution than a randomly chosen initial guess. This step is iterated over multiple (improved) guesses labelled $\chi_{n}(i)$ until convergence is reached.

At a given step we define an improved (new) field profile through

[TABLE]

Then we use a part of this new $\chi$ as the field profile for our next relaxation iteration:

[TABLE]

where $0<\omega\leqslant 1$ is a weight parameter with, in principle, a problem-dependent optimal value.

We employ two intuitive convergence diagnostics, where at each step we terminate the iteration if a certain parameter is within a predefined threshold. The first parameter is the residual function:

[TABLE]

and the second one is the all-mesh average of the fractional change in the field profile.

[TABLE]

To validate our integrator and convergence thresholds we compare the numerical solution to a known analytic solution. In our case, this known solution is an exact $\tanh(r)$ field profile, for which the corresponding density profile was recovered using Eq. 24.

When solving for the density profiles plotted in Fig. 1, we numerically evaluate the equation of motion in the range $[0,2]$ for $r/R(t)$, where the density profile for $r\geq R(t)$ is assumed to be constant. We make sure that the arbitrary choice of the upper limit has no effect on our results by testing larger values.

IV.2 Splashback

Once the symmetron field is found as a function of time, the present-day phase-space distribution of recently accreted material can be obtained by integrating numerically the equation of motion (4) with added fifth force (23) for different collapse times.

We find that after imposing self-similarity, the collapse equations can be written only as a function of three dimensionless symmetron parameters: the redshift of symmetry breaking $z_{\mathrm{ssb}}$, the dimensionless coupling $\beta$, and the ratio $\lambda_{0}/R(t_{0})$ between the vacuum Compton wavelength $\lambda_{0}$ and the present-day turn-around radius $R(t_{0})$. An important combination of these parameters is

[TABLE]

which explicitly sets the strength of the symmetron force according to Eq. 15.

From our testing, we found that values $\lambda_{0}/R(t_{0})\in[0.02,0.1]$ offer non-trivial cases. For $\lambda\sim R(t_{0})$ we always obtain thin-shell-like solutions, while for $\lambda\ll R(t_{0})$ the field is heavy and simply relaxes onto the minimum of the potential $\tilde{V}(\chi)$ in Eq. 18.

In Fig. 2 we illustrate our method and show how the symmetron force modifies the phase-space configuration of the latest accreted orbits (left-side plot). We find that the splashback position is significantly affected for parameter values $f\sim 1$, $z_{\mathrm{ssb}}\sim 2$ and $\lambda_{0}/R(t_{0})\sim 0.1$. These values imply $M\lesssim 10^{-3}M_{p}$, which is in agreement with local tests of gravity Hinterbichler and Khoury (2010).

From the same figure (right-side plot), it is clear that the innermost regions of the overdensity are screened from the effects of the fifth force at all times and this becomes relevant in the outer regions only for $z\ll z_{\mathrm{ssb}}$. Past this point, the force profile slowly transitions from a thick-shell to a thin-shell like behaviour, where the force gets progressively concentrated around the surface of the screened region Taddei et al. (2014). Due to the sudden drop in density associated with splashback, this surface is delimited by the splashback radius.

A systematic exploration of the symmetron effects on this feature as a function of all parameters is presented in Fig. 3, which represents our main result.

A clear trend with $z_{\mathrm{ssb}}$ is visible. Notice that the fractional change on the splashback position has an optimal peak as a function of $z_{\mathrm{ssb}}$ that is independent of $f$. If we call $z_{\mathrm{sp}}$ the accretion redshift of the shell currently sitting at the splashback position after its first pericenter, i.e. the splashback shell, we see that the effect is maximized when $z_{\mathrm{sp}}\simeq z_{\mathrm{ssb}}$. This is easily explained by studying the profile of the fifth force over time. For $z_{\mathrm{sp}}\gg z_{\mathrm{ssb}}$, the selected shell collapses when the symmetron is in its symmetric phase and the material spends the rest of its trajectory in a screened region, away from the effects of the fifth-force; for $z_{\mathrm{sp}}\ll z_{\mathrm{ssb}}$, the thin shell has had time to form before $z_{\mathrm{sp}}$ and the shell feels the effects of the fifth force only during a small fraction of its trajectory. Between these two limiting cases there is an efficient $z_{\mathrm{ssb}}$ for which the splashback shell has time to follow the formation of the thin shell and it is optimally positioned near the peak of the force profile for most of its trajectory. In our figure we show how this peak still has a dependence on $\lambda_{0}$, introduced by the presence of this factor on the symmetron equation of motion (24).

To conclude this section, we point out that the smoothness of the density profile as plotted in Fig. 1 has little impact on our results and no impact on the trends discussed above. Differences between the two prescriptions exist only for $\lambda_{0}\ll R(t_{0})$, when the field profile becomes susceptible to the small-scale features of the profile. However, since we expect the sharp caustic to be smoothed by gravitational instabilities, for our main results we chose not to use the discontinuous profile and assumed instead its smoothed version. Notice also that considering such high-resolution scenarios would introduce additional caveats (e.g. the presence of sub-structure) that are not the focus of this work.

V Discussion and conclusion

In this work we have explored how symmetron gravity affects the splashback feature at the edges of cosmological haloes. In our approach, we assume a self-similar mass distribution motivated by spherical collapse in an EdS Universe, where the shape of the spherically symmetric matter distribution is assumed to be only a function of $r/R(t)$. This allows us to easily solve for the corresponding symmetron fifth force and estimate its effects on the splashback feature by studying the changed phase-space distribution of recently accreted shells.

The main limitation of our study is the lack of a fully consistent framework where the density profile, the turn-around physics and the phase-space distribution are solved for in conjunction with the newly introduced symmetron equation of motion. As an example, we would expect a consistent framework to take into account the back-reaction of the scalar field on the density profile.

While deriving self-consistent solutions is outside the scope of this paper and more suited to $N$-body simulation studies, we find it useful to discuss the impact of our assumptions on the results. Changes to the turn-around physics are commonly studied through the use of different approximations, like a scale dependent Newton’s constant Schaefer and Koyama (2008); Brax et al. (2010); Hu et al. (2018); Nojiri et al. (2018); Lopes et al. (2018). In our case, if we maintain the assumptions of self-similarity and power-law accretion in Eq. 1, the main change to our formalism will come in the form of upgrading the numerical constant appearing in Eq. 2 to a function of the perturbation scale and cosmic time.

Previous works have estimated these corrections to be of the order of a few percentage points at $z\simeq 0$; see Taddei et al. (2014) for results in symmetron gravity and Lopes et al. (2018) for similar results in $f(R)$ theory. In particular, we expect our assumption to first break at a redshift $z$ such that the condition $F_{\varphi}(r)\sim F_{N}(r)$ is satisfied at the turnaround radius $r=R(t)$. In our analysis, however, we have seen that the effects on splashback are maximized when the collapse redshift of the splashback shell $z_{\mathrm{sp}}$ is equivalent to this transition reshift. After this point, the splashback shell is confined in the inner region and we expect its trajectory to be unaffected by the turn-around physics. Therefore, we consider our results around the peak of Fig. 3 to be robust against this assumption. For the same reason, however, we expect to lose predictability for higher values of $z_{\mathrm{ssb}}$, since the initial condition of the splashback shell will differ from what we have assumed.

Notice that the argument presented above also implies that our results can be extended to a standard $\Lambda$CDM scenario. The present-day splashback shell is expected to have collapsed in the matter-dominated era and to have followed a trajectory mostly unaffected by the late-time expansion, especially for low values of the accretion rate $s$ like the one considered here Shi (2016).

Effects of modified gravity on the structure of dark matter haloes are usually presented in the form of changes in the small-scale power spectra Cui et al. (2010); Davis et al. (2012); Li et al. (2013); Brax et al. (2012). In this analysis we focused instead on a particular scale, the splashback radius, and showed that up to a $10\%$ change can be induced (Fig. 3). It should be pointed out that Adhikari et al. (2018) was the first work to explore how modified gravity affects the splashback position. We stress, however, that our work differs from theirs in three major aspects. First, here we focus on symmetron gravity which displays a different screening mechanism from the chameleon or k-mouflage explored in Adhikari et al. (2018) . Second, while their results based on $N$-body simulations represent more realistic predictions, they do not allow for a simple exploration of the theory parameter space. Third, with our semi-analytical approach, we are able to gain insight by obtaining quantitative results as a function of multiple theory parameters and provide an explanation for the visible trends.

Observationally, splashback can be measured predominantly around galaxy clusters,å for which the present-day turn-around radius $R(t_{0})$ is of the order of a few Mpc. Our results therefore imply that this feature can be used to constrain fifth forces with vacuum Compton wavelength $\lambda_{0}$ just below the Mpc scale. Because measurements of splashback in the galaxy distribution around clusters have already achieved a precision below the size of our predicted effect More et al. (2016); Baxter et al. (2017); Chang et al. (2018); Shin et al. (2018), we expect to soon be able to constrain not only the symmetron, but other fifth force models on similar scales.

Note in particular that, while other works have explored the possibility of constraining symmetron gravity on Mpc scales Hammami and Mota (2017); Gronke et al. (2015), the range considered here for $\lambda_{0}$ is unconstrained for this model. Thus we expect a measurement based on splashback to naturally complement other results based on laboratory experiments Burrage et al. (2016); Brax et al. (2018), stellar and compact astrophysical objects Jain et al. (2013); Brax et al. (2014) or galactic disks and stellar clusters O’Hare and Burrage (2018); Llinares (2018); Desmond et al. (2018).

As the physics of splashback matures into a new cosmological observable, we expect it to play a powerful role on testing modifications of gravity, complementary to already established techniques such as those for large scale structure.

Acknowledgements.

OC and VV are supported by a de Sitter Fellowship of the Netherlands Organization for Scientific Research (NWO). AS acknowledges support from the NWO and the Dutch Ministry of Education, Culture and Science (OCW), and also from the D-ITP consortium, a program of the NWO that is funded by the OCW