Proca-stinated Cosmology II: Matter, Halo, and Lensing Statistics in the vector Galileon

The generalised Proca (GP) theory is a modified gravity model in which the acceleration of the cosmic expansion rate can be explained by self interactions of a cosmological vector field. In this paper we study a particular sub-class of the GP theory, with up to cubic order Lagrangian, known as the cubic vector Galileon (cvG) model. This model is similar to the cubic scalar Galileon (csG) in many aspects, including a fifth force and the Vainshtein screening mechanism, but with the additional flexibility that the strength of the fifth force depends on an extra parameter -- interpolating between zero and the full strength of the csG model -- while the background expansion history is independent of this parameter. It offers an interesting alternative to LambdaCDM in explaining the cosmic acceleration, as well as a solution to the tension between early- and late-time measurements of the Hubble constant H_0. To identify the best ways to test this model, in this paper we conduct a comprehensive study of the phenomenology of this model in the nonlinear regime of large-scale structure formation, using a suite of N-body simulations run with the modified gravity code ECOSMOG. By inspecting thirteen statistics of the dark matter field, dark matter haloes and weak lensing maps, we find that the fifth force in this model can have particularly significant effects on the large-scale velocity field and lensing potential at late times, which suggest that redshift-space distortions and weak lensing can place strong constraints on it.

1 Introduction Understanding the laws of physics that govern cosmic structure formation is indispensable for probing into the true nature of gravity, because gravity is the dominant one of the four fundamental forces on cosmological scales. Ever since its establishment, General Relativity (GR) has been a cornerstone of modern cosmology. Even though the predictions of GR have been validated against many tests, these tests are usually limited to small scales such as the solar system [1], leaving the cosmological scales underexplored. The current observational results of these latter scales, which trace the dynamics of luminous and dark matter such as stars, galaxies, galaxy clusters, and extended filaments surrounding enormous voids, are generally in good agreement with the current concordance model of cosmology, ΛCDM, despite the fact that in recent years a number of tensions between the cosmological parameter estimates from different observational probes have emerged (e.g., [2][3][4][5][6]]). However, there is currently no compelling explanation of the smallness of the cosmological constant in this model, which is why alternative models to explain the cosmic acceleration, such as dynamical dark energy and modified gravity (MG), have been widely considered. In particular, in most alternative theories of gravity, the time evolution of large-scale structures can be significantly influenced, so that the observational data in cosmology may allow accurate tests of such models on large scales (for a recent review see [7]). The last decades have seen many attempts to modify GR. According to the Lovelock theorem, GR is the only theory with second-order local equations of motion for the metric field, which is derivable from a 4-dimensional action [7], and therefore modifications to GR often involve new dynamical the ones which are most sensitive to the effect of the fifth force in the GP theory. This paper is arranged as follows. In Section 2 we introduce the GP theory and the particular instances of it that we will focus on in this work. In Section 3 we describe the set up of the N -body and ray-tracing simulations on which all following results are based. This is followed by presentations of the main results of the dark matter field (Section 4), haloes (Section 5), and weak lensing (Section 6). Finally, we summarise and discuss in Section 7.

The Generalised Proca (GP) theory
In this work, we study the generalised Proca theory of gravity, the most general vector-tensor theories with second-order equations of motion, which contains Lagrangian operators up to cubic order of the Proca field. The action of this model is given by where g denotes the determinant of the metric tensor g µν , L m is the matter Lagrangian density, L F,2,3 are the Lagrangian operators introduced by the Proca field, A µ , and the last operator is the standard Einstein-Hilbert term with the Planck mass, M −2 Pl = 8πG, G is Newton's constant, and R is the Ricci scalar. The Proca field can be decomposed as where ϕ is the temporal component of the vector field, B i is its transverse mode which is divergence free, ∇ i B i = 0, and χ is the longitudinal scalar (which can also be referred to as the scalaron field) The matter Lagrangian density is related to the energy-momentum tensor of a perfect fluid as, which, assuming that matter is minimally coupled to gravity, satisfies the standard conservation equation ∇ µ T (m) µν = 0, (2.4) where ∇ µ denotes the covariant derivative compatible with g µν . Introducing the first derivative of the vector field, B µν = ∇ µ A ν , we can build the anti-symmetric Faraday tensor as F µν ≡ B µν − B νµ . The kinetic term of the Proca Lagrangian, L F , can be described as, 5) and the self-interaction terms of the vector field are given by, where X ≡ 1 2 g µν A µ A ν , b 2 ≡ m 2 is the mass-squared of the vector field that characterises the onset of the acceleration epoch, and b 3 , p 2 , p 3 are parameters of mass dimension zero in natural units. The choice is generic enough, leaving a viable parameter space in which the theory is free of ghost and Laplacian instabilities [22]. Importantly, due to the derivative self-interaction of the vector field in L 3 , the gravitational effect of the field can be screened in dense regions as required by solar system tests. The screening mechanism in this model is analogous to the Vainshtein mechanism [31]. In this work we set p 2 = p 3 = 1 as a working example to study the qualitative behaviour of the Proca field and refer to it from now on as cubic vector Galileon (cvG). With this choice, the GP theory behaves as the standard cubic scalar Galileon model (csG) in certain limits [29].
When deriving the equation of motions (EOM), we consider the perturbed Friedmann-Robertson-Walker metric in the Newtonian gauge where a(t) is the time-dependent scale factor which is normalised to a(t 0 ) = 1 at the present day, and δ ij = diag(+1, +1, +1) represents the spatial sector of the background metric that is taken here to be flat, k = 0. As shown in [29], we expect the 'back-reaction' of B i on the evolution of χ and Φ to be very small, justifying the neglect of the B i field in the simulations. To perform cosmological simulations for this model, we rewrite all required equations in ECOSMOG's code units, which we indicate as tilded quantities (details in [29]). The equations are then rescaled through χ = 3β sDGP 2βχ , (2.8) to make an educated choice of the cvG model parameter possible, by comparing it with the well studied sDGP model (for more details see [29]). To lighten our notation, we will drop the prime iñ χ . The modified Friedman equation, which depends on the EoM of ϕ at the background level, given by where H(a) is the Hubble expansion rate at a, H 0 = H(a = 1), Ω m is the matter density parameter, and Ω P the Proca field density parameters today, We have considered only non-relativistic matter; the inclusion of radiation and massive neutrinos is straightforward. Therefore, the background expansion history in this model is completely determined by H 0 and Ω m .
The modified Poisson equation, rescaled by Eq. (2.8), under the quasi-static approximation and in the weak-field limit takes the following form in code units, whereρ is the matter density in code unit, β sDGP is the coupling strength between matter and the brane-bending mode in the sDGP model, and α and β are two time-dependent functions given by and β(a) = 1 2 which are both fully fixed by specifying Ω m and the coupling constant b 3 redefined asβ 3 ≡ b 3 (8πGH 2 0 )/c 6 [29]. Finally, the EOM for the longitudinal mode of the Proca field, χ, in the weak-field limit and rescaled by Eq. (2.8), where the source term on the right-hand side is identical to that in the sDGP equation [32], and we have defined a new time-dependent function with the following dimensionless and time-dependent function Thus, given a matter density field, we can solve for the scalaron field χ from Eq. (2.15) and plug it into the modified Poisson equation Eq. (2.12) to solve forΦ. OnceΦ is at hand, we can use finitedifference to calculate the modified gravitational force, which determines how the particles move subsequently. Note that, in this model,Φ not only determines the geodesics of massive particles, but also those of massless particles such as photons -in other words, the lensing potential is also modified. Asβ 3 is the only 'free' parameter that enters in all three key equations, it is practical to use it as the model parameter.

