The Shape Dependence of Vainshtein Screening in the Cosmic Matter Bispectrum

TL;DR

This paper investigates how the shape-dependent Vainshtein screening affects the cosmic matter bispectrum in modified gravity theories, providing analytic expressions and quantifying the effect's scale and magnitude.

Contribution

It derives the shape dependence signature of Vainshtein screening on the matter bispectrum using a second order perturbative analysis in a cubic Galileon model.

Findings

Shape dependence enters the bispectrum with a scale-factor dependence of a^{3/2}.

The effect magnitude today can reach up to 2%.

Deviations from ΛCDM linear growth rate can be up to 5%.

Abstract

One of the most pressing questions in modified gravity is how deviations from general relativity can manifest in upcoming galaxy surveys. This is especially relevant for theories exhibiting Vainshtein screening, where such deviations are efficiently suppressed within a (typically large) Vainshtein radius. However, Vainshtein screening is known to be shape dependent: it is most effective around spherical sources, weaker around cylindrical objects and completely absent for planar sources. The Cosmic Web therefore offers a testing ground, as it displays many shapes in the form of clusters, filaments and walls. In this work, we explicitly derive the signature of the shape dependence of Vainshtein screening on the matter bispectrum, by considering a cubic Galileon model with a conformal coupling to matter and a cosmological constant. We perform a second order perturbative analysis, deriving analytic, integral expressions for the bispectrum, which we integrate using hi_class. We find that the shape dependence of Vainshtein screening enters the bispectrum with a unique scale-factor dependence of $\propto a^{3/2}$. The magnitude of the effect today is up to 2 % for a model whose linear growth rate deviates up to 5 % from $\Lambda$CDM.

Tables

Table 1. Table 1 : Definition of the Galileon models we use. The value of C 3 subscript 𝐶 3 C_{3} is obtained by means of a shooting algorithm in order to match the required Ω π , eff , 0 subscript Ω 𝜋 eff 0 \Omega_{\pi,\text{eff},0} . Gal 1, …, Gal 3 test the effects of Ω π , eff , 0 subscript Ω 𝜋 eff 0 \Omega_{\pi,\text{eff},0} , and Gal 4 and Gal 5 enable us to study the influences of C 2 subscript 𝐶 2 C_{2} while keeping Ω π , eff , 0 subscript Ω 𝜋 eff 0 \Omega_{\pi,\text{eff},0} constant.

Model name	${\Omega }_{\pi ,\text{eff},0}$	$C_3\left[{\mathrm{Mpc}}^2\right]$	$C_2$
Gal 1	$-0.01$	$\mathrm{4.276\hspace{0.17em}747}\text{×}{10}^9$	$0$
Gal 2	$-0.02$	$\mathrm{1.063\hspace{0.17em}205}\text{×}{10}^9$	$0$
Gal 3	$-0.03$	$\mathrm{4.699\hspace{0.17em}956}\text{×}{10}^8$	$0$
Gal 4	$-0.01$	$\mathrm{6.127\hspace{0.17em}518}\text{×}{10}^9$	$-\mathrm{5.925\hspace{0.17em}926}$
Gal 5	$-0.01$	$\mathrm{1.918\hspace{0.17em}030}\text{×}{10}^9$	$\mathrm{5.925\hspace{0.17em}926}$
$\Lambda$CDM	$0$	$0$	$0$

Equations

$$S=\int {\mathrm{d}}^{4}x\left[M_p^2\sqrt{-g}\left(\frac{R}{2}-\Lambda +{\mathcal{L}}_{\pi }\right)+{\mathcal{L}}_m\left[\left(1+\pi \right)g_{\mu \nu },{\psi }_i\right]\right],$$

$${\mathcal{L}}_{\pi }=-\frac{C_2}{2}(\nabla \pi {)}^2-\frac{C_3}{2}(\square \pi )(\nabla \pi {)}^2,$$

$$\frac{{\mathrm{d}}^2\vec{x}}{\mathrm{d}{t}^2}=-\vec{\nabla }\varphi +{\vec{F}}_5,\text{with}{\vec{F}}_5:=-\frac{\vec{\nabla }\pi }{2(1+\pi )},$$

$$C_2\square \pi +C_3\left({\left(\square \pi \right)}^2-R_{\mu \nu }{\nabla }^{\mu }\pi {\nabla }^{\nu }\pi -\left({\nabla }_{\mu }{\nabla }_{\nu }\pi \right)\left({\nabla }^{\mu }{\nabla }^{\nu }\pi \right)\right)=\frac{3{\rho }_m}{2(1+\pi )}.$$

$$G_{\mu \nu }=M_p^{-2}\left(T_{\mu \nu }^m+T_{\mu \nu }^{\pi }\right)-\Lambda g_{\mu \nu },$$

$$\frac{T_{\mu \nu }^{\pi }}{M_p^2}$$

$${\nabla }^{\mu }\left(T_{\mu \nu }^m+T_{\mu \nu }^{\pi }\right)=0.$$

$$C_2\Delta \pi +C_3\left((\Delta \pi {)}^2-({\partial }_i{\partial }_j\pi )({\partial }^i{\partial }^j\pi )\right)=\frac{3{\rho }_m}{2},$$

$${\partial }_z\pi =\frac{3{\rho }_0{z}_0}{2C_2},\Rightarrow \frac{|{\vec{F}}_5|}{|{\vec{F}}_G|}=\frac{M_p^2}{2C_2}.$$

$${\partial }_r\pi =\frac{C_{2}r}{2C_3}\left(\sqrt{1+\frac{r_v^2}{r^2}}-1\right)\text{with}{r}_v:=\sqrt{\frac{3r_0^2{\rho }_0{C}_3}{C_2^2}}.$$

$$\frac{|{\vec{F}}_5|}{|{\vec{F}}_G|}=\frac{M_p^2}{C_2}\frac{r}{r_v}.$$

$${\partial }_r\pi =\frac{C_{2}r}{4C_3}\left(\sqrt{1+\frac{r_v^3}{r^3}}-1\right)\text{with}{r}_v:={\left(\frac{4r_0^3{\rho }_0{C}_3}{C_2^2}\right)}^{1/3}.$$

$$\frac{|{\vec{F}}_5|}{|{\vec{F}}_G|}=\frac{M_p^2}{C_2}\frac{r^{3/2}}{r_v^{3/2}}.$$

$$\mathrm{d}{s}^2=a(\tau {)}^2\left(-\mathrm{d}{\tau }^2+{\delta }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j\right).$$

$$\frac{{\mathcal{H}}^2}{a^2}={\rho }_m+{\rho }_{\pi }+\frac{{\Lambda }_c}{3},$$

$$\frac{1}{a^2}\left({\mathcal{H}}^2+2{\mathcal{H}}^{\prime }\right)=-p_{\pi }+\frac{{\Lambda }_c}{3},$$

$${\rho }_{\pi }$$

$$p_{\pi }$$

