{'fraction_non_alphanumeric': 0.03992372197638315, 'fraction_numerical': 0.019534019509571426, 'mean_word_length': 4.40629130319852, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical implementation of the Cubic Galileon model in PINOCCHIO

We present a perturbative treatment of nonlinear galaxy clustering in the context of the cubic Galileon modified gravity model, in terms of 2nd order Lagrangian Perturbation theory and an extension of ellipsoidal collapse that includes Vainshtein screening. We numerically implement such prescriptions in the approximate PINOCCHIO code, and use it to generate realisations of the matter density field and halo catalogues with different prescriptions for ellipsoidal collapse. We investigate the impact of three different approximations in the computation of collapse times on the halo mass function, halo bias and matter power spectrum. In the halo mass function, both the modified gravity effect and the screening effect are significant in the high mass end, similar to what is found for other MG models. We perform a comparison with N-body simulations to assess the validity of our approach, and show that we can reproduce the same trend observed in simulations for all quantities considered. With a simple modification to the grouping algorithm of PINOCCHIO to take into account the gravity modification, and without the need to re-calibrate the algorithm, we show that we can reproduce the linear halo bias and the mildly-nonlinear matter power spectrum of simulations with good accuracy, especially for the implementation with Vainshtein screening. We stress that, while approximate, our method is orders of magnitude faster than a full N-body simulation, making it an optimal tool for the quick generation of large sets of halo catalogues for cosmological observables.


INTRODUCTION
Within the standard ΛCDM cosmological model, the observed accelerated expansion of the Universe is ascribed to the presence of dark energy (DE) in the form of a cosmological constant Λ: a perfect fluid with constant energy and equation of state, whose negative pressure P Λ = −ρ Λ is responsible for the accelerated expansion. DE makes up most of the energy density of the Universe today together with cold dark matter (CDM), however, its nature remains elusive. In the standard framework, the cosmological constant Λ is thought to be related to vacuum energy: an interpretation that, despite the model fitting extremely well most cosmological observations (Alam et al. 2017;Aghanim et al. 2020;DES Collaboration et al. 2021), poses severe theoretical problems (Carroll 2001). Additionally, tensions in some of the cosmological parameters as measured by early and late time probes such as the Cosmic Microwave Background and redshift-space distortions of galaxies (Aghanim et al. 2020;Heymans et al. 2021;DES Collaboration et al. 2021) have prompted the study of structure formation in the context of beyond-ΛCDM models. Many alternatives have been proposed, from exotic DE models to modified gravity (MG, see Bull et al. (2016); Koyama (2018); E-mail: bhu@bnu.edu.cn Ishak (2019) for recent reviews). While the former are mainly characterized by a dynamically evolving equation of state, the latter focus on the possibility that General Relativity (GR) is not the correct theory to describe gravity on cosmological scales, and introduce an additional fifth force that drives cosmic acceleration.
On the other hand GR has successfully passed several stringent tests, therefore any viable MG model should be able to evade such constraints and reduce to GR on small scales. In order to be consistent with both small-and large-scale observations, the additional MG fifth force has to be shielded in high density regions. This is achieved by means of a so-called screening mechanism. Screening mechanisms can be classified in three general categories (Lombriser 2016), based on the condition that the Newtonian gravitational potential Ψ N or its derivatives exceed a certain threshold Λ T : (i) screening at large field values, such as the chameleon (Khoury & Weltman 2004) or symmetron (Hinterbichler & Khoury 2010) models. These are very similar to the large field inflation model: the screening effect turns on in regions where the Newtonian gravitational potential exceeds a certain threshold, |Ψ N | > Λ T . A high density environment plays a key role in this mechanism; (ii) screening with first derivatives, such as in k-mouflage models (Babichev et al. 2009). This screening effect operates when the local gravitational acceleration is above a certain threshold, |∇Ψ N | > Λ 2 T ; (iii) screening with second derivatives, such as Vainshtein screening (Vainshtein 1972). This screening mechanism activates when the local curvature is large, |∇ 2 Ψ N | > Λ 3 T . Unlike the chameleon mechanism, Vainshtein screening does not rely on the environment, the screening radius is jointly determined by the Schwarzschild radius of the object and the Hubble radius.
The investigation of alternatives to the standard cosmological model is indeed one of the key targets of modern cosmology. Upcoming Stage-IV galaxy surveys, such as Euclid 1 , LSST 2 , WFIRST 3 , DESI 4 , J-PAS 5 , CSST 6 , will probe the clustering of galaxies to high precision, allowing to test different MG models in the linear and mildly nonlinear regime of structure formation (Amendola et al. 2013;Alam et al. 2020).
It is therefore key that accurate theoretical modelling beyond GR is prepared in advance in order to compare to observations, with a particular focus on the nonlinear and mildly nonlinear regimes. The most reliable tool to trace the growth of structures deep into the nonlinear regime are N-body simulations. However, N-body simulations are computationally extremely expensive, in particular when they include MG.
This paper studies the modelling of nonlinearities in the framework of perturbation theory (PT) in the context of the cubic Galileon model (Nicolis et al. 2009). In particular, we focus on a numerical implementation in the PINOCCHIO code (PINpointing Orbit-Crossing Collapsed Hierarchical Objects, Monaco et al. (2002a,b); Taffoni et al. (2002); Monaco et al. (2013); Munari et al. (2017)), an algorithm to quickly generate simulated dark matter halo catalogues based on Lagrangian PT and ellipsoidal collapse.
We organize the papers as follows: in section 2 we describe the background, linear and nonlinear perturbation evolution of the cubic Galileon model. In section 3, we present the 1st and 2nd order Lagrangian PT for the cubic Galileon model. In section 4, we give the prescription of the "G3-PINOCCHIO" algorithm, including the extension of ellipsoidal collapse. In section 5, we show the nonlinear matter power spectrum and mass function obtained with the code, as well as a prediction for the linear bias. Finally, we present our conclusions in Section 6.

CUBIC GALILEON
Galileon gravity (Nicolis et al. 2009) is proposed in the inspiration of flat space quantum field theories. The latter are invariant under the transformation φ → φ + b µ x µ + c, where b µ and c are a constant vector and scalar in flat space respectively. The existence of b µ and c indicates respectively the Galilean and shift symmetry. To avoid the Ostrogradski instability, the field equation includes up to second order time derivatives. When extending to curved spacetime, the non-minimal coupling between the Galileon field and the metric breaks the Galilean symmetry, but keeps the shift symmetry preserved. Such model, called covariant Galileon model, has a La-grangian of the form where M pl is the Planck mass M −2 pl = 8πG, g, R, G µν are respectively the determinant of the metric, the Ricci scalar and the Einstein tensor, c 2 , c 3 , c 4 , c 5 are dimensionless constants and M 3 = M pl H 2 0 , with H 0 the present Hubble parameter. X = φ ;µ φ ;µ is the kinetic term of the scalar field. The semicolon in Eq. (1) represents the covariant derivative. There are three branches of covariant Galileon: the G3 branch, known as cubic Galileon, with c 3 = 0, c 4 = c 5 = 0, the G4 branch, also called quartic Galileon, with c 3 , c 4 = 0, c 5 = 0 and the G5 branch, also known as quintic Galileon, with c 3 , c 4 , c 5 = 0. Some recent studies that focus on constraining the covariant Galileon model by means of various cosmological observations include Barreira et al. (2012Barreira et al. ( , 2013bBarreira et al. ( , 2014b; Neveu et al. (2013); Peirone et al. (2018); Renk et al. (2017);. Although this model is disfavored by current CMB and galaxy clustering data, its extension can fit current cosmological observations (Peirone et al. 2019;. The breakthrough discovery of gravitational waves (GW) from the merger of a neutron star binary (Abbott et al. 2017a) by the LIGO and Virgo collaborations puts a tight constraint on the GW speed (Abbott et al. 2017b), effectively ruling out several MG models. In particular, the G4 and G5 branches are ruled out (Creminelli & Vernizzi 2017;Ezquiaga & Zumalacárregui 2017;Baker et al. 2017) with a high confidence level. In the G3 branch however, the scalar field is minimally coupled to gravity, and the GW speed is unaltered. Hence, G3 is still a viable MG model even according to GW costraints.
An additional reason to investigate the cubic Galileon model is related to the screening mechanism it features. Since the G3 field drives cosmic acceleration via the non-canonical kinetic energy, the extra gravitational force is screened via the Vainshtein mechanism. While chameleon screening has been extensively discussed during the past decades, in particular in the context of f (R) gravity (for example in Li et al. (2012); Puchwein et al. (2013); Llinares et al. (2014)), Vainshtein screening is usually studied as feature of nDGP modified gravity (e.g. Winther et al. (2015); Barreira et al. (2015); Winther et al. (2017); Hernández-Aguayo et al. (2021)). Studies of this type of screening mechanism that focus on cubic Galileon are instead more rare; some examples include the work of Schmidt (2009); Barreira et al. (2013a); Li et al. (2013).
A first step to include MG in the PINOCCHIO code focused indeed on f (R) gravity, which features scale-dependent growth and chameleon screening (Moretti et al. 2020). In this work, we present an implementation of scale-independent MG models with Vainshtein screening, focusing in particular on the G3 model, and study its performance with PINOCCHIO.

Background evolution under cubic Galileon
The Einstein equation and the scalar field equation can be obtained by varying the action with respect to the metric and scalar field respectively. In the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric, one can get the 1st-and 2nd-Friedmann where κ = 8πG,P r = 1 3ρ r . Quantities with a bar are background quantities and the over dot represents the time derivative. The expressions forP φ ,ρ φ arē The matter component satisfies the continuity equationρ m + 3Hρ m = 0, the radiation component satisfiesρ r + 4Hρ r = 0, and the scalar field component (or the DE component) satisfiesρ φ + 3H(ρ φ +P φ ) = 0. For the scalar field, we use the tracker solution of De Felice & Tsujikawa (2010): The meaning of tracker denotes that different initial conditions of the background Galileon field give rise to different time evolution that eventually merge into a common trajectory. Hereφ is the background field, H is the Hubble parameter and ξ is a dimensionless constant. Eq. (6) expresses the solution of the background scalar field in terms of the Hubble parameter. Naively, it would seem that there are three extra parameters in the G3 model compared with ΛCDM, namely {c 2 , c 3 , ξ }. However, that is not the case: in what follows we show that the number of parameters in G3 is actually the same as those in ΛCDM. By substituting the tracker solution Eq. (6) into the 1st Friedmann equation, one can get where E is the dimensionless Hubble parameter E = H/H 0 , Ω m,0 = ρ m,0 /(3M 2 pl H 2 0 ) and Ω r,0 =ρ r,0 /(3M 2 pl H 2 0 ) are the present matter density and radiation density respectively and a is the scale factor. Defining we can get a parameterised Hubble parameter To avoid the scaling degeneracy and without loss of generality (Barreira et al. 2014b), one can assume c 2 = −1. Combining the scalar field equation and the tracker solution, we can get a constraining equation on the G3 parameters: From Eq. (8) and (10), we have One can see that, once Ω Λ,0 is given, all the G3 model parameters are fixed and there are no additional parameters with respect to ΛCDM. Through the 2nd Friedmann equation, we obtain the acceleration parameteṙ (12) The above equations fully describe the background cosmology in the G3 model. Figure 1 shows the ratio of the Hubble parameter between G3 and ΛCDM as a function of the scale factor a: one can see that the two are the same at early times. Starting from a 0.1, the ratio decreases, reaches a minimum around a 0.5, then bounces back and approaches unity at present time.

Linear perturbations in cubic Galileon
In the context of linear perturbation theory and in the framework of GR, the matter anisotropic stress can be ignored. The Weyl potential Φ + = (Φ + Ψ)/2 felt by relativistic particles is equal to the Newtonian potential Ψ felt by non-relativistic particles, namely Φ + = Φ = Ψ. Here we denote the Newtonian potential as Ψ and the spatial curvature perturbation as Φ. Generally, the equality Φ = Ψ does not hold in MG models due to the existence of the Compton wavelength of the extra scalar field: the gravitational force below and above this wavelength is different. This phenomenon can be parameterised in the 00 component of the Einstein equation by means of a function µ L7 , that in general is time-and scale-dependent: where ∆ = δ + 3aHv/k, with δ = δ ρ m /ρ m the matter density contrast in the Newtonian conformal gauge, and v the irrotational part of the peculiar velocity. One can read the µ L function as the ratio between the effective gravitational coupling G eff and the Newton constant G N , with the superscript "L" denoting linear level. In general, chameleon models feature a k-dependent µ L function, which translates in scale dependent growth even at linear level. For the kessence type of MG models however, cosmic acceleration is driven by the non-canonical kinetic energy, the scalar field is effectively massless and the corresponding Compton wavelength is on or above the Hubble horizon scale. The effective gravitational coupling at linear level may not be equal to the Newton constant, but it is constant on all the scales relevant to the linear regime (Peirone et al. 2018). Hence, for this type of MG models (that include the G3 model we are considering in this work), Φ = Ψ is still valid and µ L is a function of time only.
The linear perturbation regime of MG models has been extensively studied in the past few years. In particular, it has been shown both theoretically and numerically that single field models can be re-expressed in the language of the Effective Field Theory of Dark Energy (EFTofDE, Gubitosi et al. (2013) (2020)). The action in the EFTofDE reads where δ g 00 , δ K µν , δ K, δ R (3) are perturbations of the time-time component of the metric, the external curvature and its trace, and the 3-dimensional Ricci scalar in the constant time hypersurface, S m is the minimally coupled term for all matter fields χ i with metric g µν , and M i ,M i , andM i are the EFT functions. The perturbation of the scalar field can be expressed by the infinitesimal time diffeomorphism, t → t + π(x µ ), where π is the perturbation of the scalar field. Following Pogosian & Silvestri (2016), we combine the Einstein equation and the scalar field equation under the quasi-static approximation, which applies only to scales below the sound horizon of the scalar field. Under this approximation, we ignore the time derivatives of the gravitational potential and the scalar field. One can then derive the following expression for µ L (a, k) in the EFTofDE framework: with In the upper panel of Figure 2, we show the linear part of the µ function as a function of redshift z and wavenumber k. The values are encoded in the color bar. On top of that, we show contours with constant µ L values. One can see that the contours are vertically distributed, in line with the fact that the modification of the gravitational force is scale independent. Furthermore, the µ values are significantly larger at low redshifts: the gravity enhancement reaches a value of around 2 at present time, while it approaches unity for redshifts z > 1. In the lower panel of Figure 2 we show the nonlinear µ NL function that includes Vainshtein screening, discussed in the next section, as a function of both redshift z (horizontal axis) and the local density contrast δ (vertical axis). One can see that with higher local density, the extra gravitational force is shielded for all redshifts considered.

Nonlinear clustering in cubic Galileon
We now consider the 2nd order perturbation contribution. Under the quasi-static approximation, after ignoring the terms suppressed by the time derivatives and Hubble expansion, the 00 component of the Einstein equation reads (Frusciante & Pace 2020) where ∂ denotes derivatives with respect to comoving spacial coordinates. The equation for the scalar field is Since Ψ = Φ is valid in both the linear and nonlinear regime, combining Eqs. (17) and (18) one gets where Assuming a spherically symmetric profile for δ φ , Eq. (19) becomes a second order ordinary differential equation. After integrating once, it becomes where m(r) = 4π r 0 r 2ρ m (r )δ m (r )dr is the mass enclosed in a sphere of radius r. Considering a top-hat density profile with radius R, when r ≤ R the physical solution of Eq. (21) reads ) where r V is the Vainshtein radius of the enclosed mass, given by where µ NL (a, R) includes the nonlinear contribution In the above equation, When the matter perturbations are large such that R R V , µ NL → 1 recovering GR, while for small perturbations when R R V , µ NL → µ L , recovering the linear result. In the lower panel of Figure 2, we plot µ NL as a function of redshift z and the top-hat object density δ . One can see that for low redshifts z < 1, it deviates from GR, same as the linear µ L function. When the density contrast is about 0.01∼0.1, µ NL is similar to µ L , roughly depending on time only.
For large values of δ > 1000, µ NL → 1, recovering GR. At fixed redshift, the value of µ NL for a high density contrast is smaller than that of a low density contrast.

LAGRANGIAN PERTURBATION THEORY FOR CUBIC GALILEON
In this section we briefly review the linear and 2nd order Lagrangian perturbation theory (LPT) for the G3 model Song & Hu (2021). The equation of motion for the displacement field ( S) in Lagrangian coordinates are Notice that the differential operation ∇ x is in Eulerian coordinates.
Transforming to Lagrangian coordinates one can write where J is the Jacobian of the coordinate transformation, and S i,i , S i, j are now spatial derivatives with respect to Lagrangian coordinates. Substituting Eq. (31) in Eq. (27), we can write the displacement field equation in Lagrangian coordinateŝ Expanding δ , δ φ , S with respect to a small parameter ε δ = εδ (1) + ε 2 δ (2) + ..., one can split the above equations into serial differential equations according to their orders. At 1st order, the displacement field equation reads Since (T − A(a)) only depends on time, S i,i (t, q) can be separated into time-and space-dependent components. Transforming the above equation to Fourier space, we obtain where the 1st order growth factor D 1 (t) satisfies Like in standard GR, the linear growth in G3 only has temporal dependence. For 2nd order in PT, the displacement field equation readŝ

T S
(2) Through Eq. (19) we obtain Plugging Eq. (40), (37) into Eq. (39), we can write The 2nd order displacement field then reads where D 2 ( k 1 , k 2 ,t) is a normalized second order growth factor that satisfies Since [T − A(a)] depends only on time, D 2 ( k 1 , k 2 ,t) can also be separated into the product of time component D 2 (t) and space compo- ]. Thus Eq. (42) can be rewritten as Redefining −2D 2 1 (t) D 2 (t) as D 2 (t), then the second order growth factor D 2 (t) satisfies Again, D 2 is independent of k, as it is in the standard GR case. To get the second order displacement field S (2) i,i ( k,t), one just needs to convolve the initial linear density fields with a kernel as in Eq. (44). The difference between G3 and ΛCDM lays in the temporal evolution of the linear and 2nd order growth factors. When A(a) = 4πGρ m , B(a) = 0,C(a) = 2πGρ m , the equations for the 1stand 2nd-order displacements recover ΛCDM (Song & Hu 2021).
We choose the initial conditions for the D 1 differential equation to match an Einstein-de Sitter universe, i.e. D 1 (a) = a, dD 1 da = 1. In the top panel of Fig. 3 we show D 1 (a) in G3 (red solid line) and ΛCDM (blue dashed line), with their ratio in the second panel, as a function of the scale factor a. One can see that D 1 in G3 is the same as in ΛCDM at early times. Starting from a 0.2, D 1 in G3 increases faster than in ΛCDM, with their ratio reaching ∼ 15% at present times. The initial conditions for D 2 are assumed to be D 2 (a) = 3 7 a 2 , dD 2 da = 6 7 a, corresponding to matter domination. In the third panel of Figure 3 we show D 2 (a) in G3 (red solid line) and ΛCDM (blue dashed line), with their ratio in the bottom panel, as a function again of the scale factor a. Similarly to D 1 , at early times D 2 for G3 and ΛCDM are the same. Moreover, at a 0.2, D 2 in G3 increases much faster than in ΛCDM. The difference is more sizable than in D 1 , about 60% for a = 1. Since the growth functions are enhanced, cubic Galileon predicts stronger galaxy clustering than ΛCDM. In Figure 4, we show the ratio of D 2 /D 2 1 . This ratio is exactly 3/7 in the Einstein-de Sitter universe. The extra clustering power in G3 comes from a larger effective gravitational constant G eff compared to the standard Newton constant G N .  2017) is a semi-analytical algorithm for generating realisations of hierarchical formation history of dark matter halos. It is based on LPT and ellipsoidal collapse. The code evolves the initial linear density field of a given primordial power spectrum on a regular grid in Fourier space. Then, particles are displaced according to Lagrangian PT. For cubic Galileon, the 1st and 2nd order initial displacement fields in Fourier space are given by The initial time a 0 is arbitrary and related to the normalisation choice for the linear power spectrum, PINOCCHIO chooses to normalise at a 0 = 1 so that D 1 (a 0 ) = 1. For convenience, one can define nor- .
Then the positions of dark matter particles then are given by where S (i) ( q, a) are the Fourier transforms of the displacements of Eq. 46. Ellipsoidal collapse is used to predict the collapsing time of each dark matter particle. The dynamics of ellipsoidal collapse can be completely described by following the evolution of the three principal axes of the ellipsoid (Bond & Myers 1996). The physical length of ellipsoid semi axes is r i = a i q, where q is the comoving radius of the "Lagrangian sphere" (a sphere concentric with and sharing the same mass as the ellipsoid, but with density equal to the background density) and a i denotes the time evolution of each axis. By enforcing mass conservation we can write a 3 q 3ρ m = a 1 a 2 a 3 q 3 ρ m,e , where ρ m,e is the density within the ellipsoid, so that the nonlinear overdensity δ is Following the approach of Bond & Myers (1996), we can derive the evolution for the a i for the G3 model: where and the deviation from ΛCDM is encoded in the µ function discussed in the previous sections. Eq. (51) is an integro-differential equation, which should in principle be solved for each particle, resulting in a computationally expensive prescription. In its standard version, PINOCCHIO relies on an approximation based on LPT to compute collapse times. Such an approximation however does not hold in the case of models where growth is substantially different than ΛCDM. For our purposes, we adopt the description of the ellipsoid's evolution presented in Nadkarni-Ghosh & Singhal (2016). This is equivalent to the Bond & Myers (1996) approach, but avoids the integral of Eq. 52 by resorting to a set of nine dimensionless parameters: Here λ a,i correspond to the eigenvalues of the deformation tensor and characterize the shape of the ellipsoid, λ v,i describe the deviation of the velocity of the ith axis from the background Hubble flow and λ d,i correspond to the eigenvalues of the tensor of second derivatives of the gravitational potential. Taking the time derivative of λ a,i , λ v,i , λ d,i , we obtain the second order ordinary differential equations where δ = λ d,1 + λ d,2 + λ d,3 . The initial conditions for this set of equations are i, j ( k, a 0 ), and D 1,ini is the normalized growth factor D 1 at an initial time a ∼ 10 −5 . At this initial time the Zel'dovich approximation is accurate enough, hence we can compute the 1st order initial displacement field S (1) i, j ( k, a 0 ) from the realisation of the density contrast. Then we diagonalise S (1) i, j ( k, a 0 ) to get its eigenvalues. Finally, by solving Eqs. (56) we can get the collapse time for each particle. The collapse time is defined as the moment of collapse of the first axis of the ellipsoid collapsed, λ a,i → 1.
Collapsed particles may become part of dark matter halos or of the filament network that connects them. Halos are constructed with an algorithm that mimics their hierarchical formation: (i) for each collapsing particle the code checks their 6 neighbours in Lagrangian space; (ii) a particle that collapses before its neighbours becomes a 1-particle halo; (iii) a collapsing particle may accrete on a halo if it "touches" it in Lagrangian space. The accretion condition is satisfied if, after the halo and the particle are moved using Eq. (49), their distance is less than a threshold distance (discussed below); (v) collapsed particles that don't accrete on halos become filaments, and may accrete later if a neighbour is accreted to a halo; (vi) halo mergers are checked each time a collapsing particle "touches" two halos, the merger takes place if the halos, moved again with Eq. (49), get nearer than a threshold distance.
In the PINOCCHIO code, the threshold distance for the collapsed particles being accreted into halos is defined as where R = (M h ) 1/3 , and M h is the halo mass. The threshold distance for merging between halos is defined as where R lar = (M h,lar ) 1/3 , with M h,lar is the mass of the larger halo, σ is the variance of the linear density contrast on the grid, and is a function of time. σ is a free parameter that controls the change of the "virial radius", and {e, f a , f m , f ra , f rm , f 200 } are additional free parameters that have been calibrated with N-body simulation in ΛCDM. Due to the enhancement of gravitational force in G3, the accretion and merging conditions are easier to be satisfied in G3, even for the same halo mass. Via the Poisson equation, one can absorb the gravitational enhancement into the rescaled density and eventually the mass: Then one can simply replace the effective masses computed with the above formula in the expression for R and R lar . On linear scales, µ L is independent of radius. From Eq. (61), one can read M h,eff = µ L (a)M h . On the nonlinear scale, µ NL (a, δ ) also depends on the radius or equivalently on the local density. For simplification, we fix δ = 200, then M h,eff = µ NL (a, 200)M h . As a consequence, the 5th force is screened in virialised objects. The output halo mass is computed as N × m p , where N is the number of particles within the halo and m p = 4.62 × 10 9 M /h is the particle mass.

RESULTS
We extend the standard PINOCCHIO algorithm according to what described in the previous sections. The code, dubbed as "G3-PINOCCHIO", will soon be available in the official repository. In order to compare our results with full N-body simulations, we choose the box size, number of particles and cosmological parameters as in Barreira et al. (2013aBarreira et al. ( , 2014a. We run "G3-PINOCCHIO" in a box with size 200 Mpc/h with 512 3 particles. The cosmological parameters we adopt are Ω m,0 = 0.279, Ω Λ0 = 0.721, h = 0.731, σ 8 = 0.997, n s = 0.953, both for G3 and ΛCDM. In what follows, we describe three different prescriptions for the G3 model, namely "lin-earG3", "vainG3" and "grG3". For all three G3 models, we compute the displacement fields of dark matter particles with 2LPT as presented in Section 3. For "linearG3", we compute the collapse times with the linear expression for the gravitational slip function, µ L (a), with no screening. For "vainG3", we use the Vainshtein screened gravitational slip function, µ NL (a, δ ). In "grG3" we compute ellipsoidal collapse as in GR: we choose this case as benchmark to demonstrate the effectiveness of the screening mechanism; in fact, the "grG3" case can be seen as an extreme case of screening. A similar approach is taken in Barreira et al. (2013aBarreira et al. ( , 2014a, where the authors define a "Linear model", equivalent to our "linearG3", and a "Full model", same as our "vainG3". However, their "QCDM model" differs from our "grG3": the former only modifies the background evolution, while in the latter we modify the 2LPT displacements. As a consequence, our "grG3" run already includes the power spectrum enhancement predicted by linear theory, i.e. the amplitude of the power spectrum on large scales is the same in "grG3" and "linearG3", "vainG3". In Fig. 5 we plot the distribution of collapse times for each particle in a smaller simulation with 200 3 particles and 500 Mpc/h box size. The upper panel shows the number of collapsed particles for each redshift bin, while the lower panel shows the cumulative numbers. The red, pink, cyan and purple shaded regions denote respectively ΛCDM, "grG3", "linearG3" and "vainG3" models. From the zoomin sub-panel, we see that at redshift z = 0 the cumulative collapsed particle numbers in ΛCDM and "vainG3" are almost the same. This is because we normalize the linear power spectra of both ΛCDM and G3 models with the same linear perturbation amplitude at z = 0. From the top panel of Fig. 5, one can see that the ΛCDM collapsing rate at low redshift is lower than the one of all G3 models. On the contrary, the ΛCDM model has more collapsed particles above redshift ∼ 1.8. This is because D 1,ini for ΛCDM is larger than that for G3 by 10 − 15% (as shown in Fig. 3), thus the ellipsoidal collapse for ΛCDM begins with larger density perturbations. This is again an artifact due to the chosen normalization setup. In the same figure, the "vainG3" collapsing rate is lower than the one of "linearG3" between redshifts 0.5 and 0. This is due to the µ NL going back to unity at late times, resulting in a weaker gravitational force with respect to the cases without screening. As an extreme screening limit, "grG3" further suppresses gravitational collapse.
The input linear power spectrum can be computed with an Einstein-Boltzmann solver such as CAMB/CLASS (Lewis et al. 2000;Lesgourgues 2011) or using the analytic fitting formula of Eisenstein & Hu (1998). To suppress sample variance in the power spectrum, the moduli of Fourier-space modes of the linear density field are fixed to their expectation value given by the power spectrum, while phases are randomly distributed between 0 and 2π. Since the linear growth rate in G3 is larger than the one of ΛCDM, to reach the same final amplitude the G3 model has to start from a more uniform density initial condition. For this reason, the high redshift matter power spectra in G3 are lower than those in ΛCDM cases.
In Fig. 6 we show the relative differences in the cumulative mass function of "linearG3" (green lines) and "vainG3" (blue lines) with respect to our benchmark, "grG3". The role of the Vainshtein mechanism can clearly be seen in the high mass end, where the fifth force is screened and the "vainG3" case shows a smaller deviation from "grG3" than "linearG3".
Additionally, we compute the linear halo bias. We use the PowerI4 package 8 to read the halo catalogues and compute the halo power spectra at different redshifts, then we compute the square root of the ratio between the halo power spectrum and the linear matter power spectrum for k < 0.2, namely the linear halo bias.
In Fig. 7 we compare our results for the halo bias with the N-body simulations of Barreira et al. (2014a). The cross symbols are our results, the solid curves are the best-fitting Sheth-Tormen formalism to the cumulative halo mass function data from simulations, and the triangles with error bars are the measurement of halo bias directly from the simulations. Different colors mark the different approximations for collapse times, as detailed in the legend. It can be seen that our results, despite being computed with a completely different method, match remarkably well the general trend of the N-body simulations.
Finally, we plot our result for the matter power spectrum in Fig.  8. For PINOCCHIO, this is obtained by displacing with 2LPT particles that do not belong to halos, while halo particles are distributed around the halo center of mass assuming an NFW profile (Navarro et al. 1997(Navarro et al. , 1996. These power spectra are computed with the Pylians package 9 . For comparison, we also show the simulation results with solid curves. The relative differences are with respect to the "grG3" model in our case, and with respect to the "QCDM" model in the case of the simulations. One can see that, at large scales, the two results are discrepant: in particular, the Barreira et al. (2013a) ratios also show the amplitude difference in the power 8 https://github.com/sefusatti/PowerI4 9 https://pylians3.readthedocs.io/en/master/ . Relative differences in the cumulative mass function with respect to the "grG3" case at a = 1 (top panel) and a = 0.8 (bottom panel) for our two approaches to ellipsoidal collapse, "linearG3" and "vainG3", shown with green and blue lines respectively. spectra that is due to the different linear evolution between G3 and ΛCDM, while ours do not. As mentioned previously, this is due to the fact that our "grG3" is not exactly the same as "QCDM", since we modify both the background and the Lagrangian perturbations up to the 2nd order while the "QCDM" model in Barreira et al. (2013a) only modifies the background evolution. For this reason, we cannot do a quantitative comparison with Barreira et al. (2013a), although it can be seen that our results reproduce the general trend of the simulations. We report in Appendix A a more detailed analysis of our halo catalogues. We compare the halo properties of matched and unmatched catalogues, including the halo mass distributions, the halo mass functions and halo power spectra.

CONCLUSIONS
This work presents 1st and 2nd order Lagrangian perturbation theory for the cubic Galileon model, as well as its ellipsoidal collapse dynamics. As in the case of ΛCDM, both the 1st and 2nd order growth factors of the cubic Galileon model are scale independent, since the extra scalar field which drives the accelerated expansion is effectively massless. Its Compton wavelength is indeed on the   , as a function of the halo mass. We mark our predictions with cross symbols, and compare to the results of Barreira et al. (2014a) showing their simulations results with triangles with errorbars, and their prediction based on the Sheth-Tormen formula with solid lines. We show the results for the "grG3" case, the "linearG3" case and the "vainG3" case with red, green and blue lines and symbols respectively.
Hubble horizon scale: below this scale, the modifications of gravity keep constant in space and vary slowly with time. Since both growth functions are larger that their ΛCDM counterparts, as shown in Fig. 3, clustering in enhanced. In particular, the maximum enhancement of gravity happens at redshift zero, when the effective gravitational constant is about twice the Newton constant. Once we include the Vainshtein screening mechanism in the high local curvature regime, the gravitational interaction restores to the standard GR case.
Using both the extension of standard 2LPT and ellipsoidal collapse to the G3 model, we create a new branch of the PINOCCHIO code, dubbed as G3-PINOCCHIO: our extension provides a fast tool to generate approximated dark matter halo catalogues with the G3 modified gravity model. We run the G3-PINOCCHIO code in a box with size 200 Mpc/h and 512 3 particles and study the properties of the obtained halo catalogues at different redshifts. To illustrate the effect of Vainshtein screening, we run three different implementations of ellipsoidal collapse, namely "linearG3", "vainG3" and "grG3". From these realisations we compute the cumulative mass  The green lines mark the ratio of the "linearG3" simulations to the "grG3" case, while we show the "vainG3" case in blue. Our results, measured from PINOCCHIO snapshots, are plotted as dashed curves. The shadowed areas mark the 1σ errors on the measurements. We also show the results from Barreira et al. (2013a) with solid lines, remarking that these ratios are obtained to their QCDM simulations, which does not correspond to our "grG3" case exactly since it only includes modifications to the background evolution.
function, the linear halo bias and the matter power spectrum. We find that, as for other types of modified gravity, the Vainshtein screening mechanism in cubic Galileon also suppresses the extra gravitational force in the high mass end of the mass function. We compute the linear halo bias and the matter power spectrum and perform a qualitative comparison to N-body simulations, showing that we can reproduce all trends remarkably well. Given the significantly reduced computational time required to generate halo catalogues with PINOCCHIO with respect to more computationally expensive N-body simulations, our implementation provides an optimal tool for the fast generation of large sets of realisations. Additionally, the code can be readily extended to include other MG models that feature Vainshtein screening.

APPENDIX A: HALO CATALOGUES
We describe here a more detailed analysis of the halo catalogues. Since we use the same initial conditions for the different implementations of the modified gravity model, namely "linearG3" and "vainG3", we can match the IDs of halos that are formed in both simulations. We call them "matched" halos, while halos that are only present in one of the two realisations are called "unmatched". In Fig. A1 we show the halo number density for the matched and unmatched halos. One can see that most halos appear in both realisations, while only ∼ 1% of the halos in "linearG3" are unmatched and ∼ 10% halos in the "vainG3" are unmatched. In Fig. A2 we show the mass distribution of the matched halos for the "linearG3" and "vainG3" prescriptions. Because the modification to the gravitational slip is stronger in the linear case with respect to the nonlinear "vainG3" one, the halo masses shift to the massive direction in the "linearG3" case (i.e. they mostly sit below the black solid line). In Fig. A3 we show the halo power spectra computed for the whole catalog, and for the catalogue split in matched and unmatched halos. The halo power spectra are computed with PowerI4, considering all halos with at least 100 particles. It can clearly be seen that the halo power spectra of the matched halos are very similar for both "lin-earG3" and 'vainG3", while they are significantly different for the unmatched halos. One can conclude that the differences in the final . We show results for the 'linearG3" and 'vainG3" cases, for the whole halo catalogue(marked as all, shown with green and blue lines) as well as for the catalogue split in matched (purple and cyan lines) and unmatched (red and orange lines) halos.
halo power spectra actually come from the unmatched halos. Finally, in Fig. A4 we show the relative differences in the halo power spectra with respect to the "grG3" model. At a = 1 we find a difference of ∼ 20% on all scales for "linearG3". Such difference is strongly reduced for "vainG3", highlighting once again the effectiveness of the Vainshtein mechanism in screening the MG fifth force. Similar results can also be seen for a = 0.8, with differences of ∼ 10% and ∼ 2% respectively for "linearG3" and "vainG3". This paper has been typeset from a T E X/L A T E X file prepared by the author.   )/P hh (k) grG3 a = 0.8,2LPT linearG3 vainG3 Figure A4. Relative differences in the halo power spectra at a = 1 and a = 0.8, shown for the "linearG3" prescription (green line) and the "vainG3" prescription (blue line). In both cases the ratio is computed with respect to the "grG3" case.