Cosmological simulations
In this section we present the set of dark-matter-only simulations for five different cosmologies which we use to investigate the phenomenology of the cvG model. Four of these take different values of the model parameter of the cvG model,β 3 = [10 −6 , 10 0 , 10 1 , 10 2 ], and one is their QCDM counterpart 1 , in the simulation. It is equivalent to the limitβ 3 → ∞ [29]. To study the cvG effects on the weak lensing (WL) signal, we extended the N -body code developed in the previous work [29] by adding an independent set of ray-tracing modules taken from Ray-Ramses [30]. This allows us to calculate the WL signal 'on-the-fly' as proposed by [33,34], while taking full advantage of the time and spatial resolution available in the N -body simulation.
We construct a light-cone for each cosmology by tiling a set of five simulation boxes, all having an edge-length of L box = 500h −1 Mpc, as shown in Fig. 1. The simulations treats dark matter as collisionless particles described by a phase-space distribution function f (x, p, t) that satisfies the Vlasov equation where p = a 2 m 0 ∂x/∂t, m 0 is the particle mass, and Ψ is the modified Newtonian potential given by Eq. (2.12). Note that as we do not include matter species such as photons and neutrinos the two Bardeen potentials are equivalent, Ψ = Φ. Hence to solve Ψ, and prior to it the longitudinal Proca mode, via Eq. (2.15), they are discretised and evaluated on meshes using the nonlinear Gauss-Seidel relaxation method [32]. The domain grid -which is the coarsest uniform grid that covers the entire simulation box -consists of N grid = 512 3 cells, which is equal to the number of tracer particles, N p . ECOSMOG is based on the adaptive-mesh-refinement code RAMSES [35], which allows mesh cells in the domain grid to be hierarchically refined -split into 8 child cells -when some refinement criterion is satisfied. In our simulations, a cell is refined whenever the effective number of particles inside it exceeds 8. This gives a higher force resolution in dense non-linear regions, where the Vainshtein screening becomes important. The Gauss-Seidel algorithm is run until the difference of the two sides of the PDE, d h , is smaller than a predefined threshold . We verified that for a value of = 10 −9 > |d h |, the solution of the PDE no longer changes significantly when is further reduced. We use the same set of five different initial conditions (ICs), for each of the five simulations that make up a light-cone for a given cosmology are different, for the different cosmologies. The ICs were generated using 2LPTic [36], with cosmological parameters taken from the Planck Collaboration [37], h = 0.6774, Ω Λ = 0.6911, Ω m = 0.389, Ω B = 0.0223, σ 8 = 0.8159. (3. 2) The linear matter power spectrum used to generate the ICs is obtained with CAMB [38]. The simulation starts at a relatively low initial redshift z ini = 49, or a ini = 0.02, justifying the use of second-order Lagrangian perturbation theory codes such as 2LPTic. One possible concern may be that, at this scale factor, differences of matter clustering are already present. However, judging from our experience [29], at this time the difference between the growth factors of the cvG model with ΛCDMis well below sub-percent level, so that modified effects on the initial matter clustering can be neglected. The light-cone, outlined by solid blue lines in Fig. 1, is constructed by positioning the five simulation boxes, outlined by solid red lines in Fig. 1, relative w.r.t. the observer. The geometrical set-up was constructed to place the sources at z s = 1, which is the starting point when the growth rate of matter density perturbations becomes higher than in ΛCDM [29]. The field-of-view (FOV) is set to 10 × 10 deg 2 (so that the wide end of the light-cone is still narrow enough to fit in the simulation box), within which 2048×2048 rays are followed by Ray-Ramses to compute quantities of interest. Ray-Ramses is an on-the-fly ray-tracing code. The rays are initialised when a given simulated box reaches a defined redshift (for the closest and furthest box to the observer the initialisation redshift is respectively z i = 0.17 and z i = 1.0), and end after they have traveled the covered length of the box, meaning 500h −1 Mpc. As here we are interested in the lensing convergence, κ, the quantity that is  computed along the rays is the two-dimensional Laplacian of the lensing potential, where 1, 2 denote the two directions on the sky perpendicular to the line of sight (LOS). The values of these two-dimensional derivatives of Φ lens,2D can be obtained from its values at the centre of the AMR cells via finite differencing and some geometrical considerations (see Refs. [30,34]). Integrating this quantity as lens,2D (χ, β(χ))dχ, (3.4) where c is the speed of light, χ is the comoving distance, χ s the comoving distance to the lensing source, and β(χ) indicates that the integral is performed along the perturbed path of the photon (χ is not to be confused with the longitudinal mode of the Proca field). The integral is split into the contribution from each AMR cell that is crossed by a ray, which ensures that the ray integration takes full advantage of the (time and spatial) resolutions attained by the N -body run. For the WL signal we wish to study in this paper, we employ the Born approximation, in which the lensing signal is accumulated along unperturbed ray trajectories. We will make further notes on the calculations in Sec. 6.1.

Matter field statistics
In this section we present the results of various dark matter statistics of the different cvG models and compare them with the predictions by QCDM, to study the impact of the Proca field on these key observables. We start with an analysis of the power spectra in Sect. 4.1. In Section 4.2, we consider the leading non-Gaussian statistic in large-scale structure clustering, the bispectrum, which is thus sensitive to deviations from linear evolved perturbations from single field inflation.
To support the analysis and interpretation of the results, we will compare the results of the Nbody simulations to Eulerian standard perturbation theory (SPT), and limit the comparison only to the tree-level statistics. In SPT, the energy and momentum conservation equations can be solved order by order to obtain higher-order corrections to the quantities of interest. The expansion in powers of the linear density field is a simple time dependent scaling of the initial density field (in the Einstein de Sitter approximation), for which the n-th order solution is with the conformal time τ = dt/a, k 1...n ≡ k 1 + ... + k n , the density contrast δ = ρ/ρ, δ (D) the 3D Dirac delta function, and F n the SPT fundamental mode coupling kernel [40,41].
When comparing a cvG model to the QCDM counterpart, we do so through their relative difference which we write in short hand as with A a placeholder of the summary statistics, and X will be one of the four cvG models. We calculate ∆A/A QCDM for each of the five pairs of cvG and QCDM simulations that share the same initial conditions to find its average and 1σ uncertainty. Taking this ratio removes contributions from cosmic variance, and so its uncertainty is not a direct indicator of how sensitive the various summary statistics are to differences between the cvG models. To provide an estimate of this sensitivity given a survey volume as large as our simulation box, we calculate the signal-to-noise ratio (SNR) of the difference between cvG models and their QCDM counterpart for some summary statistics using the expression where ∆A is the average and σ is the standard deviation of the five simulations per cosmological model. However, we note that the SNR values obtained in this way are subject to sample noise, owing to the small number of realisations.