$$-\frac{C_2}{a^2}\left({\pi }^{\prime \prime }+2\mathcal{H}{\pi }^{\prime }\right)+\frac{3C_3}{a^4}{\pi }^{\prime }\left(2\mathcal{H}{\pi }^{\prime \prime }+{\mathcal{H}}^{\prime }{\pi }^{\prime }\right)=\frac{3{\rho }_m}{2(1+\pi )}.$$

$$-{\rho }_m^{\prime }-3\mathcal{H}{\rho }_m+\frac{{\rho }_m{\pi }^{\prime }}{2(1+\pi )}=0.$$

$${\pi }_{1/2}^{\prime }=\frac{C_2{a}^2}{6C_3\mathcal{H}}\left(1\pm \sqrt{1+{\lambda }_v(a)}\right),\text{with}{\lambda }_v(a):=18H{\rho }_{m}t\frac{C_3}{C_2^2},$$

$$\mathrm{d}{s}^2=a^2\left[-(1+2\psi )\mathrm{d}{\tau }^2+(1-2\varphi ){\delta }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j\right].$$

$${\rho }_m=\bar{\rho }(1+\delta ),u^{\mu }=\frac{1}{a}\left[{\delta }_0^{\mu }+v^{\mu }\right],v^{\mu }=(v^0,\nabla v{)}^{\top }.$$

$$\delta \pi ={\pi }^{(1)}+\frac{{\pi }^{(2)}}{2}+\dots$$

$${\delta }^{(1)\prime \prime }+\left(\mathcal{H}+\frac{{\pi }^{\prime }}{2(1+\pi )}\right){\delta }^{(1)\prime }={\alpha }_{\delta }(\tau ){\delta }^{(1)},$$

$${\alpha }_{\delta }$$

$$A(\tau )$$

$$\left|\frac{\square \delta \pi (\tau ,\vec{x})}{\square \pi (\tau )}\right|\ll 1$$

$$\frac{\Delta {\pi }^{(1)}}{({\pi }^{\prime \prime }+2\mathcal{H}{\pi }^{\prime })}=\frac{3\bar{\rho }{\alpha }_{\pi }{a}^2}{2A(\tau )({\pi }^{\prime \prime }+2\mathcal{H}{\pi }^{\prime })}{\delta }^{(1)}.$$

$$\Delta {\pi }^{(1)}=\frac{3\bar{\rho }{\alpha }_{\pi }}{2A(\tau )}{\delta }^{(1)}.$$

$$\lambda (\tau ):=\left|\frac{2A(\tau )({\pi }^{\prime \prime }+2\mathcal{H}{\pi }^{\prime })}{3\bar{\rho }{\alpha }_{\pi }{a}^2}\right|\gtrsim 1,$$

$${\delta }^{(2)\prime \prime }+\left(\mathcal{H}+\frac{{\pi }^{\prime }}{2(1+\pi )}\right){\delta }^{(2)\prime }={\alpha }_{\delta }(\tau ){\delta }^{(2)}+S^{(\delta )},$$

Full text

In this work, we explicitly derive the signature of the shape dependence of Vainshtein screening on the matter bispectrum, by considering a cubic Galileon model with a conformal coupling to matter and a cosmological constant. We perform a second order perturbative analysis, deriving analytic, integral expressions for the bispectrum, which we integrate using hi_class. We find that the shape dependence of Vainshtein screening enters the bispectrum with a unique scale-factor dependence of $\propto a^{3/2}$. The magnitude of the effect today is up to $2\text{\,}\mathrm{\char 37\relax}$ for a model whose linear growth rate deviates up to $5\text{\,}\mathrm{\char 37\relax}$ from $\Lambda$CDM.

1 Introduction

The ongoing search for an explanation of the value of the cosmological constant has motivated a remarkable amount of research in modified theories of gravity [1, 2], in particular around the introduction of additional light scalar degrees of freedom. Although no conclusive solution to the cosmological constant problem has been found, alternative theories of gravity nowadays have become a field of research in their own right. The question at hand is a much broader one: what is the fundamental nature of gravity? How could deviations from general relativity (GR) appear and be observed in current and future surveys?

An extension of GR via a light scalar field will introduce a long-range fifth force whose strength could be similar in magnitude to the gravitational force. Since constraints on fifth forces are tight on Solar System scales [3], any significant modification of gravity must possess a screening mechanism suppressing the fifth force in the local environment. This can be obtained if non-linearities in the theory become important in specific regimes (for instance high density regions) whilst leaving the force unscreened on large scales.

One of the main contenders for such a screening mechanism is Vainshtein screening [4, 5] which suppresses the fifth force in regions of high second order derivatives of the scalar field. It was realized in Refs. [6, 7] that the effectiveness of Vainshtein screening is strongly dependent on the shape of the source mass. The screening is most effective around spherical sources, less effective around cylindrical bodies and non-existent in systems with planar symmetry.

The discovery of this shape-dependence motivates us to look for this effect in the cosmic web, where lots of different shapes are present, from clusters to filaments and walls. The dynamics of the cosmic web are investigated here using cosmological perturbation theory. Since Vainshtein screening is intrinsically non-linear, we have to go to at least second order in perturbation theory to observe any shape dependent effects: we will therefore use the matter bispectrum as our observable. The matter bispectrum is especially suitable to test for shape dependence as it is additionally sensitive to three wavenumbers, which form a closed triangle upon imposing background homogeneity. Different triangle shapes correspond to different shapes in real space, which can be more or less symmetric [8].

The simplest model exhibiting Vainshtein screening is the cubic Galileon [9, 10, 11]. Galileon fields are scalar fields that respect the Galilean shift symmetry $\partial_{\mu}\pi\rightarrow\partial_{\mu}\pi+a_{\mu}$ in flat space-time. Despite being characterised by higher derivative operators in the Lagrangian, the equation of motion of Galileon fields always remains second order, and therefore does not suffer from ghost instabilities. In four dimensions this Galileon model can be extended to include operators which have quartic and quintic dependence on the Galileon scalar field [9]. It is also possible to extend the theory further if we just require second order equations of motion, but do not insist on the symmetry, to the so-called Horndeski scalar-tensor theory [12, 13]. In this work we restrict ourselves to the cubic Galileon as the simplest theory to exhibit the shape dependent behaviour of interest.

