The skewness of the distance-redshift relation in ΛCDM

Starting from a recently proposed framework for the evaluation of the cosmological averages, we evaluate the higher-order moments for the distribution of a given observable. Then, we explicitly discuss the case of the Hubble-Lemaître diagram and evaluate its skewness at the leading order in the cosmological perturbative expansion of the gravitational potential. In particular, we focus on perturbations of the luminosity distance due to gravitational lensing. Finally, we discuss our findings in view of recent numerical relativistic simulations, confirming that the skewness in the Hubble-Lemaître diagram primarily originates from the late-time matter bispectrum, with other line-of-sight projection effects being sub-dominant.


Introduction
Luminosity-distance relations have been crucial in establishing the accelerating Universe [1,2] and the standard cosmological model ΛCDM.As observational cosmology has progressed, it has moved from being an order-of-magnitude science to one of percent precision.To make effective use of the wealth of unprecedented data now [3][4][5] and in the coming years [6][7][8][9][10][11], it is essential that our theoretical tools evolve accordingly.With regard to the luminosity-distance relation, a key focus is to understand how light propagates in an inhomogeneous Universe and how such effects give rise to observable features.This requires a deeper understanding of the interplay between cosmological structures and the propagation of light, enabling us to identify the subtle signatures imprinted in the observational data.
While the CMB is very well described by linear physics, the study of the Large-Scale-Structure (LSS) of the Universe requires going beyond the linear approximation to extend the range of scales that we can use to constrain cosmological parameters.This can be achieved either by analytical perturbative approaches or non-perturbatively by numerical simulations.In both directions there has been a tremendous improvement in the last decade.Although numerical simulations can, in principle, handle fully non-linear processes, analytical approaches are computationally much simpler and provide valuable insights into the underlying phenomena.In particular, the two approaches can be compared within the range of validity of perturbation theory.When describing real observables, it is crucial to consider that discrete sources (galaxies, supernovae, quasars, etc.) are observed through the collection of photons emitted, travelled and detected in a clumpy Universe.Over the last few decades, considerable effort has been devoted to providing a comprehensive relativistic framework for describing LSS within perturbation theory.This began with the pioneering work in the linear approximation [12][13][14][15] and then extended to higher orders in perturbation theory [16][17][18][19][20][21].In the meantime, progress has been made in numerical simulations, including full relativistic simulations [22], as well as the adoption of appropriate gauges to reinterpret Newtonian simulations [23,24].
When interpreting observations, it is necessary to have a robust theoretical framework to describe the statistical properties of the observed structures.This framework must include not only their distribution across the sky (i.e. the geometric mean), but also their intrinsic stochastic origin, which ultimately shapes the late-time structures.In the last few decades, an active line of research has emerged, aiming at establishing a well-posed mathematical prescription for the averaging of scalar quantities concerning both space-like hypersurfaces [25][26][27] and, more recently, generalized also to light-like ones [28][29][30][31][32][33].
Non-Gaussian distributions are characterised by non-vanishing higher moments.In our work, we introduce the formalism for evaluating the skewness of any cosmological observable on the light cone at the leading order in perturbation theory.In applying this formalism to the luminosity distance we are motivated by the numerical results presented in Ref. [62], which highlighted significant deviations from Gaussianity in the luminosity distance PDF.We will show that the leading order terms solely involves second-order perturbations of the luminosity distance.However, when dealing with 1-point functions, a direct comparison is not possible due to the limitations within the validity of perturbation theory.To remove the non-linear effects, small scales must be smoothed out in the analysis.Another key question is whether these non-Gaussian features are a direct evidence of the expected late-time non-Gaussianities or they are somehow sourced by spurious effects in the simulations, such as the finite sampling, limited sky-coverage or redshift binning.
Despite these limitations, our results are in line with those obtained from the numerical light tracing in Ref. [62], especially as we approach higher redshifts, where the effect of nonlinearities diminishes at fixed scale.While a more direct comparison, where the numerical and perturbative approaches are smoothed at the same physical scale, remains a prospect for future work, our study confirms that the dominant contribution to the skewness comes from the matter bispectrum, which affects the lensing propagation.Meanwhile, other relativistic effects play a sub-dominant role in the observed skewness.In summary, in this paper, motivated by the results obtained in [62], we are interested in characterizing the distribution of the luminosity distance in order to quantify non-Gaussianities.Our work aims to develop a general method to evaluate analytically non-Gaussianities for any desired cosmological observables.Then, we will apply our general results to the Hubble-Lemaître diagram.In this regard, bearing in mind the general average prescription in cosmology over the past light-cone [32], we will focus on the skewness of the distribution of a generic scalar observable, like the luminosity distance-redshift relation, averaged over a region of space-time.
This paper is organised as follows.In Sect.2, we review the general features of the evaluation of the average and dispersion at the leading order in perturbation theory.In Sect.3, we apply the same formalism for the evaluation of the average and dispersion to the computation of the skewness at the leading order.Sect. 4 is devoted to the explicit evaluation of the skewness for the distance-redshift relation in the perturbed FLRW Universe.Here, we limit our computation to the lensing terms.Afterwards, we specialise our science case to the number count weighted measure in the infinitesimal redshift bin limit.In Sect.5, we explicitly evaluate the leading terms to the skewness and show how they are connected to the matter power spectrum and to the matter bispectrum.Numerical results for a fiducial ΛCDM model are presented in Sect.6, whereas Sect.7 is devoted to the (possible) comparison with numerical simulations.Finally, our main conclusions are outlined in Sect.8. Appendices A and B report explicit technical evaluations needed to derive our results.
In this work we adopt the following 3-D Fourier transform convention 2 Leading-order terms review: average and dispersion In this section, we will discuss how to evaluate the moments of a given observable.The specific choice for the measure adopted for their evaluation will be left for later.In particular, we will adopt the covariant prescriptions for different kinds of light-cone averages presented in [32].Let us start then by considering a generic observable S. Its average over a given portion of spacetime is given, in complete generality, by where µ represents the measure weighting the average procedure.Eq. (2.1) is exact but, for our purposes, can be expanded perturbatively.We work in a perturbed FLRW metric for scalar perturbations where1 Φ ≡ i ϕ (i) and Ψ ≡ i ψ (i) and i refers to the order of the perturbative expansion.So by expanding to second order the observable S and the measure µ, we have S ≃ S (0) 1 + σ (1) + σ (2) , dµ ≃ dµ (0) 1 + µ (1) + µ (2) . (2. 3) The fundamental requirement to apply Eq. (2.1) is that S transforms as a scalar field under a generic coordinate transformation.With this proviso, Eq. (2.1) can be applied to the generic α-th power of S as well.
According to the perturbative scheme outlined in Eqs.(2.3), we can then rescale our observable to its background value, namely apply Eq. (2.1) to S/S (0) rather than S. On one hand, this choice simplifies a lot the equations and is in line with previous estimations provided in literature (see [34]).On the other hand, and more importantly, this is the result that has been provided by numerical simulation (see, for instance, the histogram in Fig. 2 of [62]).We then have (1) ] +α I[µ (1) σ (1) ] + α I[σ (2) ] + α (α − 1) 2 I σ (1) 2 − α I[µ (1) ] I[σ (1) ] , where we have defined The average meant by ⟨. . .⟩ is a geometrical procedure which takes into the smoothing of an arbitrary inhomogeneous manifold but does not specify the nature of these inhomogeneities.This nature is given by a further ensemble average a . . .a which provides the stochastic properties of the inhomogeneities.Having in mind that our perturbations will be sourced by the fluctuations of the gravitational potentials Φ and Ψ, we just assume for now that for the linear gravitational potentials ϕ = ψ = 0 and that higher-order powers of Φ and Ψ can have a non-vanishining ensemble average.We will provide later the mathematical properties for these operations.We then get 1) ] I[σ (1) ] ,(2.6) where we used the abovementioned properties that the ensemble average of linear perturbations vanishes.Eq. (2.6) can be readily applied to the evaluation of the average of S/S (0) : for α = 1 indeed we have m ≡ S S (0) = 1 + I[µ (1) σ (1) ] + I[σ (2) ] − I[µ (1) ] I[σ (1) ] . (2.7) This result already provides an important intermediate step to evaluate the generic α-th moment µ α of the distribution of S. In fact, this can be written as 1) σ (1) ] +α I[σ (2) ] − α I[µ (1) ] I[σ (1) ] + Eq. (2.8) is useful also for the evaluation of the variance, i.e. σ 2 ≡ µ 2 , where we have (2.9) However, all the higher moments µ α with α > 2 vanish.This is due to the fact that in this section we have limited our non-linearities to be at most of second-order in perturbation theory.In the next section, we will show how the next-to-leading order terms will impact the analytic estimation of higher-order multipoles, such as µ 3 .
3 Next-to-leading order term: skewness In this section, we will provide the analytic evaluation for the skewness of a generic observable S. As said in the previous section, we recall that we remove the background contribution for the observable by studying S/S (0) rather than just S. In order to proceed to our end, we first need to introduce, in complete generality, the standardized moments κ α as However, as we have pointed out in the previous section, second-order corrections are not enough to evaluate higher-order moments.Naively, one would expect that each µ α demands α-th order perturbations to get the leading term.This estimation is accurate when α label even moments.For what concerns the odd moments, however, leading terms are of α + 1-th order in perturbation theory (as an instance, this occurs for the average and the dispersion).This peculiar hierarchy is due to the fact the leading term for each moment µ α is I[σ (1) α ].Therefore, since σ (1) ∼ ψ (1) , we have that the α-th order term of µ α is related to the α-point correlation function of ψ (1) : as a consequence of Wick theorem, by assuming Gaussian initial conditions, the latter is non-null only when α is even.
According to what discussed so far, the evaluation of higher order moments would require the estimation of perturbations for both the measure µ and the observable S beyond the second order at least for the leading terms.These quantities are not available in literature and evaluating them is a non-trivial task which would be an interesting result per se.Fortunately, it can be shown that the expansion in Eqs.(2.3) still catches all terms for what concerns third moment at leading order even when perturbations of S and dµ up to the fourth order are consistently taken into account.The derivation is quite long but straightforward and the interested reader can refer to Appendix A for this proof.Hence, we get µ 3 = I[σ (1) 3 ] + I[σ (1) 3 µ (1) ] − I[σ (1) 3 ]I[µ (1) ] + 3I[σ (1) 2 σ (2) ] +3 I[σ (1) 2 ] I[µ (1) ] I[σ (1) ] − 3I[σ (1) 2 ] I[µ (1) σ (1) ] − 3I[σ (1) 2 ] I[σ (2) ] = I[σ (1)3 µ (1) ] + 3I[σ (1) 2 σ (2) ] − I[σ (1) 3 ]I[µ (1) ] − 3 σ 2 (m − 1) . (3.2) where m and σ 2 are respectively given by Eqs.(2.7) and (2.9).We remark again the disappearing of first term in the first line of Eq. (3.2): this would be the only third-order contribution but vanishes as a consequence of Wick theorem, as anticipated before.This result is completely general in regard of the kind of observable S and the chosen prescription of the spatial average µ.In the next section, we will apply them to the specific cases of the distance-redshift relations.
4 Higher-order moments for the distance-redshift relations In this section, we want to apply the general formalism previously developed to the explicit case of the luminosity distance-redshift relation d L (z).As recently pointed out in Ref. [62], the probability distribution of the luminosity distance d L (z) exhibits a significant non-Gaussian contribution.Our objective is to investigate the extent to which this behavior can be captured within perturbation theory.Linear and non-linear perturbations for the distance-redshift relation have been discussed in several papers [35,37,[63][64][65][66].For the purposes of this work, we refer to the leading lensing terms as given in Eqs.(B.1) and (B.2) of [40].Moreover, we will assume that there is no anisotropic stress at linear order, i.e. ϕ = ψ.Then, we can express the luminosity distance to second order in perturbation theory as follows with and where d L (z) is the luminosity distance in a FLRW background, r s is the comoving distance to the source, ∆ 2 is the dimensionless angular Laplacian on the 2-sphere, and ψ(r) = ψ(η o − r, r).We have also introduced where γab 0 = diag(1, sin −2 θ), and γ ab 0 = r −2 γab 0 .In order to evaluate the skewness of the PDF, we also need to choose a measure µ, according to Eq. (3.2).Several possibilities have been discussed in [32] for measures, which are general covariant and gauge invariant.In this work, we have chosen to utilize the galaxy number count weighted measure for computing the averages.This measure has been discussed for the backreaction of stochastic inhomogeneities to the mean value of the Hubble diagram in the limit case of small width of the redshift bin [31,32,41].However, it can be applied directly also for the case of finite redshift bin, preserving the property of general covariance [32].
In regard to the comparison with numerical simulations, the adoption of the average weighted with galaxy number count appears to be the most appropriate choice.This weighting scheme naturally assigns more significance to regions with higher density contrasts, allowing for a more accurate and meaningful comparison.The background term in the measure is then given by where dΩ is the infinitesimal solid angle, d A is the background angular distance, ρ (0) is the background density and the denominator comes from the background expansion of the number counts.We are then left with reporting the linear order galaxy number counts [12][13][14]68].The leading terms in the weak field expansion are given by [18] µ (1) where δ is the matter linear perturbation (with the galaxy bias set to b 1 = 1) and H −1 ∂ r v ∥ accounts for the redshift space distortion.With Eqs.(4.4) and (4.5), we can compute µ 3 .In details, from Eq. (3.2), we explicitly get where leading order term for σ 2 is given by Eq. (2.9), and the average is given by In the derivation of Eq. (4.6) we made use of I[σ (1) ] = 0.This is because we have considered only lensing terms in Eqs.(4.2) and terms sourced by a Laplacian vanish when averaged over spatial directions.The same result holds also in the derivation of Eq. (4.7) in regard of the term I σ LSS .In a similar manner, also terms as ).This is because the ensemble average for these terms always correlates the monopole of δ or H −1 ∂ r v ∥ with the Laplacian of σ (1) , selecting then its null eigenvalue.Sub-leading relativistic corrections to the distance-redshift relation, such as Doppler effect, can contribute with a non-vanishing monopole, especially if sources at small redshifts are taken into account.However, for the purposes of this pioneering work, we have not considered them, even thought they can be of interest for the analysis of close sources.We will investigate their presence in a forthcoming work.