Matter and velocity power spectra
To gain insights into the differences of matter clustering and peculiar velocities on linear and nonlinear scales among the various models in this work, we begin our study of dark matter phenomenology by considering the auto power spectra of the matter over-density, δ, given by Cosmic structure formation is driven by the spatially fluctuating part of the gravitational potential, Φ(x, t), in Eq. (2.7), induced by the density fluctuation δ. In cvG cosmologies we expect an additional boost to the standard gravitational potential with respect to its QCDM counterpart, induced by χ described by Eq. (2.12), in regions where the fifth force is not screened by the Vainshtein mechanism. Thus, clustering will be enhanced in the cvG models on some scales, which can be captured by P δδ . The top row of Fig. 2 compares the linear matter power spectra (black dotted lines) with the simulation results of each cosmology (coloured lines with shaded regions), at a = 0.6 (outer left), a = 0.7 (inner left), a = 0.8 (inner right), and a = 1.0 (outer right). The linear power spectrum,  Centre: Relative differences between the matter power spectra of the cvG and QCDM models. A Savitzky-Golay filter has been used to smooth ∆P δδ (k)/P δδ,QCDM (k) for k > 0.2 h/Mpc. Each panel compares linear perturbation theory (black dotted), to results obtained from full simulation (coloured solid). The vertical grey shaded region in each panel indicates where k > k Ny where k Ny is the Nyquist frequency. Bottom: The signal-to-noise ratio of the difference between the cvG models and their QCDM counterpart. P (11) δδ (k; z), is obtained by multiplying the initial matter power spectrum at z ini = 49, P δδ (k; z ini ), with [D(z)/D(z ini )] 2 . The nonlinear matter power spectra are measured from particle snapshots using the POWMES 2 code [42]. The mean P δδ of the five realisations per cosmology is shown as a coloured line while the standard deviation is indicated as shaded region. The standard deviation is largest at large scales (k 0.1 h/Mpc) due to cosmic variance and the limited simulation size. The vertical shaded region near the right edge of each panel indicates the regime of k beyond the Nyquist frequency 3 .
The centre row of Fig. 2 shows the relative differences, Eq. (4.3), of the matter power spectra. The relative difference has been smoothed to remove noise at scales k > 0.2 h/Mpc, using a Savitzky-Golay filter of second order with a kernel width of 13 data-points [43]. The power spectrum results agree with the results found in Ref. [29] and extend them by including larger scales and measurement uncertainties.
The bottom panel of Fig. 2 shows the SNR of the difference between cvG cosmologies and their QCDM counterpart. From it we can conclude that the SNR is proportional to k while it is inversely related toβ 3 .
The real-space position of tracers of the matter distribution are not directly measurable, preventing us from comparing P δδ to observations, which rely on the redshift measurement to infer distances. The reason is that peculiar velocities (i.e., additional velocities to the Hubble flow) of the tracers distort the redshift signal along the line of sight. Thus, P δδ is different from its counterpart in redshift space, P s δδ , which becomes anisotropic despite the statistical istropy of the Universe; on large scales the two are related by the linear Kaiser formula where µ is the angle between the wavevector and the LOS, and f is the linear growth rate defined as The Kaiser formula can be improved down to quasi linear scales with additional information about the auto power spectrum of the velocity divergence 4 , θ = ∇ · v, denoted as P θθ , as well as their cross spectrum P δθ , since the velocity field is more sensitive to tidal gravitational fields compared to the density field on large scales [44][45][46].
The first row of Fig. 3 compares the linear velocity divergence power spectrum (black dotted lines) and measured nonlinear (coloured) simulations, at a = 0.6 (outer left), a = 0.7 (inner left), a = 0.8 (inner right), and a = 1.0 (outer right). The linear power spectrum P (11) θθ (k; z) can be related to P (11) δδ (k; z) through the zeroth-moment of Eq. (3.1), yielding the continuity equation, On linear scales we can assume that the quadratic terms in Eq. (4.7) vanish leaving use with Thus, the linear power spectrum of the velocity divergence is given by This relation is expected to fail on non-and quasi-linear scales, as velocities grow more slowly than the linear perturbation theory predicts. Therefore, any differences in P (11) θθ between the different cvG models will appear on these scales.
In order to measure the non-linear P θθ from the numerical simulations, we first use a Delaunay tessellation field estimator (DTFE 5 , [47]) to obtain the volume weighted velocity divergence field on a regular grid. This procedure constructs the Delaunay tessellation from the dark matter particle locations and interpolates the field values onto a regular grid, defined by the user, by randomly sampling the field values at a given number of sample points within the Delaunay cells and then taking the average of those values. For our 500h −1 Mpc simulation boxes, we generate a grid with 512 3 cells. From that we then measure P θθ using the public available code nbodykit 6 [48].
We can see from the top row of Fig. 3 that the results of the simulations for all models have approached the linear theory prediction on scales k 0.1 h/Mpc for all times. On these scales, the time evolution of the power spectrum of all models is scale independent and, the relative difference encapsulates the modifications to the time evolution of P (11) θθ via H and f in Eq. (4.9). On smaller scales, the formation of non-linear structures tends to slow down the coherent (curl-free) bulk flows that exist on larger scales. This leads to an overall suppression of the divergence of the velocity field compared to the field theory results for scales k 0.1 h/Mpc.
A careful look into the relative difference ∆P θθ (k)/P θθ,QCDM (k) in the bottom row of Fig. 3 also reveals a number of other interesting features on all scales. Firstly, we see that the wavenumber at which linear theory and simulation results for ∆P θθ (k)/P θθ,QCDM (k) agree, k * , depends both oñ β 3 and the scale factor. The value of k * is pushed to ever larger scales as a → 1 andβ 3 → 0. A similar observation has been made by [32] for the DGP model. Hence, this is important for the growth rate measurement from redshift distortions. Secondly, on small scales, k 1 h/Mpc, we can see how deviations from QCDM are suppressed by the screening mechanism, reflecting the fact that inside dark matter haloes the screening is very efficient. As also shown by ∆P δδ (k)/P δδ,QCDM (k), the screening mechanism becomes more effective asβ 3 → 0. Thirdly, for a → 1 andβ 3 → 0 we However, just as the transverse mode of the Proca field, ∇i × Bi, it has a much smaller magnitude than its divergence and is thus neglected in SPT. 5 The code is in the public domain, www.astro.rug.nl/ voronoi/DTFE/dtfe.html. 6 The code is in the public domain, nbodykit.readthedocs.io. see a growing peak that for the case ofβ 3 = 10 −6 protrudes above the linear theory prediction at k ∼ 0.7 h/Mpc. A similar feature was also observed by [32] for the DGP model. The difference of P θθ between the cosmological models compared to its magnitude is very small at early times, e.g., at percent level for all models when a 0.6, but increases rapidly over time, reaching 35% forβ 3 = 10 −6 at a = 1.0. This is unlike the behaviour of ∆P δδ (k)/P δδ,QCDM (k) which increases much more slowly and only reaches ∼ 5% forβ 3 = 10 −6 at a = 1.0. This difference is because the velocity field, being the first integration of the forces, responds more quickly to a rapid growth of the fifth-force magnitude than does the matter field, which is the second integration of the forces. It shows the rapid increase of the linear growth rate of the cvG model at late times (a 0.8), and suggests that redshift-space distortions (RSD) in this time window can be a strong discriminator of this model.

Matter bispectrum
As we have mentioned, even if cosmological fields are initially Gaussian, they inevitably develop non-Gaussian features as the dynamics of gravitational instability is nonlinear. Consequently, the structures found in the density field can no longer be fully described by two-point statistics alone, and higher-order correlation functions are needed in order to unlock additional information, in particular regarding the nature of gravitational interactions. To obtain first impressions of this information we use the Fourier space counterpart of the three-point correlation function, the bispectrum, which is receiving increased attention in the recent literature, not only for making more accurate predictions (see, e.g., [49][50][51][52]), but also as a probe of effects beyond ΛCDM (e.g. [53][54][55][56][57][58][59]).
We restrict ourselves to the study of the matter field in real space at z = 0, for which the bispectrum is given by with the three wave vectors forming a closed triangle. As the study of the effects on bispectrum due to modifications to GR are still in its infancy, we shall be as comprehensive as possible by considering all possible triangle configurations between the two extreme scales k min and k max , given a specific bin width ∆k 1 = |∆k 1 | for each side. A detection of strong configuration dependence can be regarded as a compelling motivation to further investigate higher-order statistics. It would allow us to disentangle the modified gravity signal from other potential cosmological effects, which might be degenerate in two-point statistics and other alternative measures. The top panel of Fig. 4 compares the bispectrum of equilateral triangles at the tree-level (dotted line), to the measurements (solid line). It furthermore contains the measured bispectrum of squeezed triangles (long dashed), folded triangles (short dashed), and all other triangle configurations (scattered dots). Vertical lines are spaced ∆k = |∆k| apart. As we assume a primordial Gaussian random field, we can apply the Wick theorem to write the bispectrum as products of power spectra summed over all possible pairings. Thus, the lowest-order bispectrum that is able to capture non-Gaussian features at late times has to expand one of the fields in the correlator of three Fourier modes to second order, yielding (4.11) where δ (n) is given in Eq. (4.2), the primed ensemble average indicates that we have dropped the factor of (2π) 3 as well as the momentum conserving Delta function, and "cyc." stands for the two remaining permutations over k 1 , k 2 and k 3 . Note here, that we have assumed that SPT gives an appropriate description of perturbations in the cvG model and does not fail to include further mode couplings that might be introduced through the additional Proca vector field. We will see below that this is indeed an excellent approximation. The resulting bispectrum scales as square of the linear power spectrum, P (11) , and exhibits a strong configuration dependence as it is directly proportional to the second-order perturbation theory kernel, F 2 , which is given by, To measure the bispectrum from the simulations, we first use fourth-order density interpolation on two interlaced cubic grids [60] of N = 256 cells per side. Next, we measure B(k 1 , k 2 , k 3 ) using an implementation of the bispectrum estimator presented in Ref. [61]. Starting from k min = 2k f = 0.025 h/Mpc, where k f denotes the fundamental mode, we loop through all configurations satisfying k 1 > k 2 > k 3 and k 1 ≤ k 2 + k 3 (the triangle closure condition). We stop after the values of k, which are evenly spaced by ∆k = 2k f , reach the k max = 1.0 h/Mpc, up until which point the shot noise is sub-dominant. With these settings -which are chosen to keep memory consumption at bay, as it would increase rapidly otherwise -we obtain a total of 5910 distinct triangle configurations.
The top panel of Fig. 4 shows that the tree-level prediction B (211) (dotted line) for the equilateral configuration converges to the simulation measurements of B (solid line) on k ≈ 0.07 h/Mpc, which is agreement with P δδ (k) and [62]. In this panel we have also indicated the folded, squeezed and equilateral configurations by lines (see the legends). It does not come as a surprise that the measured bispectrum for equilateral triangles is consistently lower than all other configurations as in our considered range of k, the power spectrum decreases with increasing k (as can be seen in Fig. 2). The folded triangles, on the other hand, tend to have the largest amplitude, while the squeezed triangles are in between.
The middle panel of Fig. 4 shows the relative difference, Eq. (4.3), of the bispectrum of equilateral triangles at the tree-level (dotted line), and measurements (solid line); for the latter the bispectra for all triangle configurations are indicated by scattered dots. Again, the results which correspond to equilateral, squeezed and folded triangle configurations are shown by lines (the same line styles as in the top panel). We can draw the following conclusions. Firstly, as it is the case for matter and velocity divergence power spectra, the tree-level bispectrum is a good estimator on large scales (k < k * ) while the exact value of k * depends on redshift and the model parameterβ 3 . However, we can see that in general linear theory gives accurate predictions of ∆B/B QCDM at k < k * ∼ 0.1 h/Mpc for all models. Compared to the matter power spectra, the relative difference of the bispectra is roughly twice as large as ∆P δδ /P δδ,QCDM , monotonically increasing from 1% forβ 3 = 100 to ∼ 9% for β 3 = 10 −6 . Secondly, the order of triangle configurations yielding the largest signal is reversed to the top row, with the equilateral triangles yielding the largest relative difference between cosmologies with fifth force and those without, while squeezed and folded triangles seem to converge to the same relative difference for larger values ofβ 3 . This is in agreement with [56], who arrived at a similar conclusion for f (R) and DGP cosmologies.
The bottom panel of Fig. 4 shows the SNR of the difference between cvG cosmologies and their QCDM counterpart. Three general trends are revealed: Firstly, an enhancement in the bispectrum signal with increasingβ 3 relative to QCDM, as we have seen in the middle panel above. Secondly, the SNR significantly increases towards smaller, nonlinear, scales. Thirdly, there is no clear trend which triangular configuration results in the highest SNR. The median taken over the range 0.1 < k [ h/Mpc] < 1 for each cvG cosmology is: 0.88 (β 3 = 10 −6 ), 0.77 (β 3 = 1), 0.54 (β 3 = 10) and 0.22 (β 3 = 100), respectively.
A very useful statistical quantity, that isolates the configuration dependence of the triangles by removing the propagator corrections from the modified Poisson equation (contained in the nonlinear power spectrum), is the reduced bispectrum, The relative difference between the reduced bispectra for the cvG models and their QCDM counterpart is displayed in the top row of Fig. 5. We indeed see how the strong scale dependencies of ∆B/B QCDM are removed, leaving only sub-percent deviations. The SNR of the difference of Q between the cvG models and their QCDM counterpart (not shown) revealed a very weak signal on all scales for all models, with a median of ∆Q/σ 0.05. Therefore we shall not try to interpret the trends revealed by the individual cvG models, and instead conclude that Q is very weakly dependent onβ 3 .
To quantify how much extra mode coupling the cvG models have experienced compared to their QCDM counterpart beyond the leading term, F 2 (defined in Eq. (4.12)), we can divide the reduced Figure 5: Top: relative difference between cvG models and their QCDM counterpart of the reduced bispectrum measurements, Q. Bottom: relative difference between cvG models and their QCDM counterpart of the ratio between the measured reduced bispectrum and its tree-level approximation, Q (0) . Each data point corresponds to one of 5910 triangle configurations (see the text for more details). The lines represent equilateral (solid), squeezed (long dashed), and folded (short dashed) triangle configurations as in Fig. 4. bispectrum by its tree level term to define a new quantity, . (4.14) The relative difference between the R of the cvG models and their QCDM counterpart is displayed in the bottom row of Fig. 5. Again, the results are in the sub-percent level and the SNR of the difference of R between the cvG models and their QCDM counterpart (not shown) reveals a very weak signal on all scales for all models, with a median of ∆Q/σ 0.06.
The fact that for Q and R the relative difference between the cvG models and QCDM is fairly small, suggests that the fifth force in the cvG model does not produce substantial extra mode coupling corrections. This is a useful result because it means that the cvG effect mainly enters through the modified growth factors, which simplifies the modelling of the bispectrum. We stress that this does not imply that the bispectrum is incapable of placing additional constraints on the cvG models. That is because the bispectrum has a different dependence on the growth factors than the power spectrum and its configuration dependence is useful in breaking degeneracies with other parameters, e.g. parameters that describe the background model or galaxy bias, such that the combination of the two statistics can still be expected to yield significant improvements.
Finally, let us note again that here we have only looked at the bispectrum of the matter density field, rather than the halo or galaxy fields. We have tried haloes, but due to the box size and resolution in our simulations, the results are noisy and the model differences unclear. Therefore we have decided not to show them here.