The Galileon was originally introduced as a possible explanation for the current acceleration of the expansion of the universe without the need for a cosmological constant [9]. This was achieved through a mechanism known as self-acceleration [14, 15, 16]. Constraints on such self accelerating cosmologies come from observations of baryon acoustic oscillations, the cosmic microwave background, type 1a supernovae, the growth rate of structure, galaxy clustering and the integrated Sachs-Wolfe effect [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Constraints also come from bounds on the time variation of the gravitational constant in the solar system [30], and variations of particle masses and fundamental constants [31, 32]. For a general review we refer the reader to Refs.  [33, 2, 34, 35].

In the presence of the higher order Galileon operators, cosmological Galileon backgrounds also lead to differences between the speed of gravitational and electromagnetic waves [36, 37]. These theories are now very tightly constrained by the recent observation of a neutron star merger simultaneously through gravitational and electromagnetic waves [38, 39, 40, 41, 42]. Thus, we restrict ourselves to cubic Galileons in this work.

It was shown in Ref. [43] that the cubic Galileon is ruled out as the only energy component driving the late-time accelerated expansion of the Universe, as it predicts a negative integrated Sachs-Wolfe effect in conflict with current obervations. Therefore, we include a cosmological constant in our model and keep the Galileon energy density as a subdominant energy component throughout the entire evolution of the Universe. A cosmological model with such a sub-dominant cubic Galileon component was also considered in Ref. [44], where it was shown that independent initial perturbations in the Galileon field can break the usual correlation between density and velocity power spectra and lead to a form of stochastic bias.

The matter bispectrum was first studied for cubic Galileon models in Ref. [45] and, subsequently, generalised to Horndeski theories [46, 47, 48, 49]. These analyses did not include an explicit coupling to matter, only an indirect coupling due to the Galileon and matter fields coupling to the metric. Our analysis includes an explicit conformal coupling of the Galileon field to matter, which results in a stronger fifth force and stronger effects from screening.

This paper is structured as follows. We begin by introducing our Galileon model and shortly reviewing the effects of shape dependence in Section 2. Then, we derive analytic expressions for the matter bispectrum using cosmological perturbation theory up to second order in Section 3. In Section 4, we evaluate the obtained integral expressions for the matter bispectrum numerically using the latest version of the hi_class code [50, 51, 52], which is capable of integrating the exact Galileon equation of motion. We then conclude in Section 5. Some of the more involved formulas and definitions are postponed to the appendices in order to guarantee a better flow of reading.

2 Shape dependence of Vainshtein sceening

In this section we introduce our Galileon model and review the effect of shape-dependent Vainshtein screening. As a simple proxy for an extension of general relativity with Vainshtein screening we consider a cubic Galileon model similar to the one used in Ref. [45], but with a conformal coupling to the matter fields in order to have a stronger and more explicit fifth force. Furthermore, we include a cosmological constant $\Lambda$ since the cubic galileon cannot be the sole source of the late-time accelerated expansion of the Universe [43]:

[TABLE]

with the Lagrangian for the Galileon fields being:

[TABLE]

where $\Box\coloneqq\nabla_{\mu}\nabla^{\mu}$ with the covariant derivative $\nabla$. The exact form of the Lagrangian $\mathcal{L}_{m}$ of the matter fields $\psi_{i}$ is not important for the rest of the discussion; it is only relevant that the matter fields follow geodesics of the Jordan-frame metric $(1+\pi)g_{\mu\nu}$, which is called conformal coupling. In our conventions the field $\pi$ and the parameter $C_{2}$ are dimensionless and $C_{3}$ has inverse mass dimensions $2$. In the non-relativistic, weak-field limit the $(0,0)$-component of the Jordan-frame metric becomes $g_{00}=-(1+2\phi+\pi)$, with $\phi$ being the Newtonian gravitational potential. Particles following geodesics of the Jordan-frame metric are therefore subject to a fifth force mediated by the Galileon field:

[TABLE]

where $\vec{\nabla}$ is the flat-space gradient in the 3 spatial directions.

Since we will use the hi_class code for numerical analyses later on, we will use the hi_class normalization for densities, such that the Friedmann equation takes the form $H^{2}=\sum_{i}\rho_{i}$ (the same normalization will be used for pressures). For the energy-momentum tensor of the matter fields we make the standard cosmological assumption that matter is non-relativistic can thus be treated as a pressureless, perfect fluid with velocity field $u_{\mu}$. With the hi_class normalization of densities the energy-momentum tensor takes on the form $T^{m}_{\mu\nu}=3M^{2}_{p}\rho_{m}u_{\mu}u_{\nu}$.

Variation of the action Eq. (2.1) with respect to the Galileon field gives the field equation:

[TABLE]

Variation of the action with respect to the metric leads to the Einstein equations:

[TABLE]

with the energy-momentum tensor of the Galileon field being:

[TABLE]

Furthermore, the Bianchi identities give rise to the conservation equations:

[TABLE]

In order to demonstrate the shape dependence of the Vainshtein screening mechanism we examine Eq. (2.4) in the static and curvature free case and assume that $\pi\ll 1$:

[TABLE]

where $\Delta$ is the static spatial Laplacian. This equation can be solved for various different shapes of the source $\rho_{m}$ [7]:

• •

Planar Symmetry: For simplicity, we assume here that the source has the constant density $\rho_{0}$ if $z\leq z_{0}$ and zero otherwise. Any other configuration with planar symmetry qualitatively has the same result. We consider the fifth force outside of the source and its relative strength with respect to the gravitational force $\vec{F}_{G}$:

[TABLE]

For any configuration with planar symmetry the non-linear ($C_{3}$) term in Eq. (2.8) vanishes completely, there is no screening and the relative strength of the fifth force with respect to the gravitational force is constant in space.

• •

Cylindrical Symmetry: We assume the source to have constant density $\rho_{0}$ inside the radius $r_{0}$ (radius being defined as $r^{2}=x^{2}+y^{2}$ here) and zero otherwise:

[TABLE]

$r_{v}$ is called the Vainshtein radius and determines the scale below which the non-linear terms become important. Within the Vainshtein radius ($r_{0}<r\ll r_{v}$):

[TABLE]

We see that the fifth force becomes weaker than the graviational force inside the Vainshtein radius and is thus screened.

• •

Spherical Symmetry: The source is defined in the same way as for the cylindrical case except now radii are defined by $r^{2}=x^{2}+y^{2}+z^{2}$. We obtain:

[TABLE]

The relative strength of the fifth force inside the Vainshtein radius is now given by:

[TABLE]

We observe that the fifth force is screened more effectively inside the Vainshtein radius than for the cylindrical case.

To summarize, the more evenly the Galileon field depends on the three directions of space, the larger is the non-linear term in Eq. (2.8) and the more effective is Vainshtein screening. All of the above mentioned source symmetries appear on cosmological scales as walls, filaments and halos. Therefore, we will search for this shape-dependent effect in the distribution of cosmic structure.

More specifically, the shape-dependent term in Eq. (2.8) is non-linear, so that its effects can only manifest at second or higher order in cosmological perturbation theory. As a proxy for the shape dependence, we will therefore consider the matter bispectrum, which is sensitive to both non-linearities and shapes, due to its dependence on three wavevectors forming a closed triangle. Based on Ref. [7], we expect both the non-linearities and screening to be most effective for the equilateral configuration, which corresponds to the highest degree of symmetry, and least effective for the flattened configuration, where all three sides of the triangle are parallel thus containing only one-dimensional information in Fourier space.

3 Cosmological perturbation theory in the Einstein frame

In order to compute the bispectrum for our Galileon model, we have to perform a perturbative analysis up to second order in the matter density contrast. For this we start by analysing the background evolution and then proceed with a first and second order calculation. A similar analysis without an explicit coupling to matter was performed in Ref. [45].

3.1 Background evolution

At background level, the metric is given by the FLRW-metric. We choose to work with conformal time $\tau$:

[TABLE]

In the following, derivatives with respect to conformal time will be denoted by a prime. From this point on we will set $M_{p}=1$. Defining the conformal Hubble function $\mathcal{H}=a^{\prime}/a$, the Einstein equations become:

[TABLE]

where the background Galileon density and pressure are defined by:

[TABLE]

The background Galileon equation of motion takes on the form:

[TABLE]

Finally, the 0-component of the conservation equations becomes:

[TABLE]

3.2 Cosmological Vainshtein screening

Before we solve the Galileon field equation (3.6) numerically with hi_class, we can obtain useful analytical understanding by applying the simplifying assumption $\pi\ll 1$ which is typically true if the Galileon density (3.4) is subdominant on the background. This assumption enables us to approximate the continuity equation (3.7) by $\rho_{m}^{\prime}=-3\mathcal{H}\rho_{m}$ if we additionally assume the quasi-static approximation $\pi^{\prime}\sim\mathcal{H}\pi$. Both assumptions have been checked numerically for all the models we study in Section 4, and the quasi static approximation has been previously shown to hold for theories with Vainshtein screening in Ref. [53]. Under these assumptions the continuity equation has the solution $\rho_{m}=\rho_{m,0}a^{-3}$. It is now straightforward to solve the Galileon equation of motion:

[TABLE]

where $t$ is the physical time. Variants of this solution, which assumed $C_{2}$ to be negligibly small, have already been derived in Refs. [14, 44]. By solving the equation of motion numerically using hi_class we found for all models studied in Section 4 that only the negative branch of the solution is stable on a cosmological background.

The solution (3.8) is called the cosmological Vainshtein solution, because it exhibits a time-like Vainshtein screening effect. The function $\lambda_{v}$ measures the magnitude of the non-linear terms in the equation of motion against the linear terms and is the time-like equivalent of $(r_{v}/r)^{3}$ in the spatial solution around a spherical source (2.12). We will show in the next Section that the linearised fifth force will be screened at early times where $\lambda_{v}\gg 1$ compared to a situation where $\lambda_{v}\ll 1$.

3.3 Linear Perturbation Theory

We now look at scalar linear perturbations in the cosmic fluid. We choose Newtonian gauge such that the perturbed metric can be written in the following form:

[TABLE]

We parametrize the matter and velocity field perturbations by:

[TABLE]

The perturbed Galileon field is written as $\pi+\delta\pi$, where $\pi$ is the background field and $\delta\pi$ is the perturbative variable. The 0-component of the velocity, $v^{0}$, is determined by the constraint $u^{\mu}u_{\mu}=-1$. The perturbation variables $\phi$, $\delta$, $\delta\pi$, …are expanded into first order, second order, etc. perturbations, for example:

[TABLE]

The Galileon field, Einstein and Conservation equations (2.4), (2.5), (2.7) at linear order are listed in A. Combining these equations in Fourier-space using the quasi-static approximation, which we already used in the previous Section 3.2, and assuming subhorizon modes only ($\mathcal{H}^{2}\ll k^{2}$) we get a second order differential equation for the density contrast:

[TABLE]

where

[TABLE]

We see that the gravitational force of strength $M_{p}^{-2}=1$ is accompanied by a fifth force of strength $\alpha_{\pi}^{2}/(2A)$. The first term in $\alpha_{\pi}$ reflects the conformal coupling between matter and the Galileon field, and the second term is a consequence of matter coupling to the Galileon field indirectly with gravity as mediator. Inserting the cosmological Vainshtein solution (3.8) into the expression for $A(\tau)$ enables us to observe the cosmological screening of the linearised fifth force. Assuming matter domination for simplicity (radiation domination just gives different numerical factors), $A(\tau)=C_{2}(1+\lambda_{v}(a)/3)\approx C_{2}$ if $\lambda_{v}(a)\ll 1$ and $A=C_{2}(1+2\sqrt{\lambda_{v}}/3)\approx 2C_{2}\sqrt{\lambda_{v}}/3$ if $\lambda_{v}\gg 1$. We conclude that for $\lambda_{v}\ll 1$ the fifth force is unscreened with strength $\sim C_{2}^{-1}$ relative to the gravitational force which is analogous to the result (2.9), and for $\lambda_{v}\gg 1$ the fifth force is screened with strength $(C_{2}\sqrt{\lambda_{v}})^{-1}$ relative to the gravitational force, with $\lambda_{v}$ being the equivalent of $(r_{v}/r)^{3}$.

3.4 Breakdown of Perturbation Theory

In general relativity, cosmological perturbation theory is restricted to scales that obey $\delta\ll 1$. Below the non-linearity scale, the hierarchy between background, linear and non-linear terms in the Continuity and Euler equations breaks down. Similarly, for the background Galileon equation of motion (3.6) we had to assume $\left|\Box\delta\pi(\tau,\vec{x})\right|\ll\left|\Box\pi(\tau)\right|$, i.e. we neglected second derivatives of the spatially dependent part $\delta\pi$ of the Galileon field compared to second derivatives of the background field $\pi$ . This means that there is a second non-linearity scale in our system, the scale where $\left|\Box\delta\pi(\tau,\vec{x})\right|\sim\left|\Box\pi(\tau)\right|$ and a perturbative analysis breaks down.

A similar assumption had to be made when we derived the linear equation of motion (A.5). Additionally to assuming $\delta\pi\approx\pi^{(1)}$, we had to assume that terms non-linear in $\delta\pi$, like $C_{3}(\Box\delta\pi)^{2}$, are small compared to the linear term $C_{3}\Box\pi(\tau)\Box\delta\pi$. This makes the connection between this new linearity scale and the Vainhstein radius clear as both indicate the scale at which the non-linearities of the spatially dependent Galileon field become relevant in the equation of motion.

If this new non-linearity scale is larger than the standard non-linearity scale defined by $\delta\sim 1$, our analysis is even more restricted than a conventional perturbative analysis without the Galileon field. Therefore, it is important to check how the assumption

[TABLE]

compares to the assumption $\delta\ll 1$. The spatially dependent field $\delta\pi$ on the left-hand side of the condition (3.15) is in principle given by the full series $\delta\pi=\sum_{n=1}^{\infty}\pi^{(n)}/n!$ making the condition a non-perturbative statement about the perturbativity of the system. However, this makes it unfeasible to check the condition in practice. We therefore assume in the following that the order of magnitude of the non-linearity scale can be estimated if we simplify $\delta\pi\approx\pi^{(1)}$. This simplification is true on all scales up to the non-linearity scale, so the estimate for the non-linearity scale obtained in this way should be of the correct order of magnitude. Furthermore, we assume that the quasi-static approximation is valid, i.e. we can neglect the time derivatives in $\Box\delta\pi$ compared to the spatial derivatives. This assumption is valid as long as the non-linearity scale is a sub-horizon scale. With these assumptions the ratio on the left-hand side of the condition (3.15) is given by:

[TABLE]

We used that, on a linear level, the Einstein and Galileon field equations (A.3) and (A.5) can be combined to:

[TABLE]

We again made use of the quasi-static assumption. This shows that as long as:

[TABLE]

the new non-linearity scale leads to no additional restrictions of the perturbative analysis.

We will check this condition for every model we study in Section 4 and plot the quantity $\lambda$ as a function of time in Figure 2. In a simplified, analytically tractable setting, where we assume matter domination together with $C_{2}=0$, we find $\lambda=2$, i.e. the two non-linearity scales are of the same order of magnitude. What appears to be a coincidence at first sight is actually a consequence of the particular Galileon model we chose. The conformal coupling of the Galileon field with matter causes the background equation of motion to require $C_{3}(\Box\pi)^{2}\sim\bar{\rho}$ and the linear equation to enforce $C_{3}\Box\pi\Box\delta\pi\sim\bar{\rho}\delta$. Combining these two scaling relations, results in $\left|\frac{\Box\delta\pi(\tau,\vec{x})}{\Box\pi(\tau)}\right|\sim\delta$. Since it is not trivial that the two non-linearity scales will be equivalent for any Vainshtein screening model, we propose that the condition (3.15) should be checked for any perturbative study of Vainshtein screening.

3.5 Second Order Perturbation Theory

In order to compute the matter bispectrum we have to proceed to second order in perturbation theory which enables us to capture the onset of non-linear dynamics and the Vainshtein screening mechanism. Computing the Galileon field, Einstein and Conservation equations Eqs. (2.4), (2.5), (2.7) at second order and combining them appropriately, results in the inhomogeneous equation:

[TABLE]

where the inhomogeneity $S^{(\delta)}$ captures the non-linear physics. In B we state the source function $S^{(\delta)}$ in terms of the source functions $S^{(1)}$, $S^{(4)}$, $S^{(5)}$, $S^{(6)}$ and $S^{(7)}$ which were defined in the Appendix of Ref. [45] and are sufficiently long that we don’t reproduce them here. However, we would like to point out the appearance of one crucial term, which was also found in Ref. [46]:

[TABLE]

This terms has the same structure as the non-linear part of Eq. (2.8) and thus encodes the shape dependence of the Vainshtein screening mechanism.

Making use of the quasi-static approximation and a Fourier transform, $S^{(\delta)}$ can be brought into the form:

[TABLE]

where the kernel $\mathcal{K}$ is given by:

[TABLE]

with the linear growth rate $f=\mathrm{d}\ln D_{+}/\mathrm{d}\ln a$, where $D_{+}$ is the linear growth factor. We also made use of the first order result in Eq. (3.17) in order to express $\pi^{(1)}$ in terms of the density contrast $\delta^{(1)}$. The $k$-dependencies are captured by the form factors:

[TABLE]

The form factors $\alpha$ and $\beta$ are standard form factors appearing in cosmological perturbation theory in GR [54]. The additional form factor $\gamma$ originates directly as a Fourier transform of (3.20) and reflects very intuitively the shape dependence of Vainshtein screening. If the modes $\vec{k}_{1}$ and $\vec{k}_{2}$ are parallel, the term vanishes and no screening can occur. $\vec{k}_{1}$ and $\vec{k}_{2}$ being parallel means that we only capture one-dimensional information which is equivalent to a situation in real space with planar symmetry where, similarly no screening occurs, compare with Eq. (2.9). The appearance of $\gamma$ in Horndeski theories was already noted in Refs. [46, 49], but without making the connection to the shape dependence of Vainshtein screening.

Making the ansatz:

[TABLE]

we can solve the inhomogeneous differential equation (3.19) with Green’s method:

[TABLE]

The Green’s function $G$ is defined by:

[TABLE]

with the Wronskian $W$:

[TABLE]

The functions $D_{1}$ and $D_{2}$ are two independent solutions of the linear growth equation (3.12). In our numerical analysis in Section 4.2 we will associate $D_{1}$ and $D_{2}$ with the growing and decaying modes $D_{+}$ and $D_{-}$.

Assuming Gaussian initial conditions, the form factor $F_{2}$ is directly connected to the bispectrum:

[TABLE]

where $P(\tau,k)$ is the linearly evolved power spectrum. Therefore, the bispectrum depends on both the linear perturbations and the second order perturbations. It is thus convenient to introduce the reduced bispectrum:

[TABLE]

where the linear growth of the power spectra cancels out. Furthermore, the reduced bispectrum has the advantage of being mostly scale independent [54].

Using the kernel in Eq. (3.22), $F_{2}$ can be cast into the form:

[TABLE]

The time dependent functions $\mathcal{A}_{GR}$, $\mathcal{B}_{GR}$ and $\mathcal{B}_{\pi}$ are defined in A. While the functions $\mathcal{A}_{GR}$ and $\mathcal{B}_{GR}$ always appear in general relativity, $\mathcal{B}_{\pi}$ is a purely Galileon contribution describing the shape-dependent non-linearities in the Galileon equation of motion Eq. (2.8). We demonstrate in C that there is a relation between $\mathcal{A}_{GR}$ and $\mathcal{B}_{GR}$222We thank Emilio Bellini for pointing out this relation to us.:

[TABLE]

Defining $\mu$ as the cosine of the angle between $\vec{k}_{1}$ and $\vec{k}_{2}$, we thus conclude:

[TABLE]

where

[TABLE]

So we see that the shape-dependence of Vainshtein screening enters the matter bispectrum as a correction to the GR contribution $\mathcal{B}_{GR}$. The origin of $\mathcal{B}_{GR}$ can be traced back to non-linearities in the continuity and Euler equations where the flow of matter enters the total time derivative $\mathrm{d}/\mathrm{d}t=\partial/\partial t+\vec{v}\cdot\nabla$.

The contributions from Vainshtein screening to the form factor $F_{2}$, i.e. the term $\mathcal{B}_{\pi}(\tau)(1-\mu^{2})$, vanish in the flattened limit, where all the three wavevectors are approximately parallel and $\mu^{2}\approx 1$. In this case, the real-space, plane-wave density perturbations associated with these three wavevectors will only depend on one direction of space which is equivalent to the planar symmetry discussed in Section 2.

The nature of the signal imprinted by Vainshtein screening on the bispectrum may seem counter-intuitive: it vanishes where the fifth force is unscreened and is maximum where the screening is also largest. This is because to stronger screening correspond larger non-linearities, which are detected by the bispectrum. A situation of no screening corresponds to no additional non-linearities.

It is possible to compute the function $\mathcal{B}_{\pi}$ analytically in a simplified setting. During the rest of this Section we assume matter domination and that $C_{2}=0$. During matter domination, $D_{+}\propto a$ and $D_{-}\propto a^{-3/2}$ with $a=\rho_{m,0}\tau^{2}/4$, so the Green’s function takes on the form:

[TABLE]

The negative (stable) branch of the solution to the Galileon field equation Eq. (3.8) becomes:

[TABLE]

This enables us to compute:

[TABLE]

where we assumed $\pi\ll 1$. Using Eq. (3.4) we find $\rho_{\pi}/\rho_{m}=-\pi/2$ and we conclude:

[TABLE]

Finally, we can integrate for $\mathcal{B_{\pi}}$:

[TABLE]

The fraction $\rho_{\pi}/\rho_{m}$ is small due to the assumption of matter domination. During matter domination, $\mathcal{B}_{GR}$ has the standard value of $4/7$, see e.g. Ref. [54], with small corrections given by the impact of the Galileon field on the background evolution. These corrections are difficult to compute analytically but will be studied numerically in Section 4.3. We see that the contributions $\mathcal{B}_{\pi}$ from the Galileon field are small compared to $\mathcal{B}_{GR}$, but in contrast to $\mathcal{B}_{GR}$ evolve in time, which helps break the degeneracy between the two terms.

4 Numerical analysis with hi_class

In order to back up our analytic approximation in Eq. (3.40), and to generalise it beyond matter domination, we will evaluate the functions $\mathcal{B}_{GR}$ and $\mathcal{B}_{\pi}$ with hi_class333www.hiclass-code.net [50, 51], a Boltzmann solver for Horndeski-type models based on CLASS [52]. While previous public versions of hi_class [51] required one to parameterise the time evolution of the Horndeski $\alpha$ functions as defined in Ref. [55] in order to fully evolve the system, the latest version of the code [50] is able to integrate the full equation of motion of any Horndeski theory, including the Galileon.

Since hi_class works in the Jordan frame, we have to transform our action, Eq. (2.1), with the conformal transformation

[TABLE]

This transformation brings our model into the Jordan frame and allows it to be formulated in terms of the standard Horndeski functions $G_{i}(\pi,X)$ [12, 13]. The form of the Horndeski functions for our model is shown in A.

When working with both the Jordan and Einstein frame, one has to make sure to connect correctly between physical quantities in both frames. More specifically, we have to consider the density contrast in Einstein ($\delta$) and Jordan ($\tilde{\delta}$) frame which are related by:

[TABLE]

Using Eq. (3.17), we can show that:

[TABLE]

Therefore, on subhorizon scales and assuming the quasi-static approximation, the density contrast becomes approximately the same in both frames, see also Refs. [56, 57]. As we are only interested in the subhorizon scales where non-linear dynamics become important, the density contrast is effectively invariant under the conformal transformation and so is the matter bispectrum.

Other quantities like the Hubble function $\mathcal{H}$, the matter density $\rho_{m}$, the scale factor $a$, have to be transformed carefully under the conformal transformation. A summary of these transformations is given in A.

4.1 Background evolution

We now analyse the background evolution of the Universe in our Galileon model numerically using hi_class.

There is a subtlety with regards to the Galileon density $\rho_{\pi}$ when working in the Jordan frame, because the structure of the Friedmann equation (3.2) changes significantly. In the Jordan frame we have:

[TABLE]

This equation can be cast into the traditional form of a Friedmann equation, which is assumed by hi_class, by defining an effective Galileon energy density:

[TABLE]

where

[TABLE]

This effective Galileon density measures all deviations from a $\Lambda$CDM cosmology. Similarly, an effective Galileon pressure $p_{\pi,\text{eff}}$ can be defined.

The hi_class code checks for instabilities of the background by calculating the sign of the kinetic term and the sound speed of the scalar field. For all our models it turns out that stability is guaranteed if the Galileon field is negative. As a consequence, $\rho_{\pi,\text{eff}}$ is negative as well. However, the physical Galileon density, i.e. Eq. (3.4), or in terms of Jordan-frame quantities, the first line of Eq. (4.6), will always be positive.

In this analysis we will restrict ourselves to models, where $\Omega_{\pi,\text{eff}}\coloneqq\rho_{\pi,\text{eff}}/\rho_{\text{crit}}$ is small and the background evolution is at least roughly inside of current observational limits. Since the main goal of this analysis is to get a qualitative understanding of the shape-dependence effect, we postpone a thorough data analysis to future works.

We construct our background models with hi_class in the following way: Our Galileon model has two free parameters $C_{2}$ and $C_{3}$. Since we are mostly interested in the effect of shape dependence which is proportional to $C_{3}$, we consider $C_{2}$ to be either zero or subdominant on the background level. Thus, the value of $\Omega_{\pi,\text{eff},0}$ depends mostly on $C_{3}$. When we give hi_class a goal value for $\Omega_{\pi,\text{eff},0}$ as input, the code performs a shooting algorithm that fits the value of $C_{3}$ corresponding to the given $\Omega_{\pi,\text{eff},0}$. The parameter of $\Omega_{\Lambda}$ will always be used to fulfil the closure condition $1=\sum_{i}\Omega_{i,0}$. This also ensures that the sum of the energy densities driving the late time acceleration of the Universe $\Omega_{DE}\coloneqq\Omega_{\Lambda}+\Omega_{\pi,\text{eff}}$ will be identical to $\Omega_{\Lambda}$ in a purely $\Lambda$CDM model.

The Galileon models that we use throughout this work, labelled Gal 1-5, are defined in Table 1. The models Gal 1-3 enable us to study the effects of increasing $\Omega_{\pi,\text{eff},0}$, whereas the models Gal 4 and Gal 5 test the influences of the parameter $C_{2}$ while keeping $\Omega_{\pi,\text{eff},0}$ constant. We compare these models to $\Lambda$CDM, which corresponds to vanishing $\Omega_{\pi,\text{eff},0}$ achieved by the limit $C_{2}\rightarrow 0$, $C_{3}\rightarrow\infty$ (compare Eq. (3.8)). Although hi_class checks the stability of the considered models on the cosmological background, it is not guaranteed that these models will be stable on any background. However, for the sake of studying phenomenological apsects of Vainshtein screening, only stability on a cosmological background is paramount.

In order to check how much the inclusion of the Galileon field affects the background evolution, we consider an effective equation of state parameter $w_{\text{eff}}$ defined as:

[TABLE]

i.e. the equation of state parameter of all the energy components driving the late time acceleration of the Universe. The deviations of this quantity from $-1$ are plotted in Figure 1. The value of $w_{\text{eff}}$ is always smaller than $-1$ but tends towards $-1$ at late times. Deviations from $-1$ become large at early times, but dark energy only becomes dominant over the matter density at $z\approx 0.33$ for all considered models. Current bounds on the equation of state parameter of dark energy from the DES [58] indicate the values $w_{p}=-1.01^{+0.04}_{-0.04}$ and $w_{a}=-0.28^{+0.37}_{-0.48}$ for a parametrisation $w=w_{p}+w_{a}\,(a_{p}-a)$ with the pivot redshift $z_{p}=1/a_{p}-1=0.2$.

Before we compute the first and second order perturbations, we have to check the condition (3.18) for perturbativity of the Galileon equation of motion as outlined in Section 3.4. For this we plot the quantity $\lambda$ defined by Eq. (3.18) as a function of time in Figure 2. We find that the condition $\lambda\gtrsim 1$ is satisfied for all the considered models and for the entire evolution of the Universe, confirming the validity of our perturbative analysis. Figure 2 confirms our analytic prediction that $\lambda=2$ during matter domination. However, we also observe that $\lambda$ decreases over time in the late Universe, suggesting that perturbativity of the Galileon equation might be more restricted for a de-Sitter Universe. In particular, the value of $C_{2}$ appears to have significant impact on the late behaviour of $\lambda$. We leave a more detailed analysis of these phenomena for future work.

4.2 Linear growth

In this section we compute the linear growth rate numerically. In order to obtain the Green’s function in Eq. (3.28), we need two independent solutions of the linear growth equation (3.12) – let them be $D_{1}$ and $D_{2}$. To obtain them, we solve Eq. (3.12) for two different initial conditions. To establish an approximate connection between the solution $D_{1}$ and the growing mode and solution $D_{2}$ and the decaying mode, respectively, we set the initial conditions to be the solutions of the Meszaros equation [59, 60], valid during radiation and matter domination:

[TABLE]

where $y\coloneqq a/a_{eq}$ is the scale factor relative to the scale factor at radiation-matter equality $a_{eq}$. The thereby obtained solutions $D_{1}(a)$ and $D_{2}(a)$ are afterwards normalized such that $D_{1}(a=1)=1$.

The results of the integration for the model Gal 1 are shown in Figure 3. Since we assume that the Universe is dominated by matter and radiation when setting the initial conditions, $D_{2}$ deviates slightly from the true decaying mode, showing a final residual of $\sim 10^{-9}D_{1}(a=1)$. This small deviation is not a concern for the purposes of this work as all we need is the growing mode and another independent solution.

In order to quantify the linear growth in our models we consider the growth rate $f$. In Figure 4 we present deviations of $f$ from the $\Lambda$CDM result $f_{\Lambda\text{CDM}}$. We see that deviations do not exceed $5\text{\,}\mathrm{\char 37\relax}$. This is roughly within current observational bounds, which indicate order $10\text{\,}\mathrm{\char 37\relax}$ relative uncertainties for $f\sigma_{8}$ assuming a $\Lambda$CDM cosmology, where $\sigma_{8}$ is the amplitude of the power spectrum. See Ref. [61] for a compendium on past constraints on $f\sigma_{8}$ and Refs. [62, 63, 64, 65, 66, 67] for some recent developments.

4.3 The matter bispectrum

We can now integrate the functions $\mathcal{B}_{GR}$ and $\mathcal{B}_{\pi}$ in Eq. (A.7) and Eq. (A.8), where we recall that $\mathcal{B}_{\pi}$ and $\mathcal{B}_{GR}$ respectively describe the contributions to the bispectrum from the shape dependence of Vainshtein screening, and the non-linearities coming from the continuity and Euler equations.

In Figure 5, we show the relative difference between the sum of both contributions $\mathcal{B}=\mathcal{B}_{\pi}+\mathcal{B}_{GR}$ for the five Galileon models and $\Lambda$CDM. The relative difference scales as $\Omega_{\pi}\propto a^{3/2}$ during matter domination, in agreement with our analytical result derived in Eq. (3.40). The slope remains approximately the same at late times, but it appears to be sensitive to the value of $C_{2}$. In fact, we can see that the Galileon model Gal 5, characterised by $C_{2}>0$, displays a much shallower slope at $z<1$ compared to the other models.

In general, deviations in the bispectrum from $\Lambda$CDM are larger for models which also display significant modifications in the background evolution and linear growth rate, as illustrated by Figures 1 and 4. We observe the largest deviations in the bispectrum for the model Gal 3, where they are of order $2-3$\text{,}\mathrm{\char 37\relax}$$ at redshift $z=0$.

The deviations, with respect to $\Lambda$CDM, of the term $\mathcal{B}$ in the Galileon models have two different origins. First, the term $\mathcal{B}_{\pi}$ is altogether absent in $\Lambda$CDM, and second, the term $\mathcal{B}_{GR}$ depends on the evolution of linear perturbations (see Eq. (A.7) and Figure 4), which is also modified for the Galileons.

We compare these two contributions in Figure 6, to determine which is dominant. For the linear-evolution term, we display the difference between $\mathcal{B}_{GR}^{\rm(Gal)}$ and $\mathcal{B}_{GR}^{(\Lambda)}$, i.e. the term $\mathcal{B}_{GR}$ evaluated on a Galileon or a standard $\Lambda$CDM background. We can see that the shape-dependence in the Galileon equation of motion (2.8) is the dominant effect modifying the bispectrum compared to a $\Lambda$CDM cosmology.

The reason why $\mathcal{B}_{GR}^{\rm(Gal)}-\mathcal{B}_{GR}^{(\Lambda)}$ is subdominant can be understood from the definition of $\mathcal{B}_{GR}$, and generalises to other non-$\Lambda$CDM cosmologies. Eq. (A.7) shows that $\mathcal{B}_{GR}(\tau)$ is an integral over the time derivative of $D_{+}$, normalized by $D_{+}(\tau)$. This is quite intuitive since $\mathcal{B}_{GR}$ describes the non-linearities in the continuity and Euler equations due to the total time derivative $\mathrm{d}/\mathrm{d}t=\partial/\partial t+\vec{v}\cdot\nabla$ depending on the flow of matter. Since the first-order continuity equation in the quasi-static approximation reads $\Delta v^{(1)}=-\delta^{(1)\prime}$, velocities are related to the time derivative of the density fluctuations.

Since $D_{+}^{\prime}$ is normalized by $D_{+}(\tau)$ in the definition of $\mathcal{B}_{GR}(\tau)$ independent of the underlying cosmology, the difference between $D_{+}^{\rm(Gal)}$ and $D_{+}^{(\Lambda)}$ must be either growing at early times and decaying at late times or vice versa. This means that

[TABLE]

must change its sign at some time $\tilde{\tau}<\tau$. Consequently, $\mathcal{B}_{GR}^{\rm(Gal)}-\mathcal{B}_{GR}^{(\Lambda)}$ will always be small since negative and positive contributions cancel out when integrating. We conclude from this that, quite generally, intrinsically non-linear modifications of gravity, like the shape-dependent term in the Galileon equation of motion Eq. (2.8), will have a larger impact on the bispectrum than linear modifications which influence the bispectrum indirectly through $\mathcal{B}_{GR}$.

Finally, we consider the reduced bispectrum in Eq. (3.31), which is mostly independent of scale and linear growth. In Figure 7, we show the relative difference between the reduced bispectrum for the Galileon models and $\Lambda$CDM. The triangle $\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3}=0$ is parametrized by $\mu_{12}$, the cosine of the angle between $\vec{k}_{1}$ and $\vec{k}_{2}$, and the absolute values of $\vec{k}_{1}$ and $\vec{k}_{2}$. In both plots in Figure 7 we keep $k_{1}$ and $k_{2}$ constant while varying $\mu_{12}$; in the left panel ($a$), we set $k_{1}=k_{2}$, whereas in the right panel ($b$), we use $k_{1}=2\times k_{2}$.