Infinitesimal redshift bin
Before concluding this section, some comments about the integration domain are in order.
The integration over the solid angle dΩ in Eq. (4.4) covers the entire solid angle.However, for the integration over the redshift we adopt a redshift bin of width ∆z such that all the needed averages are evaluated within a certain redshift range [z s , z s + ∆z].In case of large value of the redshift bin width, all the redshift dependent terms in Eq. (4.4) needs to be integrated.However, in the limit when ∆z becomes small enough to be considered infinitesimal4 , i.e. ∆z → δz, the measure in Eq. (4.4) simplifies to and then the integrals in the definition of I[f ] in Eq. (2.5) reduce to simple angular integrals.Moreover, within this limit, all the redshift dependent quantities and δz itself factorize out of the integrals and cancel out in the ratio with the unconnected diagram in Eq. (2.5), returning then In this limit, the non-purely second-order terms in the first equality of Eq. (3.2) cancel.To this end, we anticipate in the limit of Eq. (4.9) that where we have used Eq.(2.9) and the fact that the ensemble average can not have any preferred direction due to statistical isotropy.This term exactly cancels the counterpart coming from Eq. (2.7) in Eq. (4.6).Here, we make the first interesting remark that the weight for the measure is irrelevant in the infinitesimal bin limit for the leading order terms of the skewness for the luminosity distance-redshift case study.Moreover, in this limit, the third moment reduces to the simpler expression where we have defined LSS . (4.12) We have combined Eqs.(4.6) and (4.7), adapted for an infinitesimal redshift bin by using Eq.(4.10), to cancel out terms like I[σ (1) 3 δ] and I σ (1) 3 ∂rv ∥ H , as previously mentioned.
Eq. (4.11) is entirely sourced by pure non-linear terms in the expression of the distanceredshift relation.Indeed, when considering Gaussian initial conditions, linearly evolved perturbations preserve their probability distribution.Consequently, they do not generate any non-Gaussianity in the limit of infinitesimally narrow bins, where the hypersurfaces for the averages are evaluated at constant redshift.
Hence, the terms in Eq. (4.11) are sourced by pure non-Gaussian effects and their labels follow accordingly: in Eqs.(4.12) Q, P B ans LSS respectively stand for Quadratic, Post-Born and Large-Scale-Structure since • µ Q 3 takes into account the non-Gaussianity coming from the quadratic term σ (1) 2 /2 in Eq. (4.2), • µ P B 3 contains all the relevant terms due to the Post-Born corrections to the d L (z), such as multi-lens effects, • µ LSS 3 catches the non-Gaussianities arising from the bispectrum of δ (2) δ (1) δ (1) .
Understanding how much the finite bin effect may mimic non-Gaussian behavior beyond the fundamental non-linearities is of interest per se.We will present this study in a forthcoming work.For the rest of this work, we will explicitly compute and numerically evaluate the amplitude of the expected skewness in the infinitesimal redshift bin case, where the only relevant effect is the one due to intrinsic non-Gaussianities in the inhomogeneities.