Halo statistics
This section is devoted to a detailed study of halo properties. Haloes are identified using two different algorithms, as they give complementary information about the haloes and can serve in some cases as verification. Firstly, we use the algorithm developed by [63] to find friends-of-friends groups to represent the 'main' haloes, and then run SUBFIND to identify substructures in the 'main' haloes (from now on we shall refer to the halo and subhaloes identified in this way as SUBFIND halos). Secondly, we use ROCKSTAR 7 [64] to identify FOF haloes in the 6D phase space where substructure is more easily identifiable (from now on we will refer to these as ROCKSTAR haloes). In most of this section we show results of SUBFIND haloes, although we have checked that the ROCKSTAR haloes give similar results. We use ROCKSTAR haloes to study the halo concentration mass relation, because this is directly measured by ROCKSTAR.
Note that, in principle, the unbinding procedure employed by the halo finding algorithms would need to be modified due to the presence of the fifth force induced by the Proca field. However, [65] found the effect of this modification to be quite small for chameleon models. Also, we will see below, the fifth force in the cvG models is strongly suppressed by Vainshtein screening, and so we expect its effect will be even smaller here. Thus, we use identical versions of SUBFIND and ROCKSTAR for the different cosmologies.
We compare the cvG models to their QCDM counterpart in the same way as we have done in Sec. 4 via Eq. (4.3) and Eq. (4.4).

Halo mass function
We start the analysis of the halo populations with the one-point distribution of halo masses -the halo mass function (HMF). The halo mass is defined as the mass enclosed in the spherical region of radius R 200 around the centre of the over-density, within which the mean density is 200 times the critical density ρ c at the halo redshift, In the top row of Fig. 6 we show the cumulative HMF, n(> M 200c ), which is the number density of dark matter haloes more massive than the given M 200c , at a = 0.6 (outer left), 0.7 (inner left), 0.8 (inner right) and 1.0 (outer right). The bottom-up picture of structure formation, i.e., small-scale objects collapse first and merge to form increasingly massive objects as time proceeds, is clearly visible, which follows from the fact that in our model dark matter is cold.
The bottom row of Fig. 6 shows the relative difference between the cvG models and their QCDM counterpart. The median of SNR of the differences between the models over the range shown in the figure is: 7.1 (β 3 = 10 −6 ), 6.4 (β 3 = 1), 5.5 (β 3 = 10), 2.9 (β 3 = 100). We find good agreement with [18], and have verified that the result is consistent between SUBFIND and ROCKSTAR. The fifth force enhances the abundance of dark matter haloes in the entire mass range probed by the simulations, with the enhancement stronger at late times and for high-mass haloes, which mimics the effect of the csG model [66]. This is to be expected because the strength of the fifth force increases over time [29]. Note that for massive haloes the increase in abundance is mainly due to an increase in individual halo masses, as can be seen from the top panels: we remark that more massive haloes are not necessarily more strongly screened in Vainshtein models (see, e.g., Fig. 8 of [67]), and the enhanced gravity around these massive haloes helps to bring more matter from their (matter-rich) surroundings to their vicinity, allowing them to grow larger. On the other hand, models with more efficient screening, such asβ 3 > 1, show a more restrained enhancement of the HMF.

Two-point correlation functions
The configuration-space counterpart of the matter power spectrum, P δδ , presented in Sec. 4.1, is the two-point correlation function (2PCF), ξ(r). In principle these two measures would carry the same information, but in practice this is not guaranteed since our analyses are restricted to a finite range of scales, and moreover, configuration and Fourier space statistics are impacted by different systematic effects.
For this analysis we use SUBFIND haloes, since these catalogues contain the subhaloes which can be proxies of satellite galaxies, and without which ξ(r) would decay at r 1-2h −1 Mpc due to the halo exclusion effect. We show their respective 2PCFs in the top row of Fig. 7 for a = 0.6 (outer left), a = 0.7 (inner left), a = 0.8 (inner right) and a = 1.0 (outer right). As expected, the 2PCFs drop off with halo separation, and can be well described by a power law across the entire range of scales probed here.
The relative difference between the 2PCFs of the cvG models and their QCDM counterpart for SUBFIND haloes is shown in the bottom row of Fig. 7. As for the power spectrum of the matter field, Fig. 2, we see more enhanced clustering for smaller values ofβ 3 . However, the cvG enhancement for halo clustering is smaller than for matter clustering, implying slightly smaller halo biases in stronger Note that to prevent the plot from appearing cluttered we have not shown the results for the cvG models. Furthermore, we included the standard deviation as a shaded region, but it is too small to see. Bottom: The relative differences between models. The cvG model for four values ofβ 3 = (10 −6 , 1, 10, 100) are shown, indicated by a blue, green, orange and red line respectively. The shaded regions are the standard deviations among the five simulation realizations. cvG models. This is because haloes are biased tracers of the dark matter field, and their bias generally decreases over time, as structure formation progresses: the enhanced gravity in cvG models simply speeds this up. Note that the enhancement of the halo 2PCF is nearly constant down to ∼ 3h −1 Mpc, consistent with the behaviour of the matter power spectrum (cf. Fig. 2), and reflecting the fact that in the cvG model the growth factor is enhanced in a scale-independent way in the linear regime.