In panel $(a)$, deviations from $\Lambda$CDM vanish for $\mu_{12}=1$. This corresponds to the flattened limit of the triangle with $k_{3}=2\times k_{1}$, shown in the right panel. This is in agreement with our expectations that planar symmetry, corresponding to the flattened limit, will result in no Vainshtein screening – see Section 3.5. The deviations from $\Lambda$CDM are maximal for $\mu_{12}=-0.5$, which corresponds to an equilateral triangle, i.e. the most symmetric configuration. The signal is only slightly smaller in the squeezed limit $\mu_{12}=-1$. In plot $(b)$, where $k_{1}=2\times k_{2}$, $\mu_{12}=-1$ and $\mu_{12}=1$ both correspond to the flattened limit and the signal vanishes. The signal is maximal for $\mu=-0.25$ which corresponds to $k_{3}=k_{1}$, i.e. an isosceles triangle, which, for $k_{1}=2\times k_{2}$ is the most symmetric configuration. In Refs. [46, 47] a qualitatively similar behaviour of the reduced bispectrum in general Horndeski theories was observed, but without making the connection to the shape dependence of Vainshtein screening.

Summarizing, the shape dependence of Vainshtein screening as seen by the bispectrum is perfectly analogous to the shape dependence in real space outlined in Section 2: in real space, the non-linearities responsible for Vainshtein screening are larger the more evenly the field depends on all three directions of space; in the bispectrum, the non-linearities are largest for the most symmetric triangle configurations.