Analytic expressions
In this section, we provide the explicit expression for the leading order terms of the skewness in the infinitesimal redshift bin limit.First let us remark that our derivation so far is completely geometrical without any assumption on the underlying theory of gravity, with the only assumption that light propagates along null geodesics.At this point we need to use General Relativity to relate metric and matter perturbations.We first provide some general preliminaries needed to follow our derivations.Technical details are reported in the Appendix B. In order to compute the different contribution to the third moment defined in Eqs.(4.12), it is convenient to introduce the generalized Hankel transform of the matter power spectrum at two given times η 1 and η 2 , namely where is the ℓ-th order spherical Bessel function and D 1 the growth function normalized to unity today.The reason why we decide to adopt the definition (5.1) for the Hankel transform is that it can be easily generalised to the science case of non-linear matter power spectrum, as we will discuss in Sect.6.
In evaluating Eqs.(4.12), we encounter the following terms We remark that whereas second last of Eqs.(5.2) vanishes since can not depend on the direction n due to statistical isotropy.For the same reason Last equality in Eq. ( 5.2) can be understood in the following way: they always involve an odd number of angular derivatives and hence they vanish due to statistical isotropy, since they naturally introduce a preferred direction, which returns to 0 when the ensemble average acts.The 2-point correlation functions in Eqs.(5.2) are enough also to evaluate the 4-point correlation functions of interest for us.With the aim of Wick theorem, indeed, we can write (5.6) For the sake of completeness, we report in Appendix B the derivation of Eq. (5.6).In a similar manner, we obtain also that (5.7) It is interesting to notice that only two permutations survive in Eq. (5.7) as a consequence of the last two equalities in Eqs.(5.2).Looking at the structure of Eq. (5.7), we notice that only the permutations mixing ∆ 2 ψ with one of the term of ∆ 2 (∂ψ) 2 survive.Finally, last 4-point correlation function of our interest is The interesting thing about Eq. (5.8) is that only one permutation survives in the final result and this is again a consequence of last of Eqs.(5.2).The structure of Eq. (5.8) also exhibits a complete factorization of the 2-point correlation function in the form (∆ 2 ψ) 2 .This situation is opposite to what we have shown for Eq.(5.7), where only mixed terms survive in the final permutations.The impact of these differences in the ultimate evaluation of µ 3 will be to significantly reduce the final number of non-null terms for the skewness, as we will show in Sect.5.2.These analytic preliminaries are enough to provide the explicit expressions for µ Q 3 , µ P B 3 in the infinitesimal bin case.The evaluation of µ LSS 3 is more delicate and will be treated in a specific way.

