{'fraction_non_alphanumeric': 0.03784279403184129, 'fraction_numerical': 0.013920146703342502, 'mean_word_length': 4.3250480840361, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Galaxy number counts to second order and their bispectrum

We determine the number counts to second order in cosmological perturbation theory in the Poisson gauge and allowing for anisotropic stress. The calculation is performed using an innovative approach based on the recently proposed"geodesic light-cone"gauge. This allows us to determine the number counts in a purely geometric way, without using Einstein's equation. The result is valid for general dark energy models and (most) modified gravity models. We then evaluate numerically the relevant contributions to the number counts bispectrum. In particular we consider the terms involving the density, redshift space distortion and lensing.


I. INTRODUCTION
Cosmology has entered the precision era. So far, mainly via observations of the cosmic microwave background (CMB) most recently performed by the Planck satellite [1]. But also high precision observations of cosmic large scale structure (LSS), i.e. the distribution of matter in the Universe, are being performed, see e.g. [2,3], and are under construction or planning, see [4,5] 1 and references therein.
The interpretation of observations of LSS is more complicated than the CMB. First of all we observe galaxies and not the density field, this is the biasing problem, and secondly, fluctuations on small scales become large so that non-linearities become important. Nevertheless, if we can handle these problems, a detailed observation of LSS is very interesting because it contains much more information than the CMB. LSS is a three dimensional dataset, and if we can observe it from the scale of about 1h −1 Mpc out to the Hubble scale, we have in principle ≃ (3000) 3 = 2.7 × 10 10 independent modes at our disposition. This has to be compared with the CMB with about 10 7 modes (including E-polarization).
Furthermore, non-linearities in the CMB are very small and CMB fluctuations are very Gaussian. Hence the power spectrum captures nearly all the information and (reduced) higher order correlations are very small. For large scale structure, at least on intermediate to small scales, non-linearities become important. The analysis of the highly non-linear regime can probably only be done via full numerical simulations, but in the intermediate, weakly non-linear regime higher order perturbation theory can be applied. As soon as non-linearities are present, higher order correlations develop i.e. non-trivial three-point and reduced four-point functions are generated. These functions contain additional interesting information on gravitational clustering which can be very important for the Dark Energy problem or the contribution from massive neutrinos.
The present work inscribes in this direction. We first compute a basic LSS observable, namely the galaxy number counts to second order in perturbation theory and then determine the induced bispectrum. Our approach is novel in several aspects.
We describe the number counts as a function of direction and observed redshift. This has the big advantage that it is directly related to observables and that it is model independent contrary to the galaxy correlation function in real space or its power spectrum. The latter need the knowledge of a background Universe to convert measured redshifts and angles into distances and are therefore not well adapted, e.g., to estimate cosmological parameters. This method was introduced in [6] and implemented and tested numerically in [7,8]. Galaxy number counts have the advantage of being directly measured and their correlation function ξ(θ, z, z ′ ) or the power spectrum in the form of C (z, z ′ ) can be directly compared to observations.
The first order fluctuations of the number counts have been calculated in [6,[9][10][11]. In this work we go to second order. Several second order calculations have already been performed. Expressions for the redshift perturbations and cosmological distances were obtained in [12][13][14] and applications are discussed in [15][16][17]. Similar expressions are also obtained in [18,19]. Results for the Newtonian density fluctuation were obtained in [20] and second order perturbation theory of lensing is discussed in [21]. Also expressions for the second order number counts have just appeared in [22][23][24]. However, since these calculations are very involved, it is useful to have results from independent groups, obtained in different ways, which can then be compared. Here we employ a novel method to calculate the second order number counts. We use the geodesic light-cone gauge introduced in [25]. Actually, if we would know the density fluctuation exactly, our expression for the number counts would be fully non-perturbative. Then we perform a second order coordinate transformation to translate the result to the familiar Bardeen potentials, i.e. the metric perturbations in longitudinal or Poisson gauge. We also allow for anisotropic stresses. We however neglect second order vector and tensor contributions which can be easily obtained from the first order results by inserting the second order vector and tensor perturbations induced by the Bardeen potentials, which have been calculated in the past [26][27][28]. We write the final result fully in terms of the Bardeen potentials, the peculiar velocities and the density fluctuations. Of course, the final expressions remain complicated, but they can be readily used and compared to other results in the literature.
Finally, we use our expressions to compute some of the most important new terms in the bispectrum, namely the contribution from redshift space distortion and the one from lensing, and we evaluate them numerically and compare them with the usual term from weakly non-linear Newtonian gravitational clustering.
The remainder of this paper organised as follows: in Section II we lay out the basic quantities which we want to compute. In Section III we describe the geodesic light-cone gauge and its transformation to the Poisson or longitudinal gauge. In Section IV we compute the second order number counts first in geodesic light-cone gauge and then we transform it to Poisson gauge. We also present a simplified formula containing only the dominant terms. In Section V we derive and discuss some of the most relevant terms in the bispectrum for which we also plot numerical results using classgal [7]. In Section VI we conclude. Some useful relations are collected in Appendix A.

II. GALAXY NUMBER COUNTS
What we really observe in galaxy surveys is the number N of galaxies in a redshift bin dz and solid angle dΩ as a function of the observed redshift z and of the observation direction defined by the unit vector −n (in this work n denotes the direction of propagation of the photon). Directly from these observable quantities we define the fluctuation of galaxy number counts as where ⟨⋅⟩ denotes the average over directions. If we neglect bias, the observed number of galaxies can be expressed in terms of the density ρ (n, z) and the volume V (n, z) at fixed redshift according to N (n, z) = ρ (n, z) V (n, z) .

(2.2)
To determine the galaxy number counts ∆ (n, z) up to second order in perturbation theory we need to expand the redshift density perturbation and the volume to second order. The volume and the density have to be defined in function of the observed redshift and direction of observation, then their expansion to second order can be given as Here, z denotes the observed redshift, which is also perturbed, and we define the time coordinate of a fiducial background model of a source by where η o denotes the conformal time of the observer. We also note that the background conformal distance is r s . Hence to some extent, the perturbations of the volume and the density above come from the perturbations of the redshift. In the next section we shall present the second order perturbations of the redshift, obtained in [14], and Section III is devoted to the calculation of the volume and the density fluctuations.
Once these perturbations are calculated, inserted in the definition (2.1), they yield to second order in perturbation theory the following result Here we use the fact that the angular mean of first order perturbations vanishes. We then identify the first order and the second order contributions to galaxy number counts, respectively, The first order contribution (2.7) has been already computed in [6,[9][10][11]. Here we are mainly interested in computing the second order term (2.8). For this purpose we adopt the so-called geodesic light-cone (GLC) coordinates defined in [25]. In terms of GLC coordinates we determine an exact solution which we then expand in a more conventional gauge, namely the Poisson gauge, to second order.

III. FROM THE GEODESIC LIGHT-CONE TO THE POISSON GAUGE
Let us first introduce the geodesic light-cone (GLC) coordinates [25]. They consist of a timelike coordinate τ (which can be identified with the proper time of synchronous gauge and, therefore, describes a geodesic observer that is static in this gauge [29]), of a null coordinate w and of two angular coordinatesθ a (a = 1, 2). We denote them with a tilde to distinguish these exact screen space angles from the background angles in the Friedmann metric which we introduce below.
The line-element of the GLC metric takes the form It depends on six arbitrary functions (Υ, the 2-component vector U a and the symmetric 2 × 2 matrix γ ab ). In matrix form: where γ ab and its inverse γ ab lower and raise the two-dimensional indices and U 2 = γ ab U a U b . The 2-dimensional space described by the coordinates (θ 1 ,θ 2 ) is the exact screen space with metric γ ab .
In an unperturbed Friedmann background, our coordinates correspond to cosmic time τ = t, background past light-cone w = η + r, where η is conformal time so that adη = dt, and the usual polar angles θ =θ 1 and φ =θ 2 . The background metric components are Υ = a, U a = 0 and γ ab dθ a dθ b = a 2 r 2 (dθ 2 + sin 2 θdφ 2 ). Here a denotes the cosmic scale factor.
Like in synchronous gauge, also the GLC gauge has some residual gauge freedom (for details see [13]). The condition w = constant defines a null hypersurface (∂ µ w∂ µ w = 0), corresponding to the past light-cone of a given observer, hereafter chosen to be the geodesic observer which moves with 4-velocity u µ = −∂ µ τ . The geodesic equation for u µ is (∂ ν τ ) ∇ ν (∂ µ τ ) = 0. Let us also recall that, in GLC gauge, the null geodesics connecting sources and observer are characterised simply by the tangent vector k µ = −ωg µν ∂ ν w = −ωg µw = ωΥ −1 δ µ τ (where ω is an arbitrary normalization constant), meaning that photons travel at constant values of w andθ a . This renders the calculation of the redshift particularly simple in this gauge.
We now give exact, non-perturbative, expressions for the redshift and the cosmological distances, like the luminosity and angular distances, in GLC gauge. As in the previous section, subscripts "o" and "s", respectively, denote quantities evaluated at the observer and the source space-time positions. We consider a light ray emitted by a geodesic source with 4-velocity u µ s = −(g µν ∂ ν τ ) s lying on the past light-cone of the geodesic observer defined by w = w o and on the spatial hypersurface τ = τ s . The light ray will be received by the geodesic observer at τ = τ o > τ s . The exact non-perturbative expression of the redshift z s associated with this light ray is then simply given by [25] ( On the other hand, in [13] an exact expression for the so-called Jacobi Map [30] is derived in GLC gauge and the following non-perturbative solution for the luminosity distance d L and the angular diameter distance d A is obtained: where γ denotes the determinant of the 2-dimensional matrix γ ab . We finally want to express our results in Poisson gauge. Neglecting vector and tensor contributions, the Poisson gauge (PG) metric [31] (sometimes also denoted 'Newtonian gauge' or 'longitudinal gauge') takes the form where the (generalized) Bardeen potentials Φ and Ψ are defined, up to second order, as follows: and we make no assumption on the anisotropic stress, so that Ψ and Φ may differ also at first order.
To express the redshift and the luminosity distance given in Eqs. (3.3) and (3.4) in terms of standard PG variables we have to transform the GLC gauge quantities to quantities in PG. Following [14], and employing suitable boundary conditions: i) the transformation is non singular around the observer position at r = 0, and ii) the two-dimensional spatial sections r = const are locally parametrized at the observer position by standard spherical coordinates (θ, φ) ≡ (θ 1 , θ 2 ), the coordinate transformation to second order is given by [14] where (γ ab 0 ) = diag(r −2 , r −2 sin −2 θ), and η in represents an early enough time when the perturbations (or better their integrands) were negligible. We have also introduced the zeroth-order light-cone variables η ± = η±r, with corresponding partial derivatives: Furthermore, we define Note that the integrals for τ go along the conformal time axis η at fixed r, while the integrals for w andθ a are along the background light cone, η − at fixed η + . With this we can obtain the non-trivial entries of the GLC metric of Eq. (3.2) to second order in terms of the variables (η, r, θ a ): (3.14) Here X (a⋯ Y ⋯b)⋯ denotes symmetrization of a and b.
As a first example of the use of the coordinate transformation given above let us present the perturbed redshift up to second order in PG in the presence of anisotropic stress and in function of the observer's angular coordinates [14]. This is easily done starting from the non-perturbative solution (3.3) and using Eq. (3.12) together with the fact that θ a actually are the standard angular coordinate at the observer position and are constant along the line of sight.
Therefore, any perturbative quantity can be written as a function of the observer's angular coordinates, i.e. the direction of observation, by Taylor-expanding it aroundθ a = θ a o . To this purpose, and considering general quantity with an angle-independent background value, it is enough to invert Eq.(3.9) to first order. Namely, we need the expansion To obtain familiar expressions in terms also of peculiar velocities, we recall that τ plays the role of a gauge-invariant velocity potential (see [13]), so that the spatial components of the perturbed velocity v µ of the (geodesic) observer in polar coordinates in PG are given by where τ (1) and τ (2) are the first-and second-order part of the coordinate transformation τ = τ (η, r, θ a ) given in Eq.
(3.7). The unit vector n µ along the direction connecting the source to the observer can be then expanded, in polar coordinates and to first order (which is enough for our purpose), Taking its scalar product with the spatial component of the perturbed velocity we obtain For the velocity components which are perpendicular to the radial direction we have v ⊥a v a ⊥ = γ ab 0 ∂ a P ∂ b P . (3.19) Let us also define the following useful variables which define the isotropic and anisotropic part of the Bardeen potential; ψ I is also the potential of the scalar part of the Weyl tensor [32]. In what follows we will use these variables to express our perturbed quantities. Furthermore, to simplify the notation, hereafter we only indicate the conformal time η as the integration variable along the line of sight, instead of all the arguments like (η, η o − η, θ a o ), of course the result also depends on the direction fixed by (θ 1 o , θ 2 o ). Finally, following [14], we obtain the observed redshift to second order and in function of the direction of observation Let us stress that the above integrals are along the line of sight, f (η ′ ) = f (η ′ , x i (η ′ )) and partial derivatives w.r.t. η ′ are only derivatives with respect to the first dependence. Above and hereafter we neglect the terms from the gravitational potential and the velocity at the observer position. They cannot be calculated within perturbation theory. Note, however, that in addition to the monopole and the dipole terms, they now also contribute to the quadrupole e.g. via terms of the form v 2 ∥o = (n ⋅ v o ) 2 . However, differences like C 2 (z 2 ) − C 2 (z 1 ) are meaningful and can be obtained from our formula, see [8] for a discussion of C 0 (z 2 ) − C 0 (z 1 ) for the power spectrum from linear perturbation theory. Furthermore, terms of the form X o Y s which are a product of a term X o at the observer position and a term Y s at the source position can contribute to all moments in the 4th order contributions to the power spectrum. Assuming that the correlation of X o and Y s can be neglected, for two sources at positions s 1 and s 2 they give correlators of the form ⟨X 2 o ⟩⟨Y s1 Y s2 ⟩. We have checked that X o is always either the velocity or the gravitational potential at the observer position. The first can be removed by going to the CMB frame while the second is of the order 10 −10 times the power spectrum and hence negligible. Similar considerations apply for the (reduced) bispectrum.

A. Redshift density perturbation
To second order the density perturbation at fixed PG coordinates can be written as We want to evaluate this at fixed observed redshift and with respect to the observer's angular coordinates θ a o . Therefore, we define a fiducial background model with coordinates (η for which the observed redshift and the past light-cone of our observer are given by We then expand conformal time and radial PG coordinates around the coordinates of the fiducial model as η s = η s , and obtain the terms of these expansions by perturbatively solving the following system of equations: where we use Eqs. We then obtain the perturbed conformal time and radial coordinate needed so that our measured redshift corresponds to the background redshift of our fiducial model, and where to second order we are allowed to drop the superscript (0) from η (0) s . Let us now expand the second order result (4.1) around the fiducial model, to obtain the second order density fluctuation at fixed observed redshift and direction of observation, Here all the quantities have to be evaluated at (η . Finally, using Eqs. (4.5)-(4.7) and the relation dρ dz = 3ρ 1+z , we obtain the following perturbative result: at first order and at second order where Eq. (4.13) agrees with the first order result of [6].

B. Volume perturbation
In order to compute the contribution to the galaxy number counts due to the volume perturbation, we start considering the 3-dimensional volume element seen by a source with 4-velocity u µ We then express this in terms of the observed quantities (z, The volume perturbation is determined by Let us now express the above quantities in GLC coordinates, x µ = τ, w,θ a with a = 1, 2, to show how the volume perturbation can be easily given fully non-perturbatively in GLC gauge. Considering thatθ a is constant along the geodesic of the light-coneθ a = θ a o , defined through w = const, we have Using also u w = −g wτ = Υ −1 and where we made use of Eq. (3.3) to calculate dz dτ . Note that Eq. (4.19) is again a non-perturbative expression for the volume element at the source in terms of the observed redshift and the observation angles in GLC gauge. To express the above quantity in Poisson gauge up to second order we have to calculate γ and dτ dz in this gauge. We first note that γ can be obtained, inverting Eq. (3.4), from the known result for the area distance given in [12][13][14]. In particular, neglecting the terms evaluated at the observer position, we have where we have set , the first and second order scalar perturbations of the luminosity distance, are given in Sec. III of [14]. Finally, from the results of Sec. III of [14] and after some algebraic manipulations, we obtain the following useful perturbative result to first order Here ∆ 2 denotes the angular Laplacian.
To second order we have a more involved result and This second term Ξ AS vanishes if no anisotropic stresses are present. Hereafter, we use the subscripts IS to denote terms in which the first order square contributions come entirely from ψ I = (ψ + φ) 2 and AS for terms which come from the anisotropic stress, On the other hand, to compute dτ dz to second order in perturbation theory in Poisson gauge we simply derive the coordinate transformation given in Eq.(3.7) with respect to the redshift after expanding it around the fiducial background (η . Proceeding step by step we first write The second order expansion of τ around the fiducial model is given by where now all the terms are evaluated at (η . Then, using Eq.(3.7), and defining the perturbative expansion of τ (η s , r s , θ a s ) = τ s , we obtain the following results We now derive the expressions above with respect to η (0) s . After a long but straightforward calculation (see Appendix A for some useful relations), we obtain the following results and Finally, inserting the results given in Eqs. , we arrive at the following expressions for the first and second order volume perturbation where (4.39)

C. Number Counts
We have now all we need to evaluate the galaxy number counts at observed redshift and in the direction of observation. Using the relations given in Eqs. (2.7) and (2.8) and inserting the results obtained for the density, Eqs. (4.13) and (4.14), and for the volume fluctuation, Eqs. (4.36)-(4.39), at fixed observed redshift and observer direction, we obtain to first order which agrees with [6], where it was derived directly in Poisson gauge.
To second order we obtain of course a much more involved expression: The contributions proportional to ψ A which vanish for ΛCDM are given by All these terms come from density fluctuations, redshift space distortions and lensing like terms. The ΛCDM expression corresponding to Eq. (4.45) can be easily obtained substituting ψ I with ψ (in this limit ψ I = ψ and ψ A = 0). The terms indicated in bold face are the once which we evaluate numerically in Section V for ΛCDM. The other terms, which are really new physical effects (like e.g. second order lensing), can also contribute significantly as they have four transversal derivatives We shall study them in detail in a forthcoming publication [33].

V. THE BISPECTRUM AND ITS NUMERICAL EVALUATION
In this section we compute two new contributions to the bispectrum of the number count fluctuations and evaluate them numerically for some configurations, namely the contribution from redshift space distortions and the one from lensing. Since we consider the truly observed galaxy number count, and not fluctuations on some unobservable spatial hypersurface, the 3-point function is in general a function of 3 directions 2 and 3 redshifts, The leading non-vanishing terms are in the form plus the permutations with respect to the galaxy positions (n i , z i ).
Expanding the direction dependence of ∆ in spherical harmonics, we can write where B m1m2m3 Statistical isotropy requires that B m1m2m3 where G m1,m2,m3 is the Gaunt integral given by On the second line we have expressed the Gaunt integral in terms of the Wigner 3j symbols, see e.g. [34]. This integral vanishes if 3 ∈ [ 1 − 2 , 1 + 2 ], if the sum 1 + 2 + 3 is odd, or if m 1 + m 2 + m 3 ≠ 0. The quantity b 1 , 2 , 3 (z 1 , z 2 , z 3 ) is the reduced bispectrum of the number counts. For Gaussian initial conditions it vanishes in linear perturbation theory and its first contributions are terms of the form ⟨∆ (2) (n 1 , z 1 ) Here we want to determine the presumably largest contributions to this expression. In subsequent work [33] we shall study it in more detail.
We write the perturbations in terms of Fourier modes and we express them in terms of the power spectrum of the primordial curvature perturbation R in (k). In this section, we assume simple adiabatic Gaussian initial perturbations from inflation which are given by the curvature power spectrum, For a given variable A we define the transfer function T A (η, k) by We will use also the angular power spectra defined as where P R (k) = k 3 2π 2 P R (k) is the dimensionless primordial power spectrum. ∆ A (z, k) denote the transfer functions in angular and redshift space for different sources. We will define the transfer functions explicitly for the cases analyzed in this work. The simplest relation is the one for the density fluctuations δ where j (x) is the spherical Bessel function of order . The largest terms in ∆ (1) are the density fluctuation δ (1) ρ , the redshift space distortion ∂ r v ∥s H s , and the lensing term In the following subsections we calculate some of their contributions to the bispectrum and compare the numerical results for some combinations of value of i and z i .

A. Density
We start by computing the leading term in the bispectrum, namely where the ⟨⟩ c denotes the connected part, and with (5.14) This expression for f 2 is valid for Ω M = 1 and Newtonian gravity [20]. We aim to estimate the order of magnitude of the different effects, leaving a detailed analysis to a future project [33]. So we adopt the approximation f 2 (k 1 , k 2 ) ∼ 1. This is of course not a very good approximation for certain angles of for very asymmetric situations, e.g. k 1 ≫ k 2 , but in this first analysis we just want to gain insight to the order of magnitude of the different terms. Approximating f 2 ≃ 1 significantly simplifies the calculation and still gives the correct order of magnitude for many cases. We defer the study of asymmetrical cases, like the squeezed limit, in future work [33]. Within this approximation, the second-order density perturbation is then given by Evaluating (5.12) we find where we have used the Wick theorem to compute the term into square brackets, We now expand the Fourier mode in spherical harmonics and Bessel functions Integrating over the anglesk 2 andk 3 and applying the orthogonality of spherical harmonics we find the three-point function The double integral in (5.20) is simply a product of two single integrals (due to the approximation f 2 (k 1 , k 2 ) ≃ 1 which we have adopted), so that we can simplify the result to b δδδ Here c δδ (z i , z j ) is the contribution to the number count angular power spectrum from density fluctuations. It can be calculated with the publicly available code CLASSgal described in [7]. We set the bias b between matter and galaxies equal to one, but it is easy to add a linear and scale-independent bias as outlined in CLASSgal.

C. Lensing
Finally, we also compute the contribution in the second order term which is given by the quadratic term combining density and lensing, see second boldface term in Eq. (4.43).
Repeating the steps for the previous terms and using that leads to With this we obtain the following expression for the reduced bispectrum containing one lensing term: This can be rewritten as where we have used (5.11) and (5.34).
Terms of a similar amplitude as the ones calculated here would be ⟨∆ (1RSD) (n 1 , z 1 ) δ (1) (n 2 , z 2 ) δ (2) (n 3 , z 3 )⟩ c + permutations and ⟨∆ (1L) (n 1 , z 1 ) δ (1) (n 2 , z 2 ) δ (2) (n 3 , z 3 )⟩ c + permutations, but for briefty, and since we expect them to yield similar results to the terms already evaluated here, we do not discuss them in this work. However, they have to be added when one wants a good estimate for the total amplitude of the bispectrum.

D. Numerical Results
We now evaluate and compare the amplitude of different configurations for the density, redshift space distortion and lensing terms. We use a ΛCDM model with cosmological parameters consistent with Planck [1] and compute angular power spectra with CLASSgal [7].
In Fig. 1 we plot the contribution to the bispectrum given by the density term b δδδ 1 2 3 , redshift space distortions b vδδ 1 2 3 , and lensing b Lδδ 1 2 3 by fixing two redshifts z 1 = z 2 = 0.8 while varying the third one around this value. Results are shown for different values of the multipole scale in the configuration = 1 = 2 = 3 2 (note that in particular with this choice 1 + 2 + 3 = even, as required by the Gaunt integral). Density and redshift space distortions show a peak when z 3 is equal to the other two redshifts. The density peak is symmetric, and also redshift space distortions give a nearly symmetric contribution (deviations are not visible by eye). In both cases, correlations vanish on large scales (large redshift differences) and can also lead to a negative bispectrum. At higher multipoles, the terms oscillate more and pass through zero for the first time already at smaller redshift difference ∆z. In Fig. 1 only the = 3 term is sufficiently well resolved to trust its oscillations. Note that the lensing term is typically two orders of magnitude smaller than the density and redshift space distortion terms.
On the other hand, lensing (which involves an integral along the line of sight) shows an asymmetric contribution with respect to z 1 = z 2 . Fig. 1 suggests that the higher values of the lensing term in the bispectrum are obtained for configurations where a third galaxies is placed at a redshift z 3 much lower than the other two galaxies. This comes from the fact that −c Lδ (z 1 , z 3 ) decreases with the redshift z 3 , simply since the density fluctuation decreases. Note however a small peak around z 3 = 0.8 = z 1 = z 2 , due to redshift auto-correlations of density-density contribution to b Lδδ 1 2 3 . Lensing gives a negative contribution, since its sign is given by the anti-correlation with density perturbations as already found, e.g., in [6,8]. Finally, from Fig. 1 it is already clear that redshift space distortions contributes to the bispectrum as much as density, while in this configuration lensing is subleading. Fig. 2 shows the different contributions to the bispectrum when considering three galaxies at the same redshift z while varying the multipoles = 1 = 2 = 3 2 scale. All contributions decrease with z. This is expected since we only consider deviation from a zero bispectrum coming only from non-linearities, which increase with time. There is a small turnover visible in the lensing bispectrum for ≥ 100 which comes from the suppression of lensing at small redshift. For = 3 this is not visible as the signal is dominated by very large scale structures. Not that this is the pure b Lδδ , not multiplied with any factors of which enhances the low signal.
In Fig. 3 we plot the contributions to the bispectrum again by fixing z 1 = z 2 = z 3 , but now as a function of the multipole scale 1 = 2 = 3 2 = . For the density and redshift space distortions the -dependence of the bispectrum is similar to the one of the power spectrum. A rapid increase up to ≃ 100 followed by a nearly scale invariant bispectrum, 2 b (z, z, z) ≃ constant. The lensing contribution however decays for ≳ 50. The oscillations in the lensing contribution (bottom panel) are numerical. They can be removed by enhancing the redshift resolution. Since this is numerically quite costly we have refrained from doing it in this order of magnitude discussion.
Any real experiment can only count galaxies within a certain redshift bin. Hence to compare our theoretical prediction with real observations we need to integrate the theoretical reduced bispectrum b 1 2 3 (z 1 , z 2 , z 3 ) over three window functions describing the redshift bins. Therefore we define the observed reduced bispectrum as where W (z, z ′ ) denotes the window function with center at z, normalized to the unity in the redshift bin, as a function of z ′ . Broad window functions smear out the signal along the line of sight. This typically leads to a reduction of the amplitude of the bispectrum. This reduction affects strongly the density and even more the redshift space distortions, if integrated over large redshift bins. The lensing, instead, which is a coherent signal integrated along the line of sight, is almost insensitive to the presence of the window functions. This implies that the signal coming from the lensing can become comparable to the density signal for wide window functions.
In Fig. 4 we shows the different contributions to the bispectrum integrated along a very large redshift bin, 0.2 < z < 3. We fix one multipole value 3 = 3 and vary = 1 = 2 . As expected, density and redshift space distortions signals decrease, and lensing is no longer a negligible term. This is consistent with what was found, e.g., in [8]. The redshift space distortion contribution to the bispectrum is less reduced than the one to the power spectrum since in the signal considered here we integrate over the product δ

VI. CONCLUSIONS
In this work we have computed the second order contribution to the perturbed galaxy number counts as a function of the observed redshift and the direction of observation. We have first expressed our result in the geodesic lightcone gauge [25] where it is simple, and exact for a given (exact) density fluctuation. Nevertheless, since so far the cosmological perturbation equations in this gauge are not known, we have translated our expression to the longitudinal or Poisson gauge. Our main final result is given in Eqs. (4.41), (4.43) and (4.44). Furthermore, in Eq. (4.45) we also present a simplified formula containing only the dominant terms. Since we have not used Einstein's equation in our derivation, the result is also valid for arbitrary models of non-interacting dark energy and for modified gravity theories as long as photons and dark matter particles move along geodesics.
We have then used our result to compute the probably dominant contributions to the angular bispectrum of the number counts, B 1 2 3 (z 1 , z 2 , z 3 ). In particular, we have calculated the terms coming from second order density, b (δδδ) 1 2 3 (z 1 , z 2 , z 3 ), from redshift space distortion combined with first order density terms, b (vδδ) 1 2 3 (z 1 , z 2 , z 3 ), and from lensing combined with first order density terms, b (Lδδ) 1 2 3 (z 1 , z 2 , z 3 ). As already with the power spectrum, b (δδδ) and b (vδδ) are typically of the same order, namely 2 b (δδδ) (z, z, z) ≃ 10 −4 − 10 −3 for small redshift z ≲ 0.5 and ≳ 100. The lensing term is about four orders of magnitude smaller. This results hold for a very narrow window function ∆z → 0 which we approximate as a Dirac δ-function. For broader window functions the density and redshift space distortion signals decrease while the lensing signal is virtually unchanged. If we integrate over redshift from z min = 0.2 to z max = 3 the lensing signal contributes close to 20% to the total and cannot be neglected.
This work mainly presents the derivation of the second order number counts expression and some preliminary evaluation of the bispectrum. In the future [33] we plan to analyze and compare the other terms in the bispectrum and to study whether the signal calculated here is measurable by planned surveys like Euclid. Especially, the new terms with four transverse derivatives are very promising and may well be at least as large as the lensing term discussed here. We shall also include the effects of biasing and magnification bias which have been neglected in this first discussion. The comparison of this bispectrum, which includes only the non-linear effects of gravity, with a possible primordial bispectrum will be necessary for any conclusions on primordial non-Gaussianities.

Acknowledgments
FM has enjoyed the hospitality of the Institute for Theoretical Physics at the University of Heidelberg during the completion of this work.
ED, RD and FM acknowledge support from the Swiss National Science Foundation. GM is partially supported by the Marie Curie IEF, Project NeBRiC -"Non-linear effects and backreaction in classical and quantum cosmology".

Appendix A: Useful relations
We give here some additional relations that have been used in the evaluation of the observed redshift and of the galaxy number counts to second order.
In particular to compute the second order redshift perturbation we need (see [14]): as well as