5 Summary and Conclusion

In this work, we investigated the effects of the shape dependence of Vainshtein screening on the cosmic matter bispectrum. The dependence of Vainshtein screening on the shape of the source mass was first found in Ref. [7], which showed that more symmetric sources led to a more effective suppression of the fifth force. If Vainshtein screening is present in nature, we can then expect that it will leave an imprint on cosmic structures, given that the cosmic web is characterised by many differently shaped objects like walls, filaments and clusters.

The simplest modified gravity model displaying the Vainshtein screening mechanism is the cubic Galileon, which we used as a proxy to test the effects of shape-dependence on the matter bispectrum. Unlike previous work on the bispectrum in Galileon theories, we assumed the Galileon field to be conformally coupled to matter, so as to make the fifth force mediated by the Galileon field explicit. We also restricted ourselves to models where the Galileon energy density is subdominant throughout the entire evolution of the Universe, to reflect constraints ruling out the Galileon as the single component driving the accelerated expansion of the Universe.

Previous analyses [45, 46, 47] seemed to indicate that there was no qualitatively new information in the bispectrum that was not already present in the power spectrum, which is also easier to measure. However, in this work we have shown that a unique signature emerges in the bispectrum, imprinted by the shape dependence of Vainshtein screening, that would not be observable in the power spectrum alone. We performed an analytic, perturbative analysis of our coupled Galileon model in the Einstein frame, checking explicitly that the emergence of a new non-linearity scale associated to the Vainshtein radius does not lead to a breakdown of perturbation theory.