µ Q
3 : quadratic terms The starting point for the computation of µ Q 3 is its expression in Eqs.(4.12).To our hand, we first evaluate the variance in the small bin limit.By making use of Eqs.(2.9), (4.2), and (5.2), we get that (5.9) In the same way, thanks to Eq. (5.6), we have that (5.10) Hence, the combination of Eqs.(4.2), (5.9) and (5.10) we get the quite simple result It is worth noticing what happens for the leading term of the standardized third moment.For the quadratic terms coming from Eq. (5.11), we get that namely the skewness of the distance-redshift relation due to the quadratic corrections is proportional to the dispersion of d L (z).This result already allows us to estimate the amplitude of κ Q 3 .Indeed, in [36] the dispersion for the distance-redshift relation has been estimated for the lensing contribution at higher redshift to be ∼ 1%.This evaluation takes into account the non-linear power spectrum for the gravitational potential as provided by the Halofit model [69].

µ P B
3 : Post-Born corrections For the evaluation of the Post-Born terms µ P B 3 in Eq. (4.12), we first look at the expression of Σ (2) in Eq. (4.3) and treat its terms separately.We divide mix , where we defined such the total contribution to µ P B 3 is given by and (5.16) Hence the total contribution of the Post-Born corrections to the third moment is and the skewness is (5.18) The structure of µ P B 3 reveals a sequence of five nested line-of-sight integrals, consistent with our expectations from Post-Born corrections.These corrections effectively incorporate the non-linearities arising from multi-lens effects.To address numerical challenges and minimize the number of integrals, we propose an approximation based on the observation that L(r 1 , r 2 ) exhibits a strong peak near r 1 ≈ r 2 .Thus, we can approximate L(r 1 , r 3 ) as ∼ δ D (r 1 − r 3 ), simplifying the calculations.We then have where f (r 1 ) is a function to be determined in a numerical way and Θ(x) is the Heaviside step function.This means that For what concerns f (r 1 ), since the latter is not a function of r and the integrand L contributes to the integral in Eq. (5.19)only around r 3 ≈ r 1 , we have that integral (5.19) is independent of the value of r as long as r > r 1 .To this end, then, we can evaluate the function f once for all by choosing an r equal or larger than a given comoving distance r * well-beyond the highest redshift that we investigate.For our purposes, we have chosen r * = r(z = 4) but the final results are clearly independent of this choice.In this way, we have that Numerical results within this approximation will be discussed in Sect.6.