Mean halo pairwise velocity
As outlined earlier, it is quintessential to develop a theoretical model of the pairwise velocity statistics as well as the real-space correlation function for cosmological analyses with redshift surveys, such as Euclid and DESI. Although we do not strive to actually test the cosmological models investigated here, we measure the relevant quantities to gain an intuition of how they are affected by the cvG model and to aid future work.
For this analysis we use SUBFIND haloes, as they contain the smallest haloes and subhaloes and thus can enable measurements to smaller scales, including the virial motions of subhaloes inside main haloes. We show the measured mean pairwise velocities for the different models in the top row of Fig. 8, comparing linear estimates (dotted lines) to the simulation results (solid lines) at a = 0.6 (outer left), a = 0.7 (inner left), a = 0.8 (inner right) and a = 1.0 (outer right). The linear mean pairwise velocity, v ij , is intimately related to the 2PCF of the matter field, ξ(r), through the pair conservation equation, Eq. (5.2), just as P θθ is to P δδ (see Sec. 4.1) through the continuity equation, We can replace the 2PCF in Eq. (5.2) with its Fourier space counterpart in first order, P (11) δδ , using the first-order Bessel function j 1 , and obtain the linear theory prediction of v ij expressed as where b is the linear bias of halos, f is the linear growth rate and j 1 is the spherical Bessel function of order 1 [69]. To get the bias values used in the linear theory prediction for Fig. 8, cf. Eq. (5.3), we compute the halo power spectrum, P hh , divide it by the matter power spectrum, b 2 ≈ P hh /P δδ . Due to the sparseness of haloes, the shot-noise becomes sub-Poisson on larger scales than it does for dark matter particles. Therefore we restrict the calculation of b to scales where the relation stays approximately constant, 0.025 < k h −1 Mpc < 0.1. We find that at each scale factor, the different cosmological models have the same fitted value of b (averaged over all 5 simulation realisations) up to the second decimal. Beyond the second decimal b indeed increases withβ 3 as expected from the relation of ξ(r) and P δδ (k).    The relative difference between v ij of the cvG models and their QCDM counterpart is shown in the bottom row of Fig. 7, which seems to have converged to a constant value for all cvG models at scales r > 10 h −1 Mpc. As an example, forβ 3 = 10 −6 the relative difference settles on ∼ 0.15 for large scales, which is approximately half of ∆P θθ (k)/P θθ,QCDM (k) shown in Fig. 3, partially due to the fact that P θθ ∝ f 2 . If ROCKSTAR-halos are considered the same qualitative trend is found on the larger scales.

Redshift space clustering
Motivated by the results of the real space clustering and mean pairwise velocity, we carry on to study the halo 2PCF in redshift space. In real observations, instead of their radial distances, we measure the redshifts of galaxies. The conversion from redshift space to real-space galaxy coordinates is not only determined by the Hubble expansion, but also affected by the peculiar velocities of galaxies. This induces anisotropies on what would be an isotropic galaxy correlation function, known as redshiftspace distortions (RSD). RSD is a useful probe of the peculiar velocity field, and consequently the growth rate of matter. In particular, the quadrupole of the redshift-space galaxy correlation function is sensitive to the galaxy (or halo) pairwise infall velocity, which we have seen above can be strongly enhanced by the fifth force in the cvG model. We use haloes (subhaloes) as proxies of galaxies in this study.
The mapping of the halo coordinates from real space to redshift space is given by, whereẑ is the unit vector in the line of sight direction which we have chosen to be along the zaxis of the simulation box, assuming that the galaxies are far away from the observer (plane-parallel approximation). Thus, the anisotropic correlation function is given by where s is the halo separation vector, s its magnitude, s the halo separation along the line of sight direction, and µ = cos(s /s) is the cosine of the angle between s and the LOS. We measure ξ s (s, µ), using SUBFIND-halos for the same reason stated in the previous section, over 40 bins of µ = [0, 1] and 40 bins of s = [0, 40] h −1 Mpc. In order to increase the SNR ratio, it is helpful to project ξ s (s, µ) onto a one-dimensional object which depends on s only. Therefore, we decompose the measured ξ s (s, µ) into multipole moments using its Legendre expansion, where is the order of the multipole and L (µ) is the Legendre polynomial at the -th order. Inverting Eq. (5.6) and integrating over µ, we find As the redshift space correlation function is symmetric in µ, only even values of give a non-zero contributions. Of these, we study the two lowest multipoles: the monopole ( = 0), and the quadrupole ( = 2). We omit higher order multipoles (l ≥ 4), as they do not have a big impact on the estimation of the correlation function and are noisier than the monopole and quadrupole [70].
In the top row of Fig. 9, we show the monopole, ξ s 0 , and quadrupole, ξ s 2 , moments of the QCDM model, at a = 0.6 (outer left), 0.7 (inner left), 0.8 (inner right) and 1.0 (outer right). We limit Figure 9: Top: the monopole, ξ s 0 , and quadrupole, ξ s 2 , moments of the 2PCF in redshift space. The results are obtained by averaging over the five simulations for each cosmology (solid lines) and shaded region show the standard deviation over these realization, which we show only for QCDM to maintain clearness. We have not shown the cvG results to prevent the plot from appearing cluttered. Central and bottom: the relative differences of ξ s 0 and ξ s 2 respectively. Each column shows the results for a different scale factor: outer left: a = 0.6, inner left: a = 0.7, inner right: a = 0.8, outer right: a = 1.0.
the study to scales < 40h −1 Mpc which is roughly 1/10 of the simulation box size. We know, however, that the peak position of the baryon acoustic oscillations (BAO) will be affected by the cvG model, asβ 3 → ∞ converges to QCDM andβ 3 → 0 converges to the cosmology of the csG, both being different from ΛCDM. The csG model is known to be unable to reproduce the BAO position [25,28,71] (see however [72]).
The central and bottom rows of Fig. 9 show the relative differences between the cvG models and their QCDM counterpart, for the monopole and quadrupole, respectively. The quadrupole moment encodes the anisotropies induced by redshift distortions, and as it has been the case for ξ and v ij , the relative difference of the cvG model to its QCDM counterpart increases with a decreasing value of β 3 especially on scales > 20h −1 Mpc. This implies that with decreasingβ 3 the contours of the twodimensional 2PCF in redshift space, ξ s (s , s ⊥ ), are more squashed, which is a direct consequence of the enhanced growth rate and stronger matter fluctuations as could already be anticipated from the results shown in Fig. 2. The values of ∆ξ 2 /ξ 2,QCDM converge on large scales for each cvG model to approximately the same values as for ∆v ij /v ij ,QCDM . The median SNR at a = 1 (outer right panel), taken over the range 20 < s/(h −1 Mpc) < 40, is approximately equal up to 7.2 for the monopole and 3.5 for the quadrupole for the strongest cvG modelβ 3 = 10 −6 . Although the relative difference is larger in the quadrupole, the SNR values are larger for the monopole, which is because the quadrupole is sensitive to the pairwise infall velocity v ij , which has a larger scatter than the real-space correlation function (see Figs. 8 and 7) that dominates the monopole signal. The RSD quadrupole can be a more promising probe to constrain the cvG model if the statistical uncertainties can be reduced by large amount of data.