The non-linearities of Vainshtein screening leave an imprint on the form factor $F_{2}(\vec{k}_{i},\vec{k}_{j})$ of the reduced bispectrum, adding a time-dependent term $\mathcal{B}_{\pi}(\tau)(1-\mu^{2})$ varying with the cosine $\mu$ of the angle between two wavevectors $\vec{k}_{i},\vec{k}_{j}$. This $\mu$ dependence reflects the shape dependence found in [7]: in fact, we observe that the non-linearities are largest for the most symmetric configuration of the bispectrum triangle – i.e. the equilateral one – whereas they vanish for the flattened limit, which corresponds to planar symmetry in real space.

A simplified, analytic computation of $\mathcal{B}_{\pi}$ assuming matter domination revealed that $\mathcal{B}_{\pi}$ scales with the Galileon density, which was confirmed by numerical analysis using the hi_class code. The effect of the shape dependence on the matter bispectrum is found to be at percent level today for Galileon models where the Galileon energy density $\Omega_{\pi,\text{eff},0}$ is at percent level as well.

Critically, the shape dependence of Vainshtein screening displays a distinctive time dependence $\propto a^{3/2}$, which is dominant over corrections originating from modifications of the linear growth. This signature is an independent observational effect from modifications of the background evolution or the linear growth factor, which can be similar for different models of modified gravity. As such, it may be used to break degeneracies between theories of modified gravity predicting similar deviations in both.