3
: the role of the bispectrum The evaluation of µ LSS 3 is the more demanding, since involves the 3-point function where we have defined Ψ (2) ≡ 1 2 ψ (2) + ϕ (2) .To evaluate this correlation, we use again the Poisson equation to relate Ψ (2) to the second-order matter fluctuations δ (2) .We then have that where 5 C ≡ 27 8 a(r 1 )a(r 2 )a(r 3 ) , and we have introduced the (non-symmetrized) bispectrum which at tree-level is given by with (5.26) Finally, accounting for the two other permutations of Eq. ( 5.22) we need to simply replace B(k 1 , k 2 , k 3 ) with its symmetrized version (5.27) From Eq. (5.23), we can write the 3-point function as where denotes the so-called 3-j Wigner symbol, which is non-null only when m 1 + m 2 + m 3 = 0 and (ℓ 1 , ℓ 2 , ℓ 3 ) satisfy the triangular inequality.The presence of these symbols selects only specific shapes for the sums in ℓ-space.The detailed derivation of Eq. (5.28) can be found in Appendix B.Here we directly report the result which are useful for us to discuss the underlying physics.Despite its compact form, Eq. (5.28) is still quite demanding for our scopes.Indeed, we still need to evaluate the contribution to the third moment as given in Eq. (4.12).This implies that three line-of-sight integrals must be performed over Eq. (5.28).In total, then, we have to numerically compute seven integrals for each non-null set of (ℓ 1 , ℓ 2 , ℓ 3 ) and this is quite unpractical from the numerical viewpoint.
In order to face the dramatic need for reducing the number of integrals, we invoke the Limber approximation for each integral in k-space in Eq. (5.28).Thanks to it, in fact, we can make the following approximation [70-72] whenever the function g does not vary too much rapidly.In this way, the approximation (5.30) can be used three times to get rid of the k-space integrals.This reduces the number of remaining integrations to four.Hence, Eq. (5.28) becomes where k i ≡ ℓ i +1/2 r i .The three delta Dirac distributions in Eq. (5.31) simplify the three lineof-sight integrals in the last relation of Eq. (4.12), leading to the final expression for the LSS contribution to the third moment where we used Eqs.(4.2) and (4.3), and we defined now k i = ℓ i +1/2 x and η x ≡ η 0 − x.Then, we have that the skewness due to the non-linear structures is (5.33)

Numerical results
In this section, we discuss the numerical integrations of the analytical results for the skewness given by Eqs.(5.12), (5.18) and (5.33).We work with the fiducial cosmology given by the ΛCDM parameters Ω c = 0.2638, Ω b = 0.04827 h = 0.67556, A s = 2.215 × 10 −9 and n s = 0.9619 in accordance with [62].