Concentration-mass relation
For dark matter haloes, the strongest effect of Vainshtein screening is perhaps in the density profiles. This is because the interiors of haloes are expected to be strongly screened, see e.g., [67,73,74]. The Vainshtein screening radius can be even larger in the csG model and cvG models withβ 3 → 0, than in the DGP model at late times [29], so we expect the screening to be strong and the internal properties of haloes protected by it from the influence of the fifth force.
The density distribution inside dark matter halos is well described by the universal Navarro-Frenk-White (NFW; [75,76]) profile, where ρ s and R s are the characteristic density and scale radius respectively, which can vary from halo to halo. Thus the halo mass, M 200c , can be obtained by integrating the NFW density profile  where v max = GM (< R max )/R max is the maximum circular velocity inside a halo, which occurs at r = R max 2.163R s for an NFW density profile. Note that we do not attempt to do a full fitting of the NFW profile Eq.  have excluded all haloes with fewer than 1000 simulation particles from this figure which, combined with the small box size of our simulations, allows us to analyse the c 200 -M 200c relationship for halo masses that span only one order of magnitude. Nevertheless, we can clearly see that the relationship follows a power law [77][78][79]. Note that the statistics is poor at large mass and early times, due to a lack of haloes.
Without the screening mechanism we would expect haloes in a Proca universe to be more concentrated than their counterparts in a QCDM cosmology, since the strength of gravity increases quickly at late times [29], which causes a faster steepening of the gravitational potential inside haloes, attracting more matter to the central region and leading to a steeper density profile [80]. However, in the cvG model in reality, just as for the csG model [66], inside haloes the Vainshtein screening is strong enough that there is little effect of the fifth force, as can be seen from the bottom panels of Fig. 10.

Weak Lensing statistics
In the final section we focus on the study of weak-lensing statistics. We start by analysing the lensing convergence field (κ) which can be used together with the matter power spectrum and bispectrum to circumvent the dependence on tracer bias (e.g., [81]), and end with an analysis of the abundances and tangential shear profiles of voids identified from WL maps [82,83].

Weak lensing convergence and peak statistics
Weak lensing (WL) is governed by the lensing potential, Φ lens , which is given by with Φ and Ψ being the two Bardeen potentials in the metric Eq. (2.7). Φ and Ψ are related to each other through the anistropic stress. At late times, since we neglect matter species such as photons and neutrinos, in the cvG and QCDM models, the anisotropic stress is negligible so that we have Φ = Ψ. Therefore, in the cvG model not only massive particles can feel deviations from GR, but also can massless particles, as the dynamical and lensing potentials are equal and can both be modified substantially in the case ofβ 3 → 0. This is in contrast to some other models of gravity, such as f (R) gravity and the DGP model. The relation between κ and Φ lens and how those quantities are solved 'on-the-fly' during the simulation run time was summarised in Sec. 3. Here we would like to be more explicit how Eq. where G is the gravitational constant and δρ the density contrast. However, as the expansion history is altered in QCDM compared to ΛCDM their κ field will not be the same. For the cvG model, where the fifth force and screening mechanism are included, the lensing potential is where β sDGP is the coupling strength between matter and the brane-bending mode in the sDGP model, and β and α are given by Eq. (2.14) and Eq. (2.13) respectively. This modification of the lensing potential will modify Eq. (3.4) in the linear regime as lens,2D (χ, β(χ))dχ, (6.4) in addition to the modified expansion history. Here χ, which is the comoving distance, should not be confused with the longitudinal Proca mode,χ. This simple rescaling does not account for the effects of the screening mechanism and can only be accurately predicted through simulations as used in this work.
It is important to note here, that we solve the integral of Eq. (6.4) between z = [0.08, 1.0], as we found that artefacts appear for theβ 3 = 10 −6 cvG model. The reason behind this might be explained through the failure of numerical computation of theχ field in under-dense regions. This is a problem which has been reported multiple times [17,18,84,85] and discussed in terms of the cvG model in [29].
The resulting κ map is shown in Fig. 11 for QCDM (left), together with the residual between QCDM and the cvG model, ∆κ = κ(β 3 ) − κ QCDM , forβ 3 = [10 −6 , 10 2 ] (centre and right respectively). All maps have been smoothed with a Gaussian kernel of width θ = 2.5 arcmin which we will abbreviate as S G . It is clearly visible how underdense and overdense regions are more pronounced forβ 3 → 0 while forβ 3 → ∞ the model approaches the behaviour of the QCDM cosmology.
In the middle panel of Fig. 11 we can see a number of 'dipole' features, where a positive-residual 'hot spot' (∆κ > 0) is aligned with a 'cold spot' (∆κ < 0). This is produced by the transverse (i.e., perpendicular to the line of sight) motion of the halo which contributes most for a given line of sight: for this case the κ peak in the left panel would have moved slightly, causing this dipole feature in the residual map. Such dipoles are harder to find in the right panel, again because forβ 3 → 0 the model behaves very similarly to QCDM, so that haloes move little compared with the latter case.
Another feature worth mentioning in the middle panel of Fig. 11 is that we can see that near the massive structures the convergence field is enhanced by over 10%. This is partly due to the increased halo masses, but most likely the dominant effect here is the fact that the Proca field can also modify the lensing potential, as mentioned above. While we shall not investigate it here, let us note that this means that weak lensing by galaxy clusters can be a potential probe to constrain this model. However, as in the case of csG [86], we expect that the constraining power of cluster lensing may be limited by Vainshtein screening in the vicinity of clusters. We shall see shortly that this strong enhancement of convergence can be detected in the convergence power spectrum (or the shear correlation function) which can probe large-scale variations of the lensing potential.
In observations, the WL signal is obtained by averaging the shearing of source galaxy shapes over a large number of source galaxies whose intrinsic ellipticity dominates over the physical tangential shear signal. This effect is known as galaxy shape noise (GSN) and is a main source of uncertainty on small angular scales. We include the GSN by modelling it as a Gaussian random field which we will denote as N GSN . Therefore we assume that N GSN is independent of the underlying κ. Furthermore, we assume that the correlation function of N GSN is a δ function, thus pixel values show no correlation. The standard deviation of the Gaussian distribution is given by where σ int is the intrinsic ellipticity dispersion of the source galaxies, θ pix is the width of each pixel,  and n gal is the measured source galaxy number density. We use σ int = 0.4 and n gal = 40 arcmin −2 , which match LSST specifications [87].
In the top row of Fig. 12 we show the results for the power spectrum of the κ maps (faint) and the κ-N GSN -S G maps (bright). We do not include the linear theory prediction, as it holds up to 10 2 and is thus outside of the range of multipoles we are able to extract from the maps. The left panel shows the absolute power spectra measurements for which we have not included the results for > 10 4 as such small angular scales are not well-resolved given our simulation resolution. In terms of the relative difference between the cvG models to their QCDM counterpart in the right panel, the curves show the expected behaviour that, on large angular scales ( < 10 4 ), the amplitude is higher Figure 12: Weak lensing statistics: lensing convergence angular power spectra (top), probability distribution function of the weak lensing convergence field (middle), weak lensing peak abundance plotted as a function of peak height (bottom). The results shown here are obtained using a 10×10 deg 2 partial sky-map for a redshift range z = [0.08, 1.0]. We show results of the κ maps (faint) and the κ maps including the galaxy shape noise map, N GSN , and smoothed with a Gaussian kernel of width θ = 2.5 arcmin (bright) for the cvG model variants (colour) and their QCDM (black) counterpart.
in the cvG models with smallerβ 3 . However, since we use a partial-sky map of 10 × 10 deg 2 , the power spectra in the left panel could suffer from a large sample variance. This, however, should not strongly affect the result of the relative difference, as it roughly cancels out. As we go to smaller angular scales, l → 10 4 , all cvG models converge toward their QCDM counterpart, which reflects the operation of the screening mechanism on small scales, e.g., inside haloes. Note that the smoothed maps behave similarly, though not identically, to the unsmoothed ones at 10 3 , while on smaller angular scales the smoothing significantly changes the model difference. This indicates a potential limitation on using the convergence power spectrum or shear two-point correlation function to test the cvG model, but we note that the large angular scales are where the model difference is most prominent anyway.
The middle row of Fig. 12 shows the one-point distribution of the κ maps (faint) and the κ-N GSN -S G maps (bright). It contains information on non-Gaussian aspects of the convergence field that are not included in the convergence power spectra. We can see that cvG models with smallerβ 3 have larger numbers of pixels with both high and low κ values. This behaviour is as expected because the fifth force in the cvG models helps to move more matter towards (from) dense (underdense) regions, as can be seen in Fig. 11. It is good to see that increasing theβ 3 parameter indeed leads a smooth transition to QCDM, which is what is needed to cure the problem of having too strong a lensing effect in the csG model. The same happens to the void γ t profiles too, as will be shown in the next subsection.
The bottom row of Fig. 12 shows the WL peak abundance for the κ-N GSN -S G maps. This result is useful on its own because WL peak statistics can be a useful cosmological probe (e.g., [88][89][90][91][92][93][94][95]) but will also be useful for the study of void identified through WL peaks in the next subsection. We identify peaks as pixels whose κ values are larger than those of their eight neighbours. For consistent definitions between the different cosmological models, we define the amplitude of ν of a map pixel as ν = κ σ GSN , (6.6) where σ GSN = 0.007 is the standard deviation of the N GSN -S G map generated using the LSST specifications given above. From the bottom panels of Fig. 12, we can see that forβ 3 → 0, there is a significant increase in the numbers of the high-amplitude peaks, which indicates that the fifth force strongly enhances the lensing signal of these pixels (note that the fifth force also increases the halo masses as found in Fig. 6, which also contributes to this). On the other hand, the abundance of small peaks (ν < 1) is reduced asβ 3 → 0, because some of the haloes that produce peaks with ν < 1 in QCDM have been able to produce peaks with ν > 1 in the cvG models. This trend agrees qualitatively with results found for the nDGP cosmology [96].