Acknowledgments

We acknowledge use of the hi_class code [50, 51, 52]. The authors thank Miguel Zumalacarregui, Emilio Bellini and Ignacy Sawicki for generously sharing the latest version of hi_class that can evolve the exact Galileon equation of motion numerically prior to its official release. The authors also thank David Parkinson and David Seery for helpful discussions. This work was supported by a Research Leadership Award from the Leverhulme Trust. CB is also supported by a Royal Society University Research fellowship.

Appendix A Formulas and Definitions

The Conservation equations, Einstein equations and the Galileon field equation at linear order:

[TABLE]

[TABLE]

The functions $\mathcal{A}_{GR}$, $\mathcal{B}_{GR}$, $\mathcal{B}_{\pi}$ used to describe $F_{2}$ are defined by:

[TABLE]

With the transformation Eq. (4.1), the action Eq. (2.1) can brought into the form of a Horndeski theory defined by:

[TABLE]

with

[TABLE]

where $2X\coloneqq-(\nabla\pi)^{2}$, $G_{\mu\nu}$ is the Einstein tensor and the $G_{i}$ are arbitrary functions of $\pi$ and $X$. The Horndeski functions for our model are given by:

[TABLE]

Transformation rules between Einstein and Jordan frame (Jordan-frame quantities are marked by a tilde):