The smoothing scale
The evaluation of the skewness, or any other 1-point function, is strongly sensitive to the nonlinear nature of structure formation.However, by working in perturbation theory we know that our prediction will fail beyond the mildly non-linear scale.Therefore, our prediction can be compared to real observations only when the non-linear effects are filtered out, by introducing a smearing scale in real space of size ρ.In this regard, we use a spherical top-hat window function W(r, ρ) in real space.This procedure introduces a window function W (k, ρ) in Fourier space given by which affects the power spectrum as With such smoothing scale, all the integrals over the momentum k converge quickly and the choice of k max become numerically irrelevant.
For the term κ LSS 3 we will consider also some non-linear fitting expression of the bispectrum, in particular BiHalofit [73].For this reason, to smooth out small scales non-linearities we apply directly the window functions as follows To include the window function we simply need to replace eq. ( 5.32) with In the following, we will investigate how suitable choices for ρ can be made: we will consider the cases of ρ equal to 5, 10 and 20 Mpc/h.

Quadratic and Post-Born terms
We start our numerical investigation by looking at the contributions to the skewness arising from the linear gravitational potential, derived in Sects.5.1 and 5.2.For this discussion, we will use the cutoff in real space on small scales as prescribed by Eq. ( 6.2).Results are shown in Fig. 1 respectively for the choices of ρ = 5, 10, 20 Mpc/h.For each case, we consider linear matter power spectrum (left panels) and Halofit model for the matter power spectrum (right panels).
With these results, we can spot common features among the different cases.First of all, we notice that quadratic terms in κ Q 3 (blue solid lines in Fig. 1) are always positive and lead to an increasing skewness with redshift.The order of magnitude of this term reaches about 10 −2 .In this regard, the positive sign was already expected just by looking at the structure of κ Q 3 in Eq. (5.12), since it is just proportional to the dispersion of the d L (z) distribution.We have also checked that our framework is in a good agreement with [36] when we push the coarse-graining scale up to 0.3 Mpc/h, roughly corresponding to the UV cutoff k U V = 10 h/Mpc used in [36].
For what concerns the Post-Born contribution to the skewness, namely κ P B 3 in Eq. (5.18), blue dashed lines in Fig. 1 show two important features.First of all, the Post-Born contribution to the skewness is always negative.This is important to be addressed since κ P B 3 in Eq. (5.18) does not show any manifest sign, contrary to κ Q 3 .As a second matter of fact, Post-Born corrections contribute to the skewness with the same order of magnitude of κ Q 3 .This is somehow in line with the structure of κ Q 3 and κ P B 3 , since they just involve different the line-of-sight integrals of the same kernels L(r, r ′ ).Hence, these two features lead to a competitive effect between κ Q 3 and κ P B 3 .Indeed, we have a neat effect for the linear gravitational potential skewness, which is still positive but attenuated (solid red lines in Fig. 1).We remark that all the features are shared regardless the value of ρ and whether linear or non-linear matter power spectrum is considered.
The actual contribution of the non-linear scales to the sum κ Q 3 + κ P B 3 is quantified in Fig. 2. Indeed, here we show the relative difference as for different smoothing scales ρ.For a value of ρ = 20 Mpc/h, the non-linearities in P ρ (k) marginally contribute to the total effect with a relative correction of ∼ 0.1%.The non-linear scales happen to be more relevant when ρ decreases, namely with a few percent and almost 10% relative correction respectively with a coarse-graining radius of 10 Mpc/h and ρ = 5 Mpc/h.The impact of non-linearities in the matter power spectrum also decreases by going to higher redshifts.Let us remark that, in principle, replacing the linear power spectrum in the perturbative expansion with Halofit is not self-consistent.However this can give us a rough estimation of the expected accuracy of our prediction based on perturbation theory for different smoothing scales.

Bispectrum
The numerical evaluation of the κ LSS 3 from Eq. (5.33) requires some subtleties to be accounted for.In fact, from Eq. ( 5.32) we have to specify the value of ℓ max to effectively compute the contribution of the bispectrum to the skewness.Since the physical cutoff is determined by the smoothing scale ρ, we only need to ensure that ℓ max ≳ r s /ρ, where the specific value is chosen to guarantee a percent precision.
We perform numerical investigations by using linear matter power spectrum, Halofit model and the BiHalofit fitting function for the bispectrum provided in [73].For different values of the coarse-graining scale ρ, the numerical results are reported in Fig. 3, where we have adopted linear power spectrum in the left panel, Halofit model in the center and BiHalofit in the right panel.At first sight, a general feature is that the absolute value of κ LSS 3 is decreasing with the redshift.In particular, whereas the values at z ∼ 1 are of order 0.1, the skewness at lower redshifts varies more and becomes O(1) when ρ becomes smaller.Moreover, the value is quite insensitive of the kind of spectrum adopted for ρ = 20 Mpc/h.This last feature is in line with the fact that the involved scales still evolve almost linearly for this case.
A similar behavior occurs also for the case ρ = 10 Mpc/h.Even in this case, the value of ρ is still such that non-linear features in the power spectrum are marginally relevant, and then the Halofit and BiHalofit models return results that are almost alongside the ones obtained by the linear power spectrum.However, we notice that the overall amplitude increases quite a lot when we lower the coarse-graining scale from 20 to 10 Mpc/h.As a quantitative instance, the comparison between red and blue curves in Fig. 3 at z = 0.1 exhibits an increasing in the absolute value of the skewness that is of ∼ 40%.This confirms that a significant amount of information is encoded in the skewness within the scales of 10 to 20 Mpc/h.However, as the prediction from linear theory does not deviate significantly from the estimations using both Halofit and BiHalofit, we are still within a regime where perturbation theory has not yet failed.
Finally, with orange curves in Fig. 3 we report the numerical results for ρ = 5 Mpc/h.For this case study, we have new emerging features.First of all, we start to appreciate a more prominent difference between the linear power regime and the non-linear one, and this difference is more evident at smaller redshifts.In fact, here we have an enhancement when the BiHalofit model is considered of ∼ 10%, as can be appreciated also in Fig. 4, where we show the analogous of Fig. 2 and compare the Halofit fitting formula against the BiHalofit one, but only for κ LSS