Cosmic voids
Cosmic voids are regions in the Universe where the densities of dark matter or tracers are low. In recent years it has been shown that voids (e.g., [97][98][99]) can be a useful probe for a variety of models (e.g., [96,[100][101][102][103][104][105][106][107][108][109][110][111][112][113]), including the test of modified gravity models that are featured by Vainshtein screening [102,104,107,110]. There are a large number of methods to find voids, and it has been argued that void identification based on WL convergence maps can lead to the better constraints of modified gravity theories [96]. This has motivated us to use voids from the two dimensional convergence field through the tunnel and watershed algorithms as the resulting void catalogues have been shown to be amongst the most promising [83]. Whilst the convergence profiles of voids allow for a simpler physical interpretation of the mass content, where positive and negative κ correspond to projected over-dense and under-dense regions, it is the tangential shear which can be measured directly in observations. Therefore, to offer a more straightforward comparison with observations, we study the void tangential shear profile γ t (r), which is related to the convergence profile through whereκ is the mean enclosed convergence within radius r.

Tunnels
The tunnel algorithm of [82,96,108] identifies voids based on a WL peaks catalogue. We will from now on refer to these voids as tunnels. We find peaks using the κ map smoothed by a compensated Gaussian kernel wither an inner kernel width of θ inner = 2.5 arcmin and a outer kernel width of θ outer = 15 arcmin, which we will abbreviate as S cG . The use of S cG instead of S G is motivated by the larger number of identified peaks, which again will results in more identified tunnels and thus better statistics. Each identified peak is placed into three categories based on Eq. (6.6): ν > [1,2,3].
For each category, a Delaunay tessellation with the peaks at the vertices is constructed. This produces a tessellation of Delaunay triangles, with a peak at the corner of each triangle, and no peaks within the triangles. Each Delaunay triangle is then used to construct its corresponding circumcircle, with the three vertices of the triangle falling on the circumcircle's circumference. This unique tessellation, by definition, produces circles which do not enclose any peaks. In order to increase the number of Figure 13: Top: the tunnel abundance as a function of their radii for the three WL peak categories: left: ν > 1, centre: ν > 2, right: ν > 3. Bottom: relative difference between the cvG cosmologies and their QCDM counterpart. tunnels, which is necessary because of the small area of our convergence maps, we use all possible tunnels, including neighbouring ones which have a large degree of overlap in our study. The top row of Fig. 13 shows the tunnel size distribution identified from peak catalogues of different significance: ν > 1 (left), ν > 2 (centre), and ν > 3 (right). The smallest tunnels are generated by the ν > 1 peak catalogue, which also produces the most tunnels, because the large number of peaks in this catalogue tends to partition the map into smaller Delaunay triangles. As the ν threshold increases, the typical tunnel size increases, however there are also fewer tunnels overall. This implies that each of the three categories should respond differently to the large scales modes of the κ map, and thus creating the tightest constraints through combined analyses. Due to our small sample size, this remains to be tested.
The bottom row of Fig. 13 shows the relative difference between the cvG models and their QCDM counterpart. It is interesting to observe, that while smaller tunnels (R v 0.2 deg) are more abundant in cvG withβ 3 → 0 than in QCDM it is vice versa for larger voids (R v 0.2 deg). This is a consequence of a higher abundance in WL peaks for the cvG cosmologies compared to their QCDM counterpart for all of our peak categories, see Fig. 12, more small voids and fewer large voids are found in cvG than in QCDM. Fig. 14 shows the tangential shear profiles, Eq. (6.7), of the three tunnel catalogues shown in Fig. 13. The profile are based on the κ-N GSN maps, as smoothing would dampen the void profiles and Figure 14: Top: Tunnel tangential shear profiles as a function of the scaled distance from the centre, r/R v , for QCDM (black) and cvG models withβ 3 = 10 −6 (blue), 10 0 (green), 10 1 (orange) and 10 2 (red). The shaded region indicates the standard deviation of all tunnels in the QCDM map (for clarity we do not show this for the other models). The shaded region indicates the standard deviation. Bottom: The relative difference between the cvG models and their QCDM counterpart. From left to right the panels are respectively for tunnels identified from peak catalogues with peak height ν > 1, 2 and 3. We do not show the standard deviation as they very large due to our small sample size. the differences between the cosmological models. We compute the γ t profiles statistics by stacking all voids in a given catalogue, weighting them depending on their size (the smaller the void, the less its statistical weight). To obtain the 1-σ error, indicated by the shaded region in the top row, we loop through 100 bootstrap resamples. We recover the typical tangential shear profile, which indicates that voids act as concave lenses. The extrema of the profile is located at r ≈ R v for all void categories and is increasing as the void sizes increase.
In the bottom row of Fig. 14 we can clearly see that the potential well get deeper asβ 3 → 0, reflecting the effects of enhanced structure formation and modified photon geodesics. We do not show the bootstrapped 1-σ error for the relative differences, as our sample size is too small.

Watershed
The watershed algorithm of [114] identifies voids based on the basins in the topographic map which is constructed from the κ map. To find the watershed basins, each pixel of the κ map is connected to its neighbour with the lowest κ value -a process that is repeated for successive neighbours until a local minimum emerges. All pixels connected to the same minimum in this way form one watershed basin, with ridges of local high κ values along the basin boundary. We could have used the WL peak catalogues to identify watershed voids, as is done for tunnels, but the results are generally very noisy [83]. To mitigate the impact of GSN, [83] found that the basin boundary should have a minimum κ Figure 15: Statistics for the watershed voids. Left: the cumulative void abundance as a function of the effective radius of the watershed voids, R v . Right: the tangential shear signal of these voids, as a function of the scaled radius from void centre, r/R v . The upper panels show the results for QCDM (black) and cvG models withβ 3 = 10 −6 (blue), 10 0 (green), 10 1 (orange) and 10 2 (red), while the lower panels show the relative (for the void abundance) and absolute (for the tangential shear profile) differences between the cvG models and their QCDM counterpart. The shaded region in the top row indicates the standard deviation of all profiles in the QCDM model. value of σ GSN /2, as it allows watershed basins that have been artificially split by spurious structures introduced by GSN to be re-merged. Unlike tunnels, the watershed voids are formed by a collection of Delaunay cells, and therefore have irregular shapes. We define the void centre to be the barycentre of all selected cells for a given watershed void, and the void radius R v as the radius of a sphere whose volume is equal to that of the void. The watershed algorithm has the advantage of simplicity from fewer free parameters in the void identification process, since no tracers are used, multiple WL peak catalogues do not need to be defined. However, Ref. [83] also find that tangential shear profile from the watershed algorithm is more susceptible to GSN than the tunnel algorithm.
The left column of Fig. 15 shows the watershed void abundance as a function of the void radius, R v , and the relative difference between the cvG models and QCDM. In contrast to tunnels, there are overall fewer watershed voids, and they never reach the large void size as tunnels do. This is because watershed voids by definition cannot overlap. Among the different models, little difference is found, apart from the large-R v end, where the cvG models produce up to ∼ 20% fewer voids than QCDM. The main reason for this is a change of void sizes, rather than a decrease in their number. This is likely due to the enhanced κ field magnitude in local overdensities residing in larger underdense regions, which means that these structures would more easily have κ > σ GSN /2 and therefore become basin boundaries in the cvG models, leading to a split of a large waterbasin into smaller ones.
The right column of Fig. 15 shows the tangential shear profiles, γ t (r), of watershed voids and their relative difference between the cvG models and their QCDM counterpart. They are smoother, wider, and shallower compared to all tunnel categories. However, both tunnels and watershed voids reach their tangential shear profile minimum at 0.9 − 1.1 R v . The error bars on the QCDM tangential shear profiles from the two algorithms are also similar in size, which suggests that both algorithms may offer similar constraining power, consistent with Ref. [83] which finds roughly similar tangential shear signal-to-noise ratios between the two algorithms. The relative differences between the cvG models and their QCDM counterpart peak at the minimum of the profile, with a 10% difference for cvG withβ 3 = 10 −6 , roughly the same as the relative difference found for tunnels in the same size range (which is the tunnel category for ν > 1).

Discussion and conclusions
In this paper, we have performed a thorough phenomenological study of a simplified version of the generalized Proca theory, the vector Galileon model (cvG). To study the impact of the cvG models free parameter,β 3 , we have run a set of five realizations of simulations forβ 3 = [10 −6 , 10 0 , 10 1 , 10 2 ] and their QCDM counterpart, resulting in a total of 25 simulations. The study relied on an adapted version of the ECOSMOG N-body code augmented with the ray-tracing modules of the Ray-Ramses algorithm. We used the five independent realisations for each model to create a light cone that covers a field of view of 10×10 deg 2 from z = 0.08 to a source redshift of z = 1 (cf. Sec. 3 and Fig. 1). This allows us to study the matter, halo and weak lensing statistics. In the following we shall summarise the results of each those three topics.
The study of dark matter field statistics finds good agreement with [29] about the matter power spectrum (P δδ , cf. Sec. 4.1 and Fig. 2), but extends the results of that paper by including larger scales and showing statistical uncertainties. In addition: • the simulation measurements of the velocity divergence power spectrum (P θθ , cf. Sec. 4.1 and Fig. 3) converge to the linear-theory prediction on scales k 0.1 h/Mpc for all times, while for k 0.1 h/Mpc we reproduce the well-known result that P θθ is suppressed compared to the linear theory results. The relative difference, ∆P θθ (k)/P θθ,QCDM (k), shows that the wavenumber at which linear theory and simulation results agree reasonably, k * , is pushed to ever larger scales as a → 1 andβ 3 → 0. Finally, for a → 1 andβ 3 → 0 we see a growing peak that for the case ofβ 3 = 10 −6 protrudes above the linear theory prediction at k ∼ 0.7 h/Mpc. A similar feature was also observed by [32] for the DGP model.
• for the matter bispectrum (B, cf. Sec. 4.2 and Figs. 4, 5), we find that the magnitudes depend on the triangle configurations, and increase in the order of equilateral, squeezed, and folded triangle configurations. However, this order is reversed when considering the relative difference. The relative difference confirms that, as it is the case for P δδ and P θθ , the tree-level bispectrum is a good estimator on large scales k < k * ∼ 0.1 h/Mpc, while the exact value of k * decreases with a → 1 andβ 3 → 0. We show that the enhancement of the bispectrum due to the fifth force is marginally stronger than in the case of power spectrum, but the reduced bispectrum shows that B/B QCDM is to a very good approximation equal to (P/P QCDM ) 2 . The scales at which we are able to measure the bispectrum do not show a strong signature of the Vainshtein screening.
The study of halo statistics is mostly based on SUBFIND cagalogues, as they contain the smallest haloes and subhaloes and thus can enable measurements to smaller scales, although where possible we have also cross-validated the results with FoF haloes. The main observations are the following: • the halo mass function (n(> M ), cf. Sec. 5.1 and Fig. 6) shows that the fifth force enhances the abundance of dark matter haloes in the entire mass range probed by the simulations, with the enhancement stronger at late times and for high-mass haloes. Models with a weaker fifth force, e.g., withβ 3 → ∞, show a more restrained enhancement of the HMF.
• the two-point correlation function (ξ(r), cf. Sec. 5.2 and Fig. 7) shows more strongly enhanced clustering for smaller values ofβ 3 , for which the fifth force is stronger. The enhancement of the halo ξ(r) is nearly constant down to ∼ 3 h −1 Mpc, consistent with P δδ , and reflecting the fact that in the cvG model the growth factor is enhanced in a scale-independent way in the linear regime. However, the enhancement in halo clustering is weaker than in matter clustering, for all models at all times.
• the relative difference of the mean halo pairwise velocity (v ij , cf. Sec. 5.3 and Fig. 8) remains constant for all cvG models at scales r > 10 h −1 Mpc, in very good agreement with linear-theory prediction. For the latter, we have measured the halos bias, b, for four different scale factors through the relation between the halo and matter correlation functions. The resulting measurements of b for the different models are similar, but show a slight decrease asβ 3 → 0, as the fifth force enhances matter clustering more than halo clustering, as mentioned above.
• the redshift space halo clustering (ξ (s), cf. Sec. 5.4 and Fig. 9) is sensitive to the halo pairwise velocity and hence the fifth force. The relative difference between cvG and QCDM can be up to ∼ 3 times larger for the quadrupole, ξ 2 (s), than for the monopole, ξ 0 (s), although its SNR is ∼ 0.5 times smaller on the range 20 < s h/Mpc < 40 due to larger statistical uncertainly in the halo velocity field. Future data of redshift space distortions should provide strong constraints onβ 3 .
• the result of the halo concentration-mass relation (c 200 , cf. Sec. 5.5 and Fig. 10) shows that in the cvG model, just as for the csG model, the Vainshtein screening is strong enough inside haloes that there is little effect of the fifth force.
Our final section concerns the properties of the weak lensing convergence, peak and void statistics, where voids are identified using the tunnel and watershed algorithms. The main results are the following: • the difference of the convergence map (κ, cf. Sec. 6.1 and Fig. 11) between QCDM and cvG forβ 3 = 10 −6 shows that around massive structures the convergence field is enhanced by over 10%. However, we caution about taking this as an indication that weak lensing by galaxy clusters can be a potential probe to constrain this model, as we have not performed an analysis of stacked weak lensing convergence profiles.
• the relative difference of the angular power spectrum (C , cf. Sec. 6.1 and Fig. 11) is largest on linear scales 3×10 2 , reaching ∼ 30% forβ 3 → 0. These scales are also where the smoothing of the map has little impact on the relative difference. For higher multipoles the model differences reduce.
• the relative difference of the probability distribution function of κ (PDF(κ), cf. Sec. 6.1 and Fig. 11) shows that cvG models withβ 3 → 0 have more pronounced under-and overdense regions.
• the relative difference of the weak lensing peak abundance (N p , cf. Sec. 6.1 and Fig. 11) shows larger (smaller) numbers of high-(low-)amplitude peaks for ν > 1 (ν < 1) in the cvG models withβ 3 → 0, because the fifth force enhances the convergence values of the peak pixels.
• the relative difference of the tunnel and watershed void abundances (N (> R v ), cf. Sec. 6.2 and Fig. 13, 15) shows fewer large-sized voids in the cvG cosmologies compared to their QCDM counterpart, since they produce more weak lensing peaks which splits large voids into smaller ones (for the tunnel case), or increase the convergence values so that the regions satisfying the chosen void definition criterion shrink in size (for the watershed case).
• the relative difference of the tangential shear profile for tunnels and watershed voids (cf. Sec. 6.2 and Fig. 14, 15) peak at approximately the void radius, with up to 10% difference for the cvG model withβ 3 = 10 −6 (similar to what has been observed in the convergence maps), and the model difference decreases asβ 3 → ∞.
Overall, we find that for the cvG model studied here, the fifth force effect is strongest on velocity and lensing statistics. The former is because velocity is the first integration of acceleration, and thus reacts quickly to the enhancement of gravity due to the fifth force, which happens only at late times; the matter density field, in contrast, reacts more slowly as the second integration of acceleration. The latter is because in the cvG model, unlike for some other MG models, photon geodesics are affected in two different ways: (1) indirectly, by the modified growth of matter fluctuations, and (2) directly, by the fifth force. This suggests that redshift space distortions and weak lensing shear correlation functions can both be promising cosmological probes to constrain theβ 3 parameter in this model. On small scales, the models are generally more difficult to constrain because the screening mechanism suppresses the fifth force effect; for example, internal properties of haloes, such as the concentrationmass relation, are insensitive to the fifth force. Another potentially useful way to constrain this model is by cross-correlating galaxies with the integrated Sachs-Wolfe effect [28], because asβ 3 → 0 the fifth force becomes stronger, causing the lensing potential to getting deeper rather than shallower [29] as suggested by observations. This possibility will be investigated in future. What is a bit surprising is that weak lensing by voids do not seem to be as promising a probe, even though the lensing potential is significantly modified in low-density regions: perhaps this is because weak lensing is a cumulative effect along the line of sight, and this strong effect in low-density regions is somehow cancelled out by the weaker effects in high-density regions.
Recently, various studies to constrain the generalized GP theory using cosmological observations have been conducted, see, e.g., [25,26,71]. These studies focused on general nonlinear functional forms for G 2,3 , because linear forms of these functions, such as the models studied here, have been found as a poor fit to observational data. However, as suggested by [29], adding massive neutrinos with significantly nonzero mass (see, e.g., [72]) may be a way to make the GP model with linear G 2,3 agree better with data. This possibility will be studied in a follow-up work, and correspondingly we hope to include massive neutrinos in future simulations.