Relativistic wide-angle galaxy bispectrum on the light-cone

Given the important role that the galaxy bispectrum has recently acquired in cosmology and the scale and precision of forthcoming galaxy clustering observations, it is timely to derive the full expression of the large-scale bispectrum going beyond approximated treatments which neglect integrated terms or higher-order bias terms or use the Limber approximation. On cosmological scales, relativistic effects that arise from observing on the past light-cone alter the observed galaxy number counts, therefore leaving their imprints on N-point correlators at all orders. In this paper we compute for the first time the bispectrum including all general relativistic, local and integrated, effects at second order, the tracers' bias at second order, geometric effects as well as the primordial non-Gaussianity contribution. This is timely considering that future surveys will probe scales comparable to the horizon where approximations widely used currently may not hold; neglecting these effects may introduce biases in estimation of cosmological parameters as well as primordial non-Gaussianity.

Until very recently, any modeling of the galaxy bispectrum has relied on small-scales approximations, neglecting relativistic effects and radial correlations, given that LSS precise measurements were available only on relatively small scales. However, galaxy surveys are now at the point in which measurements of LSS are becoming extremely precise, and will be available on very large scales, therefore reaching the regime in which, as said above, a proper GR treatment is necessary and approximations such as the plane-parallel one fail. Hence, it is timely to have a full expression of the galaxy bispectrum in GR, avoiding approximations that are, in principle, inaccurate for correlations on very large scales. To this aim GR, bias and light-cone effects must be computed beyond linear order, at least to second order in the perturbations [52,54,55,63].
Various investigations on modeling the bispectrum including GR effects on large scales have been made recently. In particular, [49,50,52,53] compute the three-point correlation function in configuration space, including only some projection terms, but not the bias and relativistic perturbation solutions at second order (i.e. results obtained in [70,71,[98][99][100][101]); Refs. [54,55] analyzed the bispectrum in Fourier space including some bias terms but not nonlocal integrated terms (which [52,53] showed can be dominant for long radial correlations). In addition, the flat-sky limit is assumed even at very large scales, where this approximation has been shown to fail for the two-point correlation case [15,33,[102][103][104][105]. Of all the possible GR effects, it is well known that, at least on first-order statistics, magnification bias and wide angle corrections are the dominant ones [18,23,25,36,105]. These are also expected to be important for higher-order statistics, but other contributions may turn out to be comparable; a full numerical investigation for different surveys is needed in order to have more details on this. Given all these considerations, we extend previous works and derive, for the first time, the full wide-angle GR second-order expression of the bispectrum, i.e. including second-order bias, non-Gaussianity, velocities, Sachs-Wolfe, integrated Sachs-Wolfe and time-delay terms and lensing distortions from convergence and shear. This work builds on Refs. [46,47,50] where SW and lensing contributions to the second-order matter overdensity are presented. We envision that the resulting full expression of the bispectrum will be important for accurate estimations of primordial non-Gaussianity and cosmological parameters, taking into account all effects at very large scales, therefore fully exploiting the ultra-large-scale correlations that will be measured for the first time in the next decade. It will also provide an additional precise test of Einstein GR on those very large scales. Beside the effects we consider here, if the probability of observing a galaxy depends on its ellipticity, there is a selection effect beyond magnification, which is affected by shear and intrinsic alignment (see, for example, [106,107]). Here for simplicity we are not considering this orientation-dependent selection effect which will produce new terms both at first and second order that depend on the large-scale tidal field. The inclusion of this effect is left to future work.
We improve upon previous works in the literature in the sense that for the first time we write down the full expression in Fourier space, and we do so using what we call the "spherical-Bessel" formalism, which was first developed by [108] and, subsequently, in redshift space by [109] for the two-point correlation function. (Let us point out another interesting approach studied in [110,111], where they compute bispectra with total-angular-momentum waves.) The expressions are inevitably lengthy and involved but whenever possible we give a physical insight on their meaning and we summarize schematically the structure of the full expression in Section II. The expressions provided here will not be directly relevant to data analyses, their value is in providing a starting point to devise useful approximations, and a quantitative evaluation of the different contributions. The next natural step is to derive a more compact expression including the dominant terms, for different geometrical configurations (shapes), which will be presented elsewhere.
The rest of the paper is organized as follows. We begin by presenting the second-order number counts on the light cone in Sec. III. This is a key ingredient in the expression for the bispectrum which we present in Sec. IV, which expression depends on several kernels that are written down explicitly in Sec. V for the first-order, and Sec. VI at second order. The second-order bias is derived in Sec.VII. We conclude in Sec.VIII. The appendices report useful relations and expressions needed to follow the derivation of our main results and link our analytical expressions to numerical codes available in the literature 2 into our expressions.

II. EXECUTIVE SUMMARY
Here we summarize in a concise way the main result of this work, showing the main expression for the spherical-Bessel bispectrum, including all general relativistic, local and integrated, wide-angle and mode-coupling terms, that we call S 3 ({k m} 1 , {k m} 2 , {k m} 3 ), where {k m} i = (k i , i , m i ) for i = 1, 2, 3. We write the galaxy bispectrum using the "spherical-Bessel" formalism first proposed by [108] and, subsequently, in redshift space by [109] (see also [113]) for the two-point correlation function, and then by [114] for the three-point statistics. The same formalism (for the power spectrum) has been applied to real data in [115] and extended to include GR effects in [17]. [See also Ref. [116] where they discuss the limit of various approximations (e.g. flat-sky, limber) which are applied to the lensing signal.] Here we generalize this formalism to include all wide-angle and relativistic terms, and write down the full expression of the galaxy bispectrum, without any approximations besides second-order perturbative expansion in the perturbations and the bias and without neglecting any contributions (apart from vector and tensor perturbations). In this section we report the general expression for S 3 ({k m} 1 , {k m} 2 , {k m} 3 ), and in the rest of the paper we show its derivation and the explicit expression for all the terms. Here and hereafter: where ∆ g is the observed galaxy fractional number overdensity (for further details, see Sec. III) and contains the contributions of all the local and integrated terms including bias, and the indexes t, u, v indicate the order, so that in the right-hand side two of them will be 1 and one will be 2, cyclically. To be precise S 3 is actually the spherical 3-point correlation function of Fourier-space galaxy overdensity (but, for simplicity, we will refer to it as the bispectrum); in the same way as for the power spectrum, the correlation of the field δ in Fourier space δ(k 1 )δ(k 2 ) is related to the power spectrum P (k) via a Dirac delta function δ(k 1 )δ(k 2 ) ∝ P (k 1 )δ D (k 1 + k 2 ).
As an example, let us consider below only the first additive term of Eq.(1), where t = 2 u, v = 1 (the other terms will be just permutations of this). Since each ∆ g(t) m (k) is a sum of several components, it will be a sum of several contributions which we indicate by the running indices a, b, c. Each of these indices will model the effect of all the physical quantities we consider: where δ refers to the matter overdensity (intrinsic clustering), "rsd" and "v" are velocity (peculiar velocities (RSD) and doppler) terms, L contains lensing terms and "pot" gravitational potentials (ISW and STD). Each of these quantities enter (alternatively) at the first and second order. Within our formalism, in this specific example, the index a is associated with a second-order quantity and refers to the terms expanded in Eq. (39), while the indices b, c are associated with first-order quantities and refer to the terms that we will make explicit in Eq. (31). Therefore, we have Explicitly (for derivation details, see Sec. IV C and Appendix D) the bispectrum building blocks Υ abc 1 m1 2 m2 3m3 are functions that contain all of the information needed to express the observed bispectrum: 2 In Appendix E 2 b, we give a possible prescription on how to insert the numerical outputs obtained in SONG [112].
where the sum is over p 1 m p 1 q 1 m q 1¯ 1m1 , and Here we have used the Gaunt integral The spherical multipole functions M contain all the contributions (at first-indicated by superscript (1)-and secondsuperscript (2)-order) for the density, velocity, lensing and gravitational potentials, and their combinations; their explicit expressions will be presented below. Here P Φ (k) is the primordial linear power spectrum of the Bardeen gravitational potential; it will be clear below why it is advantageous to express the second-order GR bispectrum for number counts in terms of a primordial quantity, but we will anticipate it briefly here. The primordial linear Bardeen potential is the only relevant spatial field that has no gauge issue, that-by construction-remains invariant in time, is statistically homogeneous and isotropic and thus for which it makes sense to perform a three-dimensional Fourier transform. In fact, for any quantity that evolves along the line of sight or may be affected by projection effects such as lensing or ISW, or for which the flat-sky approximation does not hold, there is no unambiguous three-dimensional Fourier transform. It is also useful to provide a relation between the spherical Bessel bispectrum S and the bispectrum in Fourier angular space, B 1 2 3 (k 1 , k 2 , k 3 ): and between the 3-point function and the bispectrum in Fourier angular space: after resumming over 3 As usual, Y m here denote the spherical harmonics and m 1 m 2 m 3 the Wigner 3-j symbols. In our case, the bispectrum in Fourier angular space is written as . Thus all local and integrated effects due to, for example, the bias (see Sec.VII), lensing, SW, gravitational evolution, non-Gaussianity etc., are enclosed in the spherical multipole functions M. [Precisely, the generating functions are 3 The denominator of the first term in Eq. (7) reads ℵ * (k 1 )ℵ * (k 2 )ℵ * (k 3 ) [see Eqs (71) and (80)], with the coefficients being ℵ = k (2/π) 1/2 in the formalism of [17,114], while they become ℵ (k) = 4πi in the total angular momentum formalism of [110,111] (for details, see Appendix B). { 1 m 1 , 2 m 2 , 3 m 3 }. Here and hereafterk denotes the angular position on the unit sphere of the corresponding unit vector.

III. SECOND-ORDER NUMBER COUNTS ON THE LIGHT CONE
In this section we start presenting all GR effects that have been computed previously by 4 [50]. We begin by assuming a concordance background model, and at first-order we neglect anisotropic stress, vector and tensor perturbations. In the Poisson gauge, the metric and peculiar velocity expanded to second order are [70] where η is the conformal time, a is the scale-factor, and we have omitted the superscript (1) indicating terms of the first-order expansion on familiar quantities such as the metric perturbation Φ, Newtonian potential Ψ and the galaxy peculiar velocity ∂ i v. Here and hereafter, second-order terms and indicated by the superscript (2). Hereω [70,[117][118][119][120][121][122]. Redshift-space or redshift-frame is the "cosmic laboratory" where we probe the observations. In redshift-space we use coordinates which effectively flatten our past light-cone so that the photon geodesic from an observed galaxy has the following conformal space-time coordinates [14,47,123]: Hereχ(z) is the comoving distance to the observed redshift, calculated in the background (i.e., a redshift-space quantity),n is the observed direction to the galaxy, i.e.n i =x i /χ = δ ij (∂χ/∂x j ) . Usingχ as an affine parameter in the redshift frame (at zeroth order), the total derivative along the past light-cone is d/dχ = −∂/∂η +n i ∂/∂x i . In our analysis we use only the observed redshift z rather than the background (Hubble flow) redshift. In particular all background quantities are not evaluated at the background, redshift (i.e. the redshift that would have been observed for the same source without any perturbations along the line of sight), they are instead evaluated at the observed redshift which include real world effects such as peculiar velocities. Defining x µ (χ) as the coordinates in the physical frame [see Eq. (9)], where χ is the physical comoving distance of the source, we can set up a mapping between redshift space and real space (the "physical frame") up to second order in the following way: x µ (χ) =x µ (χ) + ∆x µ(1) (χ) + ∆x µ(2) (χ)/2.
In this work, ⊥ denotes projection into the screen space (with projector P ij = δ ij −n inj ), indicates projection along the unit line-of-sight vectorn i , and we define the derivatives Now we want to study the physical number density of galaxies n g as a function of the physical comoving coordinates x µ and the magnification M. In particular we consider the cumulative physical number density sample with a flux larger than a observed limitF which can be translated in terms of a the inferred threshold luminosityL(z) (n g = N in [13,19]). The physical number density contained within a volume V is given by where n g (x α , M) is the physical number density which occupies the comoving physical volume dV, a is the scale factor, √ −ĝ = √ −g/a 4 ,ĝ µν is the comoving metric. In the redshift frame is, by definition, where the observed comoving volume is d 3x = dV. Then, relating the observed galaxy number density with the physical one, i.e. Eqs. (13), (14), we obtain the observed fractional number overdensity where and Here we have considered the corrections up to second order of the volume: the scalar factor: a(x 0 (χ)) =ā 1 + ∆ ln a (1) whereā = a(x 0 (χ)) = 1/(1 + z), and .
Finally, we used where ∆n g (x α ,L) (1) and ∆n g (x α ,L) (2) contain also the fluctuation of the luminosity distance at the observed redshift z (e.g. see [47]). In order to write explicitly all above relations, first of all, let us define at first-order the following quantities: • the scale factor correction 5 where I (1) is the integrated Sachs-Wolfe (ISW) effect at first-order, i.e.
where, here and hereafter, the prime denotes ∂/∂η; • the weak-lensing convergence term • the weak-lensing shear γ (1) ij and rotation ϑ ij terms, where where∂ i = ∂/∂x i ; • the radial displacement at first order that corresponds to the usual (Shapiro) time-delay (STD) term [13] • we define the following 3D nonlocal vector [46] S which can be split as The physical meaning of this vector can be understood by noting that the transverse part, i.e. S i(1) ⊥ , can be related with T (1) and κ (1) in the following way −2κ (1) • Finally the overdensity wheren g = n g (x 0 ,L) (0) is the background number density of sources with luminosity exceedingL and n (1) g = n g (x α ,L) (1) .
Expanding the galaxy fractional number overdensity, at linear order we find [13,14,46,47,50] ∆ (1) where b e = ∂ lnn g (ā,L) ∂ lnā + 3 is the evolution bias term related to the background number density, and H = a /a is the conformal Hubble factor and we have defined the background magnification bias as All the quantities defined at linear order above, at second order become [46,47,50] ∆ ln a (2) where A (i B j) = (A i B j + B i A j ) /2, and 1 where the density contrast at second order has been defined as and for convenience of notation we sometimes use n (2) g = n g (x α ,L) (2) . Here we have defined the first-order magnifi-cation bias in the following way Let us point out that, using the result from Ref. [50], we have indirectly assumed that galaxy velocities follow the matter velocity field, i.e. no velocity bias, and we have used the velocity equation both at first and second order (for example see Appendix A). The readers interested in the primordial non-Gaussianity contribution (or f NL for short) should keep in mind that it is implicitly enclosed in second-order primordial quantities such as Φ (2) ; see e.g., Sec. E 2 a. Now, replacing Eqs. (22) in Eq. (31) we have the galaxy density contrast at first-order that it is given by 8 terms involving δ g , first and second derivatives of the velocity, potential, convergence, ISW and STD: At second order, replacing Eqs. (22), (34) and (35) in Eq. (39), we have where we have separated the second-order contribution to ∆ (2) g in local (subscript loc) and integrated (subscript int) terms. Here the local contribution is the sum of three terms contains convergence, ISW and STD terms. In we have terms like ∂ ⊥i A (1) ∂ i ⊥ B (1) , where A (1) or B (1) can be a local or an integrated term at first-order. In we find integrated terms such as, for example, ISW and STD at second order. In the fourth term, we identify contributions in which the product between a local and integrated term (or two local terms) is within another integral along the line of sight. Finally in the fifth term we find all symmetric trace-free terms with orthogonal partial derivatives. Obviously, Eqs (44) and (48) are written according to the properties of the various local and integral terms. Hereafter for simplicity we will compress all these equations by writing where Eq.(54) is given by the sum of all terms contained in Eqs. (45), (46), (47), (49), (50), (51), (52) and (53).

IV. BISPECTRUM
The spherical Bessel representation uses a complete set of orthogonal basis functions |k m in a spherical Fourier space (for a very brief review of this formalism, see Appendix B). A scalar field like ∆ g in configuration space can be decomposed in the following way: It is important to note that by definition the monopole ∆ g 00 can be removed and set to zero for k → 0. In general, we can discard ∆ g 00 because only the meann g (z) contributes to the monopole. Therefore we will compute the spherical power spectrum only for > 0.
Here we introduce a more realistic definition of x|∆ g = ∆ g (x) which in a realistic case becomes: where we have included the radial selection function W(χ).
Then at first and second order we find where the index b runs over the terms defined in Eq. (42) and a represents the index of summation over all additive terms at second order, see Eq. (54). Finally, ∆ g(1) m (k) is the spherical Bessel transform of Eq. (42) and the second-order terms are the spherical Bessel transforms of Eq. (54).

A. First-order terms
Considering each term of the first summation on the right-hand side of Eq. (57), we have where we have already included the radial selection function. Usually, we can define a generic first-order perturbation as follows where Φ p (k) is the primordial potential set during the inflation epoch, T b (k, η) is a generalized transfer function which relates the linear primordial potential with a generic perturbation term (here labeled with b). Here k is a Fourier space vector,n is the observed galaxy unity vector on the sky, and we have used the following relations where ℵ (k) is a function that depends on and k (see Appendix B) which specific expression depends on the conventions adopted in the expansion. For each contribution in Eq. (42), we have defined a weight function W b which is a generic operator that depends onχ, η, ∂/∂χ, ∂/∂η and n . This operator encloses the physical effects due to the fact that in any observation we collect the photons emitted from a source after they have traveled through the past light cone of the observer. It is this "projection effect" on the perturbations that is captured by W; W takes slightly different explicit expressions depending on what perturbation is being considered, as labeled by its superscript index. Applying it to the spherical harmonics, we find we have removed its angular dependence. Using Eq. (60), we have (61) Then, taking into account that where δ K denotes the Kronecker delta and * denotes complex conjugate, we find where Before concluding this subsection, it is useful to calculate the spherical power spectrum (see also [17]), i.e.
where P Φ (k) is the power spectrum of the potential at initial epoch

B. Second order (scalar case)
In the same way, at second order, the spherical Bessel transform of each term contained in the last summation of Eq. (57) can be written as In general (for the scalar case), similarly to what was done at first-order, let us define the spherical multipole functions M (2) , at second order, in the following way: In addition, as we will show explicitly in Sec. VI (see also Appendix F), M can be expanded as where we used the following theorem that relates Legendre polynomials P¯ to spherical harmonics and we obtain As we will see in Sec. VI and Appendix F, M a(2) m p mp qmq¯ m (k; p, q) will be computed explicitly for all terms contained in Eq. (54). Finally, using Eqs. (63) and (70), in the next section we derive the expression for the spherical Bessel bispectrum.

C. Master relation for bispectrum (scalar case)
In order to compute the projected bispectrum in spherical Fourier space, let us start with the following relation 6 where we used See also Eq.(B18). In particular For simplicity, let us consider only the first additive term of Eq.(72) where the index a identifies all the terms at second order in Eq. (54), and the indexes b and c identify all the terms in Eq. (42). For a, b and c fixed we have Taking into account that and using where Υ [j]abc 1m1 2m2 3 m3 , for j = 1, 2, 3, are the bispectrum building blocks and contain all the information on the bispectrum in redshift space. Explicitly, we have and Υ . For symmetry reasons one may expect that In Appendix D will show explicitly that this is the case. The explicit proof offers a consistency check on our calculations and expressions. In Appendix D we will also prove that This implies that we can discard Υ [1]abc because here we consider only terms with > 0. Finally we will demonstrate that the three-point function of the coefficients defined in Eq.(72) may be factorized by isotropy 7 into where 1 + 2 + 3 =even, m 1 + m 2 + m 3 = 0 and 1 , 2 , 3 satisfy the triangle rule. From the above relation and Eq.(71) we can write one of the main results of this paper, i.e.
The rest of paper is devoted to computing explicitly all M a(2) m p mp qmq¯ m (k; p, q) in Eq. (54) and in Sec. VII we give a prescription to incorporate consistently the bias at second order within δ

VI. SPHERICAL MULTIPOLE FUNCTIONS AT SECOND ORDER
In this section we will apply the instructions contained in Sec. IV B where we proved that a generic (scalar) term at second order can be expressed in terms of M a(2) m p mp q mq¯ m (k; p, q). In the next subsections we will consider separately each additive term contained in Eqs. (44) and (48), i.e. ∆ (2) g loc−i for i runs from 1 to 3, and ∆ (2) g int−j where j runs from 1 to 5.

A.
Terms from ∆ Equation (45) contains all local contributions with second-order perturbation terms. These additive terms can be written in the following way: where W α has been introduced in Eq. (59), and the superscript α[b] means that we are taking combinations of different explicit "projection effect" expressions W, labeled by the superscript index α, with ∆ labeled by the superscript index b.
Here ∆ b(2) = δ (2) g , Φ (2) , Ψ (2) , v (2) , respectively, for b = δ (2) g , Φ (2) , Ψ (2) , v (2) and where F b(2) (p, q,k; η) are specific kernel functions (fully symmetric functions) and depend also on the parameters related to primordial non-Gaussianity. Explicit expressions for F b(2) (p, q,k; η) are reported in Appendix E 2. Using Eq. (83) we find Here, by construction, we note that, in these kernels,k is related to p, q, and the angle between them, cos(θ pq ) viã , η] and in general we can expand in Legendre polynomials the dependence of the kernel on θ pq using Hence the second-order kernels can be expanded as In Appendix F we present an alternative way to compute these terms. This new method might be important if F b(2) (p, q,k; η) is function ofk n , where n < 0.
Finally we find Below we show explicitly all the terms in Eq. (45): All local projection terms with time and partial space derivatives along the line-of-sight direction of ∆ (2) g loc−2 in Eq. (46) can be written in the following way: where the superscript ij means that we are taking combinations of different g . [Note that in Eq. (57) the index a also runs through all combinations of i, j.] Applying the weight function W j on ∆ j(1) we have where we have used Eq. (60). Starting from Eq.(66) and using Eq.(5), we find Then, explicitly, we have Replacing the indices ij by the corresponding symbol for each additive contribution in Eq.(46)-and thus making more explicit their physical meaning-below we list all M ij m p mp qmq¯ m : C. Terms from ∆ In ∆ (2) g loc−3 [see Eq. (47)], we find all local projection terms with transverse partial derivatives and we can write these terms as where here the indices r, s, and rs have the same meaning of the superscript i, j and ij of previous subsection, ∆ s(1) was already defined in the previous subsection, and Starting from (see Appendix C) where we used we find Here (2) ∇ i is the covariant derivate on the unit sphere, ð (ð) is the spin raising (lowering) operator and s Y m are spin-weighted spherical harmonics (for more details, see Appendix C). Using where we have defined we have i.e. we can sets = 0 (because s 1 + s 2 +s = 0 wheres = s 3 ). This implies that we can project ∆ a(2) (x) on |k m space and finally obtain Here we may further simplify the notation in the following way. We note that Given that p + q +¯ is even, we can use the following identity [149] p q˜ p q˜ 0 0 0 and defining [134] I mm1m2 we can rewrite Eq. (139) as and, explicitly, we have Applying Eq.(142) for each term in Eq.(47) we find D. Terms from ∆ or (ii) (147) For the case (i) we obtain and, using the Gaunt integral Eq.(5), ∆ j(2) (x)/2 projected in |k m space turns out and the transfer function can be written in the following way Likewise, for (ii), we have Explicitly, we list the terms computed from Eq. (49) E. Terms from ∆ Here we have terms like ∂ ⊥i A (1) ∂ i ⊥ B (1) , where A (1) or B (1) can be a local or an integrated term at first order. From ∆ (2) g int−2 we can write these terms as (i) or (ii) where here the index a = rs. Following the prescription in Sec. VI C we obtain Then, for (i) we find and, for (ii), where we have directly used the following identity: where Im mpmq p q has already been defined in Eq. (140). Projecting these contributions in |k m space, we find, for i), and, consequently, Finally, for (ii), we have and M rs(2) Here below we list all transfer spherical multipole functions for cases (i) and (ii) : F. Terms from ∆ In ∆ (2) g int−3 we find integrated terms as, for example, ISW and STD at second order. In particular, for the first two additive terms of Eq. (51), we can use the following expression In this case, applying the same prescription used in Sec. VI A, we find quickly and Using the same approach, the last two additive terms of Eq. (51) can be written together in the following way From Eq. (133), we note that and using Eqs. (5) and (171) we find Let us conclude by noting that, in Appendix F, we compute also these terms in a different way.
In ∆ (2) g int−4 , each contribution is an integral along the line of sight of the product between a local and an integrated term (or two local terms). Specifically, we have the following possible terms (i) (ii) and (iv) Following the approach used in the previous subsection we find, for (i), and, therefore, For (ii), we obtain and, consequently, Instead for (iii) we have and so Finally, for (iv), it turns out and hence Here below we write explicitly all terms contained in Eq. (52): H. Terms from ∆ We immediately see that ∆ (2) g int−5 contains all symmetric trace-free terms with orthogonal partial derivatives. In order to compute correctly all these relations, it is useful to rewrite∂ i ⊥∂ j ⊥ with covariant derivates (2) ∇ i and, then, the usual spin-raising and spin-lowering operators (for more details, see Appendix C). In particular, we havē where and using P ij = 2m +(i m −j) , we find Now taking into account Eq. (134) and Eq. (229) turns out (234) Using these results we can focus on the following terms: (i) and, (ii) a shear term related to post-Born term perturbations (236) For (i) and using Eqs. (25) and (232) we find where that the CS gauge is entirely appropriate to describe the matter overdensity at second order (for details, see [63]). The bias should be defined in the rest frame of CDM, which is assumed to coincide with the rest frame of galaxies on large scales and can be computed using the peak-background split approach [126]. In CS gauge, the spherical collapse model has an exact GR interpretation [70,98,127,128]. Indeed, only in this frame we can resort to the peak-background split [124,125] approach, in which, in the Press-Schecter-inspired prescription [126], halos of a given mass, M, collapse when the linearly growing local density contrast (smoothed on the corresponding mass scale) reaches a critical value. This is important for a self-consistent calculation of the so-called non-Gaussian halo bias and for precision parameter estimates, introduced by nonlinear projection effects. In ΛCDM, the CDM rest frame is defined up to second order in the CS gauge-in which the galaxy and matter overdensities are gauge invariant. The CS gauge is defined by g CS 00 = −1, g CS 0i = 0 and v i CS = 0. At second order, = 0 (where, for simplicity, we have removed the purely transverse space-dependent constant C ⊥ i in g 0i ; see [129]). Before proceeding into the main part of this section, some comments are in order.
(i) The primordial non-Gaussianity (NG) has to be considered in all second-order contributions. NG is set at primordial inflationary epochs on large scales. At later times cosmological perturbations reenter the Hubble radius during the radiation or during the matter epoch. Beyond linear order, integrated (or projection) effects couple large and small scales. Thus the NG large-scale information leaks into smaller scales. Therefore a complete general relativistic computation is required to evaluate all the observable imprint of PNG in the LSS.
(ii) The long-mode curvature perturbation modulates the matter overdensity, with an effective f GR NL = −5/3 [64], but this effect cancels out in the halo overdensity, when perturbations are evaluated at a fixed local scalex, rather than fixed global scale x. More in detail, the small-scale density at a fixed local physical scale is independent of the long-wavelength perturbation. Thus the long-wavelength mode has no effect on the small-scale variance of the density field smoothed on a fixed mass scale. In other words, the long-wavelength perturbations are not observable locally if the distribution of the primordial metric perturbation 8 ζ is a Gaussian random field and we conclude that, within δ (2) mCS , f GR NL = −5/3 is, in the strictly squeezed limit, reabsorbed via a local coordinate transformation (see [64,130,131] and Appendix E 2 a).
Taking into account the above discussion, we define the scale-independent bias in a completely general way at first and second order (considering scales down to the mildly nonlinear regime) as [56,59,60,62,65,66] depend on the conformal time 10 , ζ is the comoving curvature 8 Here the comoving curvature perturbation ζ has the same definition and notation as in e.g. Refs. [58,100,131]. Finally ζ here is −ζ in [121], and ζ = −Rc defined in [99]. 9 Let us note that if the single-field consistency conditions hold, then for example we have b 01 = 0. 2 should be also functions of space and time. Then, in Fourier space, δ gCS (k) should be written in the following way: perturbation and Let us point out that δ (2) gCS depends also on primordial non-Gaussianity (see Appendix E 2 a for an analysis at k ≤ k eq where the primordial f NL is written explicitly). Now we need to connect δ gCS in Eq. (251) with δ g defined above in the Poisson gauge. From [46,47,50] we know that Note the useful relation 11 v = χ (1) /2. In particular, Eq.(256) can be expressed as follows where we used the relation Bearing in mind Eq.(251), Eqs. (255) and (257) become Now, using this definition for the bias in Eqs. (251) and (255), the transfer function T δg , defined in Eq. (82), can be explicitly written as Finally we can also rewrite the magnification bias, which first appeared in Sec. III, as and, consequently, T Q (1) , in Eq. (128), turns out to be (262) Appendix E is devoted to write explicitly all the transfer functions T b (k, η) and the kernels F b(2) (p, q; η) for the (spatially flat) ΛCDM model (all these results can be quickly generalized for CDM+dynamical dark energy models).

VIII. CONCLUSIONS
In this paper we presented, for the first time, the full expression of the galaxy three-point function (the bispectrum) including all relativistic local and integrated terms at second order in the perturbations, in the wide-angle geometry (we never rely on the Limber and flat-sky approximations) and including the non-Gaussianity parameter f NL . We also include the galaxy bias up to second order, including a tidal term and add all magnification corrections. We believe that the calculations and expressions presented here, despite being cumbersome and appearing tedious and complicated, are a valuable and timely contribution. Future and forthcoming galaxy surveys will map the sky on ultralarge scales with unprecedented precision bringing LSS in regime where the simple Newtonian approximation is not longer sufficient. The next-to-leading-order correlation (such as the bispectrum) despite being a lower signal-tonoise quantity than the power spectrum, encloses key complementary physical and cosmological information, which should not be neglected if we are maximizing the scientific return from this observational effort. Given that the bispectrum on large scales will be measured, the correct interpretation of the measurement relies on accurate and precise theoretical modeling of the signal. It is for this reason that we believe it is of fundamental importance at this point to have the full expression of the galaxy bispectrum, and we embarked in this challenging and time consuming task, despite its apparent complexity.
Working in spherical Bessel coordinates, we derived a compact expression for the bispectrum that encompasses all the physical contributions at first and second order, in curved sky and including integrated terms for radial configurations. We found that we can write the full GR bispectrum as after resumming over { 1 m 1 , 2 m 2 , 3 m 3 }, U is a scale-dependent angular function and T SH i,mi are tripolar spherical harmonics; see Eq. (80).
In this case, the bispectrum in Fourier angular space is written as where C j is a combination of multipoles and 3j-and 6j-Wigner symbols, P Φ denotes the primordial (linear, Gaussian) power spectrum of the Bardeen gravitational potential and K a(2) It is evident that the expressions provided here will likely be impractical for most realistic applications or data analyses. Hence approximations will have to be made. The results presented here provide the starting point to devise suitable approximations and a reference and benchmark to assess the validity, accuracy, performance, advantages and disadvantages of any such approximation. For example for numerical and computational reasons, real data analyses should involve an optimal choice of which terms to include. We envision that the results presented here will be a reference for assessing which physical effects are most important and which one are ultimately negligible. The choice will depend on specific survey characteristics, on the physics to test, models considered or question at hand. Significant work remains to be done in evaluating the relative importance of the (many) different terms and physical contributions to the overall signal in different regimes. While preliminary investigation suggest that this is a viable program, we leave this effort to future work. and At the moment, ℵ (k), (k) and ‫ג‬ (k) are generic functions. Obviously, as we will see below, these functions are closely related to each other. The completeness conditions on all possible bases are defined in the following way: Here, by construction, (k) and (k) are real, i.e.
and using the completeness conditions, we find For φ m (k): Relations between ℵ (k), (k), ‫ג‬ (k) and (k): i) From x|k we find From the above relation, immediately, we obtain the following property where we have used Then we find In the literature, these coefficients are usually fixed in the following way: where in (1) see Refs. [17,114] and in (2) see Refs. [2,110,111].
In the paper it is also useful to write the inverse of Eq. (B10). Using the results obtained in this section, we immediately find Appendix C: Covariant derivative on unit sphere and the spin-weighted spherical harmonic decomposition In this appendix, we outline the decomposition of three dimensional quantities on a unity sphere in the observer tangent space and show how one can derive the relation between the covariant 2D derivative and spin-weighted spherical harmonics function. This discussion is based on Refs. [123,[132][133][134][135][136][137][138].
Starting from the orthonormal polar basis vectors {n, e θ , e ϕ }, it is possible to define the following helicity basis which are also called spin±1 unit basis vector on the unity sphere. (In the literature, e.g. Refs. [135,136], m ± is also denoted as m + =m and m − = m.) Under a right-handed rotation of the coordinate system {e θ , e ϕ } around n by an angle ψ, so that e θ → e θ = cos ψ e θ + sin ψ e ϕ e ϕ → e ϕ = − sin ψ e θ + cos ψ e ϕ , m ± transform as m ± → m ± = exp (±iψ) m ± .
In particular, m ± satisfy the following properties Given a spin-weight s complex function s f (n), it transforms as s f (n) → e −isψ s f (n). Then, defining the projected tensor as ..m is − transforms as a spin-weight s object. Specifically, if T j1...js is tensor symmetric and trace-free component of a rank s tensor, we can associate it with a s-spin weight function in the following way At this point it is important to define the covariant derivative on the unit sphere 12 1 χ (2) ∇ i T j1...js = P p i P q1 j1 ...P qs js ∂ p T q1...qs .
From the above relation we note immediately that, and, for s = 0 we have (Here 0 f = f .) Instead, with two derivates we find Now it is useful to relate (2) ∇ i to the usual spin-raising and lowering operators, i.e. ð andð, defined and analyzed for the first time in Refs. [139][140][141][142][143][144] (see also e.g. [145][146][147][148]). First of all, the spin operators are defined in the following way and have the ability to transform under the angle ψ rotation as Then the covariant derivate of T j1...js can be written as follows: where Conversely, for s ≥ 0, and we have to replace m is − with m is + for s < 0. In order to describe correctly objects of spin-weight s we need to generalize the spherical harmonic basis Y m (θ, ϕ). This new basis is called spin-weight spherical harmonics s Y m (θ, ϕ) and it can be related with the usual Y m (θ, ϕ) = 0 Y m (θ, ϕ) in the following way: and satisfies the properties Then we findð and we can rewrite Eq. (C7) by using the spin-s spherical harmonics, i.e.
Now we have all the required tools to connect quickly the covariant derivative on the unit sphere with the spin-s spherical harmonic basis. For example, and Let us conclude this section with the orthogonality and completeness expression and the generalization of the Gaunt integral where I s1s2s3 1 2 2 has already been defined in Eq. (137).
In Appendix F we will computeF a (p, q,k; η). With this new approach we will not need to compute σ 2 (p, q).

b. Prescription for SONG
SONG 16 is an open-source second-order Boltzmann code which includes all the effects of metric, CDM, baryons, photons and neutrinos, see e.g. Ref. [112]. From this code one can compute all possible kernels in Fourier space from the radiation era, both at large and small scales; the addition of the non-Gaussianity parameter f NL is straightforward following our prescription.
Taking into account that SONG is written in the Poisson gauge, i.e.