3
. This increasing of κ LSS 3 , as given by Eq. (5.33), follows the separate increasing of the two quantities in its ratio, namely µ LSS 3 in Eq. (5.32) and σ 3 as derived from Eq. (5.9).Given that, we have that the numerator and the denominator of this ratio separately increase when we decrease ρ from 20 to 5 Mpc/h.Hence, by considering the discrepancy between the linear power spectrum and BiHalofit as an indicator of the validity range of perturbation theory, we observe that for ρ = 5 Mpc/h, the accuracy of perturbation theory cannot be trusted beyond approximately 10%.We also remark that at those redshifts our results are only indicative of the order of magnitude, since other relativistic effects here neglected, such as the Doppler correction might significantly alter the result.This is indeed the case for the dispersion, as shown in Fig. 8 of [36].We plan to investigate more in details the contribution of other relativistic effects on small redshift in a forthcoming paper.On the other hand, we also point out that the result is quite stable for z > 0.4 when we decrease the value of ρ.
Two final comments are in order before concluding the discussion devoted to κ LSS

3
. First of all, we recall that our derivation is based on the Limber approximation for the lensing terms.This gave us an appreciable speed-up of the numerical evaluation, since it reduces the number of numerical integrations to one line-of-sight integral, as given in Eq. (5.32).On the other hand, for smaller redshift, the accuracy of this result might be somehow limited, especially because the discrete sum is limited to larger angular scales when the Limber approximation shows some issues.The second comment is about the use of the Halofit model to account for the non-linearities in the matter power spectrum.This is clearly a limitation, again on smaller redshifts, when the role of non-linearities in the LSS is more prominent and this is also made manifest by the fact that BiHalofit enhances by almost 10% the amplitude of obtained with Halofit at ρ = 5 Mpc/h.Both these shortcomings are less serious at higher scales since i) we weight more the smaller angular scales, where the Limber approximation is more accurate, and ii) we move the peak of the lensing kernel to higher redshifts, where the impact of non-linearities is still relevant but attenuated.
Besides all the numerical limitations due to the usage of approximations, a take-home message that emerges from this study is that κ LSS 3 is always dominant with respect to κ Q 3 and κ P B 3 , regardless of the choice of the smoothing scale ρ and the use of linear or non-linear matter power spectrum, as one can see from a direct comparison of Figs. 1 and 3. Hence, the overall sign of the skewness in the investigated range of redshift due to the lensing is always negative, in accordance with the indication of the numerical simulation of [62].However, a deeper discussion about how to compare our results with [62] is in order.This is the topic of the next section.

Comparison with numerical simulations
In this section, we will address the challenges that arise when attempting to compare our results with those obtained from ray-tracing across the numerical simulation presented in [62].We will also outline the recommended approach to ensure a valid and meaningful comparison.To begin, we provide a brief overview of the key elements employed in [62].The ray-tracing of light-like geodesic in [62] has been performed through a relativistic Universe simulated with gevolution [22], with the following cosmological parameters: the reduced Hubble constant h = 0.67556, the cold dark matter density Ω c = 0.2638, and the baryon density Ω b = 0.048275.Furthermore, the radiation density that includes massless neutrinos with N eff = 3.046 and linear initial conditions are computed with CLASS [74][75][76][77] at redshift z ini = 127, assuming a primordial power spectrum with amplitude A s = 2.215×10 −9 (at the pivot scale 0.05 Mpc −1 ) and spectral index n s = 0.9619.In this regards, the skewness has been obtained with the following remarks: 1. all the ray-traced light-like geodesics start from a point where a local overdensity has created an Halo; 2. structures evolve in a box with size 2.4 Gpc/h whose grid space is 312.5 kpc/h.Hence, we expect that structures can be investigated roughly up to a scale in Fourier space of 3 h/Mpc; 3. the distribution of the obtained distance-redshift relation is binned in four bins of width 0.5, ranging from z = 0 until z = 2.
For the sake of completeness, we also report that the ray-tracing in [62] is performed through a light-cone with a partial covering of the observed sky of 450 deg 2 .Moreover, the simulations in [62] link the density contrast to the gravitational potential beyond the Newtonian approximation through the second-order Hamiltonian constraint (see Eq. (2.19) of [78]).With this in mind, we have then extrapolated from Fig. 2 in [62] the values for the skewness in the four bins shown in Table 1.The results of Table 1  with what we have obtained in Sect.6 for the first two bins, whereas they match better at higher redshifts.However, we have to keep in mind the previous bullet points to understand how our results should (or better, could) be interpreted against [62].First of all, bullet point 1 states that only regions in the simulated Universe where the overdensity is high enough to have created Halos are spanned by the ray-traced distanceredshift distribution.This motivates well the adoption of a number-count weighted prescriptions for the averages in Eq. (4.5) since this naturally weights more regions where structures are more likely to have been created.
Moving forward with the discussion, bullet point 2 tells that no smoothing-scale procedure has been introduced in [62], beyond the grid space of the numerical simulation.Moreover, this implies that the expected results can probe scales up to the deeply non-linear regime and this is somehow problematic for the analytic investigations, since the validity of our perturbative scheme is at least questionable on those scales.However, this point could be easily overcame, since a smoothing scale procedure could be applied as well to the simulated Universe.
Finally, for what concerns bullet point 3, the bin width adopted in [62] is quite large when compared to our infinitesimal bin formalism and this can introduce spurious contamination to the extrapolated skewness in two regards: first, by introducing background effects which are not at all related to the LSS and, secondly, by introducing a plethora of other relativistic effects, such as the contamination due to cross correlations with density fluctuations and redshift-space-distortion in Eq. (4.5), which could complicate the actual analysis.This point can be treated either by slicing the redshift space with more narrow bins or by developing our analytic formalism to the finite bin case.We plan to achieve the latter task in a forthcoming work.
Given all the assumptions and the intrinsic differences that are present between our analytical findings and the numerical ones of [62], we believe that it is quite remarkable that the two results share the same signature and differ only by an order unity factor.In our opinion, our results open an interesting window to resolve the previous bullet points and provide an ultimate comparison.