[TABLE]

Appendix B The source term $S^{(\delta)}$

Here we give the full expression for the source term $S^{(\delta)}$ in terms of the source terms $S^{(1)}$, $S^{(4)}$, $S^{(5)}$, $S^{(6)}$ and $S^{(7)}$ which can be found in the appendix of [45]. Since we have a Galileon model with a conformal coupling some of the source terms computed in [45] have to be expanded for our model:

[TABLE]

The source term $S^{(\delta)}$ can now be defined as:

[TABLE]

Appendix C Simplification of the form factor $F_{2}$

In this section we want to prove that $\mathcal{A}_{GR}+\mathcal{B}_{GR}=2$ which greatly simplifies the form factor $F_{2}$ in Eq. (3.32). For this, we firstly note that the following differential equation holds for the Wronskian $W$ defined in Eq. (3.29):

[TABLE]

Now we can compute $\mathcal{A}_{GR}+\mathcal{B}_{GR}$. For simplicity of notation we will not write all of the dependencies on the integration variable $\tilde{\tau}$ explicitly, however, we will denote dependencies if they differ from $\tilde{\tau}$ or are crucial for the understanding of the equations.

[TABLE]

From Eq. (C.3) to Eq. (C.3) we used that the linear growth equation Eq. (3.12) holds for the growth function $D_{+}$. Now we will integrate parts of this integral by parts:

[TABLE]

Substituting this result back into the full expression Eq. (C.3) for $\mathcal{A}_{GR}+\mathcal{B}_{GR}$ and using the equation Eq. (C.1), we arrive at:

[TABLE]