Summary and Conclusions
In this work we have provided for the first time an analytic evaluation of the skewness for the distance-redshift distribution in ΛCDM.The primary goal has been to provide an analytic scheme to compare and understand the analogous results obtained from numerical simulations [62].To this end, we have applied the covariant and gauge-invariant formalism for the lightcone averages detailed in [32] to evaluate the higher-order moments of the distribution.
Within our framework, we have found that the second-order perturbations of the distanceredshift relation are enough to evaluate the leading-order terms of the skewness, whereas only linear-order corrections of the measure are needed.In this regard, among all the formally viable prescriptions for these averages, we have adopted the number-counts weighted one for our evaluation.This is of particular importance since it is ideally in line with the way the ray-tracing of photons across the above-mentioned simulation is done.
Furthermore, we have focused our treatment to the infinitesimal bin width for the redshift distribution.In this case, we have outlined several cancellations in the leading-order terms, see Sect.4.1.According to these, we are left only with three terms in the final formula (4.12) due to lensing, Post-Born corrections and higher-order gravitational potential.In particular, the latter sources the total skewness with the matter bispectrum, integrated along the line-ofsight.Given that, what has emerged from numerical investigation is that this also provides the biggest contribution among all.On the other hand, lensing and Post-Born terms exhibit a competitive effect which turns out to be positive and small when compared against the bispectrum one.
For what concerns a closer comparison with numerical simulations in [62], it turned out that our analytical estimation of the skewness shares the same (negative) amplitude.This is quite remarkable, given the number of approximation that we had to require to obtain our results (see Sect. 7 for a detailed discussion).From a physical viewpoint, our skewness from the bispectrum exhibited a strong dependence on small scales and this is not surprising, since it is a one-point function.What is important to remark is that, in order to get a result independent of the UV behavior, we have introduced a coarse-graining (or smoothing) scale in real space and quantified the impact of this smoothing scale on the final amplitude.In this regard, we have shown that varying this coarse-graining scale from 20 to 10 Mpc/h contributes to a 40% increasing of the total amplitude at small redshift (z ∼ 0.1).This dependence is still important even though less severe for distant sources (z ∼ 1).This seems to suggest that there is room for better agreement with the results of [62], since they have probed smaller scales than ours (see again Sect.7).
To conclude, this proves that we need only corrections up to second-order for the observable and linear in the measure perturbation to compute the high-order moments.We remark that this is not ensured for the moments of S, since the quantity S (0) − m (0) does not vanish on the background, and then the generic moment of the distribution of S includes high-order corrections for the observable and the measure.

(B.5)
It is now worth to focus on the sum over ℓ 1 appearing into Eq.(B.5).To this end, we recall that the Legendre polynomials P ℓ (cos θ) are eigenfunctions of the 2-D Laplacian with eigenvalues −ℓ(ℓ + 1) and that P ℓ (1) = 1.Hence, we have as prescribed by the general dictionary provided in [79] for the relations among the different cosmological transfer functions.

Figure 1 .
Figure 1.Quadratic (blue), Post-Born (blue dashed) terms and their sum (red) for the skewness, evaluated at with a coarse-graining scale ρ = 5 Mpc/h (top), ρ = 10 Mpc/h (center) and ρ = 20 Mpc/h (bottom).Left panels are obtained with linear power spectrum, whereas right panels consider Halofit model for the matter power spectrum.

Table 1 .
[62]es of the skewness for the distance-redshift relation distribution obtained in[62]from the ray-tracing of photons across a ΛCDM Universe simulated with gevolution.Results have been obtained by binning the dataset in four redshift bins from 0 to 2 with bin width of 0.5.