Cosmic variance of H_0 in light of forthcoming high-redshift surveys
CCosmic variance of H in light of forthcoming high-redshift surveys Giuseppe Fanizza , ∗ Bartolomeo Fiorini , † and Giovanni Marozzi , ‡ Instituto de Astrofis´ıca e Ciˆencias do Espa¸co, Faculdade de Ciˆencias da Universidade de Lisboa,Edificio C8, Campo Grande, P-1740-016, Lisbon, Portugal Institute of Cosmology & Gravitation, University of Portsmouth,Dennis Sciama Building, Burnaby Road, Portsmouth, PO1 3FX United Kingdom Dipartimento di Fisica, Universit`a di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, Pisa, Italy (Dated: February 25, 2021)Forthcoming surveys will extend the understanding of cosmological large scale structures up tounprecedented redshift. According to this perspective, we present a fully relativistic framework toevaluate the impact of stochastic inhomogeneities on the determination of the Hubble constant. Tothis aim, we work within linear perturbation theory and relate the fluctuations of the luminositydistance-redshift relation, in the Cosmic Concordance model, to the intrinsic uncertainty associatedto the measurement of H from high-redshift surveys (0 . ≤ z ≤ . H which will be more precise than the one obtained from local sources, at least in regard of theintrinsic uncertainty related to a stochastic distribution of inhomogeneities. INTRODUCTION
The last 30 years have witnessed the evolution of cos-mology from an order of magnitude description of theUniverse into an area of the hard sciences where preci-sion measurements are achievable. This change has itscornerstone in the estimation of cosmological parame-ters through the detection of Cosmic Microwave Back-ground (CMB) spectra which provide so far the mostprecise measurements of these parameters. Given thisastonishing success, forthcoming missions aim to pushforward this result in order to infer measurements of cos-mological parameters also from late-time dataset, suchas Large Scale Structure (LSS) surveys (see, for example,Euclid [1] and Vera C. Rubin Observatory Legacy Surveyof Space and Time (LSST) [2]) and Intensity Mapping(IM) (for instance, SKA [3] and HERA [4]). One of thecommon goals shared within the community is to studylate time probes in order to estimate cosmological pa-rameters at a comparable level of precision with respectto the one achieved from CMB dataset. This program ishighly motivated by several reasons. To mention a few,late time surveys provide 3-dimensional catalogs aboutthe observed Universe, hence they may help understand-ing the evolution of the Universe rather than just givinga picture of it, as 2-dimensional datasets (CMB indeed) ∗ [email protected] † bartolomeo.fi[email protected] ‡ [email protected] do . Secondly, the non-linear evolution of the structuresgets enhanced at later times. As a consequence, this al-lows to explore better non-linear scales in the distributionof matter in the Universe, distribution that contains use-ful information about the energy content of the Universeitself (for instance dark matter candidates and massiveneutrinos). Last but not least, the measurements of latetime probes furnish estimations for the cosmological pa-rameters which are (almost) independent of the adoptedmodel, differently from the CMB ones.For what concerns the comparison between differentreconstructions of cosmological parameters, last yearshave shown a discrepancy between the estimation ofthe present Hubble rate H from CMB and late timeprobes. Indeed, CMB measurements provide a value for H = 67 . ± .
54 km s − Mpc − [9], whereas late timeestimations based on local probes, such as SupernovaeIa (SnIa), returns H = 73 . ± . − Mpc − [10](see [11–15] for other estimations): a discrepancy of al-most 5 σ emerges. A priori, both measurements might bequestionable because of the following reasons. CMB esti-mation provides nowadays the most precise measurementof H but, to this end, requires (as mentioned) that a cos-mological model is chosen in order to analyze the spectra.On the other hand, SnIa catalogs for close sources do notneed any cosmological model to infer H but require tocalibrate the relative magnitude of the observed standard Here we underline that non-linearities in CMB lensing due tothe bispectrum of the gravitational potential might turn CMBspectra into 3 dimensional datasets as well [5–8]. a r X i v : . [ a s t r o - ph . C O ] F e b candles wrt some sources whose luminosity is known (forinstance the Cepheid host).This tension seems to survive also when different latetime probes are investigated [16] and this has raised theinterest of the community into the following question:is there any need for new unknown physics to cancel(or at least mild) this discrepancy? In this regard, itmight be that the sound horizon scale r s could be lowerthan the one predicted by our current knowledge of pre-recombination physics in the framework of ΛCDM model[17]. In fact, since the acoustic angular peak measuredby the CMB is proportional to r s H , this would auto-matically raise the value of H detected from the CMBtowards the one inferred from local measurements. How-ever, as discussed in [18], the modification of r s mightbe not enough to completely reabsorb the tension. Ingeneral, the search for new physics which could explainthe discrepancy questions several aspects of the currentCosmic Concordance model. In particular, in [19] a set of7 key assumptions which might be broken to explain thetension has been identified. All these reasons have mo-tivated an intense research activity during the last years(see [20, 21] for overviews, and [22–25], and referencestherein, for a partial coverage of the related literature).We mention also that a very interesting attempt to mea-sure H from CMB lensing, without invoking the knowl-edge of r s , has been done in [26]. This method returnsa r s -independent constraint of H = 73 . ± . − Mpc − , which seems to reabsorb the tension.Given this state-of-art, we will adopt an agnostic ap-proach to the problem. Indeed, within the conservativeframework of ΛCDM model, the question we aim to an-swer is the following: in view of the forthcoming LSSsurveys, is there any theoretical bias which might in-crease the standard deviation in order to alleviate thetension? This question has already been partially facedin the past. In particular, in [27] the effect of velocitydispersion of local SnIa (at redshift lower than 0.1) hasbeen studied and it has been shown that it can introducea further intrinsic error in the local estimation of H of ∼ r s -independent measurement of H is about σ H = 3 km s − Mpc − .To perform the analysis, our starting point is to derivethe analytical formula for the 2-point correlation functionof the luminosity distance-redshift relation. Hence, wewill discuss general aspects of this function, whose inter-ests goes beyond the ones concerning H . In particular, we will investigate numerically the cosmological informa-tion encrypted in the lowest angular multipoles of the2-point correlation functions. To conclude, we will final-ize our analysis by providing forecasted errors of ∼ . ∼ .
01% for LSST in regard of the measure-ment of H . These results are obtained by using linearpower spectrum. However, as we will show, non-linearscales do not dramatically change this scenario. Theircontribution enhances the forecasted errors to ∼ .
1% forboth EDS and LSST. This renders the measurement forhigh-redshift standard candles ideally much more precisethan the one so far discussed in [27] for close sources.The paper is organized as follows. In Sect. I we de-scribe the general method followed to infer the value of H , from higher LSS surveys, thanks to the knowledge ofthe linear luminosity distance-redshift relation. In Sect.II we compute in details the 2-point correlation func-tion for all the relativistic effects involved in the linearluminosity distance-redshift relation, providing generalexpressions for them. In Sect. III, we further give nu-merical details about the spectrum of lower angular mul-tipoles for lensing and doppler effect, and a general treat-ment about the monopole concerning all the effects. InSect. IV, we assume a numerical set of the cosmologicalparameters within the ΛCDM model and estimate the 2-point correlation function effect by effect. Furthermore,we discuss some technical approximations which speed upthe numerical evaluations of the lensing effect. In Sect.V we outline the consequences of our numerical analysisaccording to the specifics of EDS and LSST and discussthe impact of non-linear scales. Finally, in Sect. VI wesummarize our main conclusions. Moreover, in App. Awe furnish technical details about the calculations of the2-point correlation function. In App. B, we report ex-plicit expression for the lower multipoles of the lensingangular spectrum. App. C contains useful properties ofthe spherical Bessel functions, whereas App. D containsdetails about the multipole expansion of lensing 2-pointcorrelation function. I. COSMIC VARIANCE
The starting point for our analysis is the well-knownexpression for the luminosity distance-redshift relation d L ( z ) in the homogenous and isotropic Cosmic Concor-dance model d L ( z ) = 1 + zH (cid:90) z dz (cid:48) (cid:112) Ω m (1 + z (cid:48) ) + Ω , (1)where H is indeed the Hubble constant and Ω m andΩ are respectively the energy density for the Cold Mat-ter and Cosmological Constant today. Within the ΛCDMmodel assumptions, it is well known that Ω = 1 − Ω m such that Eq. (1) contains only two free parameters. Eq.(1) can then be inverted and provides a relation whichcan be used to infer the value of H as H = 1 + zd L ( z ) (cid:90) z dz (cid:48) (cid:112) Ω m (1 + z (cid:48) ) + 1 − Ω m . (2)In this way, given an a priori knowledge of Ω m , inde-pendent estimations for redshift and luminosity distancefrom a given sample of sources can provide a measure-ment of H .Within this framework, the question that we want toaddress is then the following: how precise can in princi-ple be this estimation if we consider the inhomogeneitiesall around our observed Universe? In other words, theobserved luminosity distance-redshift relation d L ( z ) is af-fected by the inhomogeneities and this provides an in-trinsic dispersion for d L ( z ) which is governed by the waythe cosmological structures are distributed and evolve.Hence, in this regard, we consider the observed inhomo-geneous luminosity distance-redshift relation (cid:102) d L ( z ) (cid:102) d L ( z, n ) = d L ( z ) (cid:104) δ (1) ( z, n ) + δ (2) ( z, n ) (cid:105) , (3)where n is the observed direction for the given source and δ (1) and δ (2) are linear and second order corrections tothe luminosity distance-redshift relation. Let us under-line that here we do not have to consider perturbationsin the redshift z . Indeed, the inhomogeneous (cid:102) d L ( z ) isevaluated by construction at constant observed redshifthypersurfaces. This means that redshift corrections arealready taken into account in the (cid:102) d L ( z, n ). From the geo-metrical viewpoint, this choice coincides with slicing thespace-time on constant observed redshift time-like hyper-surfaces and then the time-like gauge mode sets redshiftperturbations null by construction. Since (cid:102) d L ( z ) is an ob-servable, this choice is completely allowed and does notaffect the result.At this point, following [27],we define the inhomoge-neous value of the Hubble constant (cid:102) H as (cid:102) H ≡ z (cid:102) d L ( z ) (cid:90) z dz (cid:48) (cid:112) Ω m (1 + z (cid:48) ) + 1 − Ω m = H d L ( z ) (cid:102) d L ( z )= H (cid:20) − δ (1) − δ (2) + (cid:16) δ (1) (cid:17) (cid:21) ( z, n ) . (4)Eq. (4) contains also pure second order perturbationsof the luminosity distance-redshift relation. Because ofthat, (cid:102) H inferred from the observation of a single source isexpected to deviate from H (see [27] for the case of smallredshift surveys). To estimate this deviation, we shouldselect a prescription of the light-cone average taken allaround the observed sky at fixed redshift which is well-suited for our observables (see [30–33] for the generalclassification of the viable well-posed prescriptions for thelight-cone averages). However, if we consider only twodimensional spheres at constant redshift, the impact of such a measure is null on the estimation of the variance atthe leading order [34]. In addition, also the estimation ofthree dimensional light-cone averages over a redshift binreduces to a two dimensional average for the limit case ofsmall redshift bin (see [33] for the detailed discussion ofthis limit and also [35, 36]). Hence, also in this case theprescription for the average is irrelevant for the leadingorder term of the variance. Because of that, we simplyskip the measure in our formalism and define the light-cone average as (cid:104) . . . (cid:105) ≡ π (cid:90) d Ω( . . . ) , (5)having in mind that this is no longer valid neither forthe evaluation of next-to-leading order contribution norfor the finite-size redshift bin average. On top of that,we also denote with · · · the ensemble average over allthe possible configuration of cosmological perturbations,provided that linear perturbations have a gaussian distri-bution with null mean value . In this way, the variancerelated to the estimation of H is given by σ H ≡ (cid:104) (cid:102) H (cid:105) − (cid:104) (cid:102) H (cid:105) . (6)Hence, from Eq. (4) we get at second order (cid:104) (cid:102) H (cid:105) = H (cid:20) − (cid:104) δ (2) (cid:105) + 3 (cid:104) (cid:0) δ (1) (cid:1) (cid:105) (cid:21) (cid:104) (cid:102) H (cid:105) = H (cid:20) − (cid:104) δ (2) (cid:105) + 2 (cid:104) (cid:0) δ (1) (cid:1) (cid:105) (cid:21) . (7)It then follows from Eq. (6) σ H = H (cid:104) (cid:0) δ (1) (cid:1) (cid:105) . (8)As above-mentioned, the leading order of σ H is entirelygiven by linear perturbation theory. This is in agreementwith [34] and is a quite general result, independent ofthe chosen observable (see also [37] for the applicationof this result related to the estimation of cosmologicalparameters to other cosmological observables).Eq. (8) provides then the intrinsic uncertainty to theestimation of H given by the presence of inhomogeneitiesall around our observed Universe and it is the lowest the-oretical uncertainty we can reach, according to the sam-ple of sources we have access to. This quantity is usuallynamed cosmic variance and quantifies the error for theestimation of H from a single source placed in an idealsurvey of a large number of sources for each constantredshift hypersurface, uniformly distributed all over thesky. In practice, however, all the surveys contain a fi-nite number of sources N which can cover only a partial We will provide analytic version for these assumptions later, inEqs. (18). window of the sky. Given that, following [27], a moreobservationally oriented definition for the estimation of σ H is provided by (cid:104) . . . (cid:105) → N (cid:88) i,j ( . . . ) , (9)where the indices i, j run over all the pairs ( z i , n i ), re-spectively labeling redshift and observed position of the i -th source in the survey. Eq. (9) then adopt the follow-ing prescription for the uncertainty σ H H = 1 N (cid:88) i,j δ (1) ( z i , n i ) δ (1) ( z j , n j ) , (10)rather than Eq. (8). Eq. (10) is the variance associatedto the average value of H inferred from a finite surveyof N sources and corresponds to the locally measuredHubble parameter H from the covariance matrix of the (cid:102) d L ( z ), given an arbitrarily distributed sample of N ob-served sources at positions ( z i , n i ). Indeed, the varianceassociated to the average value of H inferred from a fi-nite survey of N sources is σ H = (cid:32)(cid:88) i (cid:102) H ( z i , n i ) N − H (cid:33) (cid:88) j (cid:102) H ( z j , n j ) N − H = 1 N (cid:88) i,j (cid:16) (cid:102) H ( z i , n i ) (cid:102) H ( z j , n j ) − H (cid:17) = H N (cid:88) i,j δ (1) ( z i , n i ) δ (1) ( z j , n j ) , (11)which precisely corresponds to Eq. (10).In the next section, these general preliminaries will beapplied to the case of linear perturbations of the lumi-nosity distance. This will provide the explicit expressionfor σ H due to all the linear relativistic corrections. II. ANALYTIC RESULTS
In the previous section, we have shown in completegenerality that the cosmic variance σ H is sourced at theleading order only by linear perturbations. To make theexplicit evaluation of all the terms needed for its estima-tion, we first need to consider all the linear relativisticcorrections involved in the δ (1) . To this aim, we onlyconsider linear scalar perturbations in the Longitudinal Gauge without anisotropic stress ds = a ( η ) (cid:8) − (1 + 2 ψ ) dη + (1 − ψ ) (cid:2) dr + r d Ω (cid:3)(cid:9) , (12)and then formally write the linear perturbation of δ (1) [40–43] as δ (1) ( z, n ) = (cid:88) E ˆ O E ψ ( η E , r E n ) , (13)where the index E denotes the sum over the linear rela-tivistic effects, the linear operators areˆ O P V = − Ξ s (cid:90) η s η in dη a ( η ) a ( η o ) ∂ r ( . . . )ˆ O SW = − (1 + Ξ s )( . . . )ˆ O ISW = − s (cid:90) η o η s dη ∂ η ( . . . )ˆ O T D = 2∆ η s (cid:90) η o η s dη ( . . . )ˆ O L = − η s (cid:90) η o η s dη η − η s η o − η ∆ ( . . . ) , (14)where η o is the present conformal time, η s is the confor-mal time of the source, η in is an initial time when pertur-bations were negligible (or, more precisely, the integrandsof the related operators were negligible), ∆ η x = η o − η x (where x can be either s or i ), ∆ is the angular Laplacianand Ξ s = (cid:18) − H s ∆ η s (cid:19) , (15)where H s = ∂ η a ( η s ) /a ( η s ) is the conformal Hubble func-tion. In Eq. (13), η E and r E in ψ depends on whichoperator acts on ψ . Indeed, for the relativistic effects in-tegrated along the observer’s past light-cone, i.e. ISW,TD and L, we have that η E = η and r E = η o − η , so bothof these variables are integrated. On the other hand,for what concerns the PV, the ψ is integrated along thesource world-line and then η E = η whereas r E = η o − η s .This means that for the PV only time is integrated whenthe ˆ O P V acts on ψ . Finally, SW is a local relativisticeffect and then, in this case η E = η s and r E = η o − η s . This assumption might look too restrictive. However in the fol-lowing sections, we will take into account sources located afterthe decoupling, where our assumption works well. The readerinterested in the general expression in presence of anisotropicstress can have look at [38] for the general non-linear expressionof luminosity distance-redshift relation and [39] for the non-linearexpression of redshift containing also the observer terms. Here the subscripts stand for Peculiar Velocity (PV), Sachs-Wolfe (SW), Integrated Sachs-Wolfe (ISW), Time Delay (TD)and Lensing (L). Hence the label E runs in the set (PV, SW,ISW, TD, L). We omit relativistic corrections due to the gravi-tational potential at the observer position. At this point, for a practical evaluation of the ensembleaverage, we move from real to k -space. We then Fouriertransform the gravitational potential as ψ ( η E , r E n ) = 1(2 π ) / (cid:90) d k e i k · n r E g ( η E ) g ( η o ) (cid:101) ψ ( k ) , (16)where g ( η ) is the standard approximated expression ofthe growth function of scalar perturbations in terms ofthe current values of the critical density parameters Ω m and Ω Λ (see e.g. [44]), namely g ( η ) = 52 g ∞ Ω m Ω / m − Ω Λ + (cid:0) Ω m (cid:1) (cid:0) Ω Λ (cid:1) Ω m = Ω m (1 + z ) Ω m (1 + z ) + Ω Λ0 , Ω Λ = Ω Λ0 Ω m (1 + z ) + Ω Λ0 , (17)where Ω m + Ω Λ0 = 1, and where g ∞ is a normalizationconstant fixed such that g ( η o ) = 1. Moreover (cid:101) ψ ( k ) aredelta-correlated functions (cid:101) ψ ( k ) = 0 , (cid:101) ψ ( k ) (cid:101) ψ ( p ) = | ψ k | δ ( k + p ) . (18)In terms of this expansion, linear perturbations in Eq.(13) can be written as δ (1) s = 1(2 π ) / (cid:90) d k (cid:88) E ˆ O E (cid:20) e i k · n r E g ( η E ) g ( η o ) (cid:101) ψ ( k ) (cid:21) . (19)In this way, the combination of Eq. (10) with Eq. (14)gives σ H H = 1 N (cid:88) i,j δ (1) i δ (1) j = 1 N (cid:88) i,j (cid:88) E,E (cid:48) ˆ O Ei ˆ O E (cid:48) j ψ ( η Ei , r Ei n i ) ψ ( η E (cid:48) j , r E (cid:48) j n j ) , (20)where the index i, j run over the sources and the index E, E (cid:48) run over all the effects for each source. Hence, byinserting the expansion (16) in Eq. (20) and using Eq.(18) we obtain σ H H = 1 N (cid:88) i,j (cid:88) E,E (cid:48) (cid:90) dkk P ψ ( k ) W Ei,E (cid:48) j , (21)where we have defined W Ei,E (cid:48) j ≡ π (cid:90) d Ω k ˆ O Ei ˆ O E (cid:48) j × (cid:20) g ( η Ei ) g ( η o ) g ( η E (cid:48) j ) g ( η o ) e i k · ( r Ei n i − r E (cid:48) j n j ) (cid:21) , (22) with Ω k to be meant as the solid angle in k -space and weused the so-called dimensionless power spectrum P ψ ( k ) ≡ k π | ψ k | . (23)The variance in Eq. (21) can be written as a sum ofthe contribution over different pairs of effects. Since wehave 5 different effects we will find 15 different pairs ofeffects. We distinguish between the contribution of theeffects of the same kind, which we refer to as pure terms and of the effects of different kinds, which we refer to as mixed terms σ H H = 1 N (cid:88) i,j (cid:88) E ξ E ( z i , z j , n i · n j )+ 1 N (cid:88) i,j (cid:88) E (cid:54) = E (cid:48) ξ EE (cid:48) ( z i , z j , n i · n j ) , (24)where we have defined ξ EE (cid:48) ( z i , z j , n i · n j ) = (cid:90) dkk P ψ ( k ) W EE (cid:48) ij ,ξ E ( z i , z j , n i · n j ) = ξ EE ( z i , z j , n i · n j ) . (25) ξ EE (cid:48) ( z i , z j , n i · n j ) are nothing but the 2-point correla-tion functions between the relativistic effects E and E (cid:48) evaluated for two different sources with redshifts z i and z j along the observed directions n i and n j . We concludethis section with the explicit expressions for all the W EE (cid:48) .In particular, we define G i = (cid:90) η i η in dη a ( η ) a ( η i ) g ( η ) g ( η o ) R ( η x , η y , ν ) = (cid:113) ∆ η x + ∆ η x − η x ∆ η y νL ( η x , η y , ν ) = ∆ η x ∆ η y νR ( η x , η y , ν ) H ( η x , η y , ν ) = ∆ η x ∆ η y √ − ν R ( η x , η y , ν ) , (26)where R is the distance between two sources and L and H are respectively the normalized scalar and (modulo ofthe) vector products between the two directions of thesources. We then have that the pure terms are W P V ij = Ξ i Ξ j G i G j k (cid:26) ∆ η i ∆ η j (1 − ν ) R j ( kR ) + ν j ( kR ) − j ( kR )] (cid:27) ( η i , η j , ν ) W SW ij =(1 + Ξ i )(1 + Ξ j ) g ( η i ) g ( η o ) g ( η j ) g ( η o ) j ( kR ( η i , η j , ν )) W ISW ij =4 Ξ i Ξ j (cid:90) η o η i dη (cid:90) η o η j dη (cid:48) ∂ η g ( η ) g ( η o ) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) j ( kR ( η, η (cid:48) , ν )) W T Dij = 4∆ η i ∆ η j (cid:90) η o η i dη (cid:90) η o η j dη (cid:48) g ( η ) g ( η o ) g ( η (cid:48) ) g ( η o ) j ( kR ( η, η (cid:48) , ν )) W Lij = 1∆ η i η j (cid:90) η o η i dη η − η i η o − η (cid:90) η o η j dη (cid:48) η (cid:48) − η j η o − η (cid:48) g ( η ) g ( η (cid:48) ) g ( η o ) (cid:104) k H j ( kR ) − k H L j ( kR )+ k (cid:0) L − H (cid:1) j ( kR ) + 4 k L j ( kR ) (cid:105) ( η, η (cid:48) , ν ) , (27)where ν ≡ n i · n j , j n are the spherical Bessel functions of n -th order and, in the same way, the mixed terms are W P V i,Lj = Ξ i ∆ η j G i (cid:90) η o η j dη η − η j η o − η g ( η ) g ( η o ) (cid:26) − k ∆ η (∆ η i − ν ∆ η )(∆ η − ν ∆ η i ) R j ( kR )+ k ∆ η (cid:18) η ∆ η i − ν ∆ ηR − ν (cid:19) j ( kR ) − k ∆ ηR (cid:2) k ∆ η (∆ η i − ν ∆ η ) − ν (cid:3) j ( kR ) (cid:27) ( η i , η, ν ) W P V i,SW j =Ξ i (1 + Ξ j ) G i g ( η j ) g ( η o ) k ( ν ∆ η j − ∆ η i ) (cid:18) j ( kR ) R (cid:19) ( η i , η j , ν ) W P V i,ISW j =2 Ξ i Ξ j G i (cid:90) η o η j dη ∂ η g ( η ) g ( η o ) k ( ν ∆ η − ∆ η i ) (cid:18) j ( kR ) R (cid:19) ( η i , η, ν ) W P V i,T Dj = − i ∆ η j G i (cid:90) η o η j dη g ( η ) g ( η o ) k ( ν ∆ η − ∆ η i ) (cid:18) j ( kR ) R (cid:19) ( η i , η, ν ) W Li,SW j = 1 + Ξ j ∆ η i g ( η j ) g ( η o ) (cid:90) η o η i dη η − η i η o − η g ( η ) g ( η o ) (cid:20) k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) (cid:21) ( η, η j , ν ) W Li,ISW j =2 Ξ j ∆ η i (cid:90) η o η i dη η − η i η o − η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) (cid:20) k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) (cid:21) ( η, η (cid:48) , ν ) W Li,T Dj = − η i ∆ η i (cid:90) η o η i dη η − η i η o − η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) (cid:20) k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) (cid:21) ( η, η (cid:48) , ν ) W SW i,ISW j =2 (1 + Ξ i )Ξ j g ( η i ) g ( η o ) (cid:90) η o η j dη ∂ η g ( η ) g ( η o ) j ( kR ( η i , η, ν )) W SW i,T Dj = − i ∆ η j g ( η i ) g ( η o ) (cid:90) η o η j dη g ( η ) g ( η o ) j ( kR ( η i , η, ν )) W ISW i,T Dj = − i ∆ η j (cid:90) η o η i dη ∂ η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) j ( kR ( η, η (cid:48) , ν )) . (28)The detailed analytic derivation for these 15 differentcontributions is reported in App. A and follows thederivation obtained in [45]. We just remark that the ker-nel W P V ij is in agreement with the one found in [27].
III. MULTIPOLES ANALYSIS
Let us now investigate the angular decomposition ofthe 2-point correlation function ξ in Eqs. (25). To thisend, we expand the angular dependence in multipoles as ξ EE (cid:48) ( z , z , ν ) = ∞ (cid:88) (cid:96) =0 C EE (cid:48) (cid:96) ( z , z ) P (cid:96) ( ν ) (29)where P (cid:96) ( x ) are the Legendre polynomials of order (cid:96) and C EE (cid:48) (cid:96) ( z , z ) = 2 (cid:96) + 12 (cid:90) − dν ξ ( z , z , ν ) P (cid:96) ( ν ) . (30)With this decomposition, we aim to investigate the be-havior of low multipoles, in order to understand how fastthe truncated version of Eq. (29) converges to the fullnumerical results. To the extent of this paper, we justlimit our analysis to the lensing and doppler terms inthe 2-point correlation function, but this analysis can beapplied to all the effects. A. Lensing
We start by analyzing in detail the 2-point correlationfunction for the lensing. Hence, we apply Eq. (29) to theterm ξ L ( z , z , ν ) = (cid:90) dkk P ψ ( k ) W L ( z , z , ν ) (31)where W L is the lensing kernel for the 2-point correla-tion function as reported in Eqs. (27). First of all, wenotice that the monopole is exactly C = 0, regardlessof the chosen redshifts (see Appendix D for details). Inthe ideal case of infinite number of sources densely dis-tributed in each redshift bin all over the sky, the statis-tical average tends to the monopole. Hence, in this case,lensing is not expected to affect the variance of H at all.It is interesting to notice that this property about the monopole stands for all the cross-correlation terms be-tween lensing and other effects, as we will prove in Eqs.(44) and (45). This means that the above-mentionedideal case is not affected at all by the expected leadingcorrection due to lensing.However, realistic surveys deals with partial sky cover-ages. As discussed before, this sky coverage is very lim-ited for realistic forthcoming surveys (see, for instance,Euclid Deep Survey and LSST). Hence, the effect ofhigher multipoles is expected to contribute to the vari-ance for realistic surveys. This is indeed due to thefact that window function introduced by the partial sky-coverage is convolved in (cid:96) -space with the higher multi-poles and then an amount of power is transferred fromhigher multipoles to the monopole itself. In particular,each multipoles can be written as an integral in k -spaceas C L(cid:96) = (cid:90) dkk P ψ ( k ) L L(cid:96) ( z , z , k ) , (32)where the kernel of the integrand L L(cid:96) is L L(cid:96) = 2 (cid:96) + 12 (cid:90) − dν W L ( z , z , ν ) P (cid:96) ( ν ) . (33)It is hard to solve Eq. (33) analytically. However, wehave outlined the following approximation scheme to dealwith it. We first perform variable changes x = η o − η and y = η o − η (cid:48) into the integrals in W L in Eq. (27). Thisleads to W Lij = 1∆ η i η j (cid:90) ∆ η i dx ∆ η i − xx (cid:90) ∆ η j dy ∆ η j − yy g ( η o − x ) g ( η o − y ) g ( η o ) (cid:104) k H j ( kR ) − k H L j ( kR )+ k (cid:0) L − H (cid:1) j ( kR ) + 4 k L j ( kR ) (cid:105) ( x, y, ν ) , (34)where the geometrical functions R , L and H now simplifyto R ( x, y, ν ) = (cid:112) x + y − xyνL ( x, y, ν ) = x y νR ( x, y, ν ) H ( x, y, ν ) = x y √ − ν R ( x, y, ν ) . (35) Hence we adopt a polynomial expansion of the growthfunction g as g ( η o − x ) = (cid:80) i =0 g i x i , where the preci-sion of the approximation is shown in Fig. 1 and thecoefficient g i are reported in the related caption. As aconsequence, Eq. (34) as well can be expanded as a sumof polynomial terms as W Lij = 1∆ η i η j (cid:88) n,m =0 g n g m g ( η o ) (cid:90) ∆ η i dx ∆ η i − xx (cid:90) ∆ η j dy ∆ η j − yy x n y m (cid:104) k H j ( kR ) − k H L j ( kR )+ k (cid:0) L − H (cid:1) j ( kR ) + 4 k L j ( kR ) (cid:105) ( x, y, ν ) , (36) - - - - - - - Z Δ g g FIG. 1. Relative error between the exact numerical solutionfor the growth function in term z and its approximated poly-nomial expression with coefficient g = 1 , g = 1 . × − , g = − . × − , g = 1 . × − , g = − . × − and g = 1 . × − . Solid lines refer topositive values and dashed lines stand for negative ones. Aswe can see from this plot, our polynomial approximation for g ( z ) is precise at 0 .
01% level.
Thanks to this trick, integrals over x, y and ν in Eq.(32) can be done analytically and the computation of the C L(cid:96) ’s eventually requires a single integration in k -spaceleft. On one side, this polynomial expansion allows to getanalytic results for the integrand of (33). On the otherside, this allows to lower the computational time for eachmultipoles and increase the numerical precision. More-over, thanks to the angular integration, integrals over x and y factorize. We show this feature for the dipole (cid:96) = 1,since the monopole (cid:96) = 0 is null, as above-mentioned.Indeed, after a long but straightforward calculation, theangular integration in Eq. (33) for (cid:96) = 1 gives L L = 12 k ∆ η i ∆ η j (cid:88) n,m =0 g n g m g ( η o ) (cid:90) ∆ η i dx ∆ η i − xx × (cid:90) ∆ η j dy ∆ η j − yy x n − y m − ( kx cos kx − sin kx ) × ( ky cos ky − sin ky )= 12∆ η i ∆ η j (cid:88) n,m =0 g n g m g ( η o ) Q n (∆ η i ) Q m (∆ η j ) , (37)where Q n ( z ) = − (cid:90) z dx ( z − x ) x n − j ( kx ) . (38) Q n are nothing but integrals of the form (cid:82) dx x α sin( kx )and (cid:82) dx x α cos( kx ) which are analytically solvable withmultiple integrations by part. The same evaluation per-formed for first 6 multipoles shows that we can generalizeEq. (37) for (cid:96) ≤ L L(cid:96) = (2 (cid:96) + 1) (cid:96) ( (cid:96) + 1) ∆ η i ∆ η j (cid:90) ∆ η i dx ∆ η i − xx × (cid:90) ∆ η j dy ∆ η j − yy g ( η o − x ) g ( η o − y ) g ( η o ) j (cid:96) ( kx ) j (cid:96) ( ky )= (2 (cid:96) + 1) (cid:96) ( (cid:96) + 1) ∆ η i ∆ η j (cid:88) n,m =0 g n g m g ( η o ) Q (cid:96)n (∆ η i ) Q (cid:96)m (∆ η j ) , (39) where Q (cid:96)n trivially generalizes Eq. (38) as Q (cid:96)n ( z ) = − (cid:90) z dx ( z − x ) x n − j (cid:96) ( kx ) . (40)In App. B we explicitly evaluate Q (cid:96)n ( z ) for the lowermultipoles (cid:96) = 1 , ,
3. Eq. (39) suggests that it mightbe generalized for every (cid:96) . So far, the general proof forevery (cid:96) is still missing, since we have only checked (cid:96) by (cid:96) its validity for the first 6 multipoles, also plotted inFig. 2. However, we notice that first line of Eq. (39)is consistent with the evaluations done in [40] about themultipoles expansion of the luminosity distance 2-pointcorrelation function.To show this point, we first remark that Eq. (39) takesinto account the whole contribution of the angular lapla-cian in the lensing kernel, which contains also sublead-ing terms in the counting of power of k (see Eqs. (A6)and (A7)). On the contrary, in [40] the lensing termscontains only leading terms according to the weight in k -space. Hence only first two terms in our Eq. (A6) aretaken into account in the term C (5) (cid:96) of [40], whereas thelast term in our expansion (A6) is accounted for in theterm C (3) (cid:96) of [40]. Beside this different classification, thismeans that our Eq. (39) should be compared with theappropriate combination of kernels of Eqs. (70), (73) and(74) of [40], according to our Eq. (A6), rather than justtheir expression for C (5) (cid:96) . Once this subtle point is takeninto account, our evaluations agree with [40]. Since theresult of [40] are valid for any (cid:96) , this supports the factthat our Eq. (39) can be considered for any (cid:96) .We remark the appearance of the prefactor (cid:96) ( (cid:96) + 1) in Eq. (39). This is due to the presence of ∆ in eachˆ O L involved in the lensing 2-point correlation function.The fact that monopole vanishes then directly followsfrom the fact that the eigenvalue of P for ˆ O L is 0 (seeAppendix D for details).Beside the technical aspects, it is easy to study thebehavior of L (cid:96) ’s themselves, such that they can be inter-preted as power spectra for each multipole. In Fig. 2we plot these power spectra for the first six multipolesfor the lensing-lensing correlation at the same redshift,which ranges from z = 0 .
15 (bluer) to z = 1 .
55 (redder). × - × - × - × - × - × - k [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) × - × - × - × - × - k [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) × - × - × - × - × - × - k [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) × - × - × - × - × - × - × - k [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) × - × - × - × - × - × - × - × - k [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) × - × - × - × - × - k [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) FIG. 2. Power spectrum for the different multipoles of the 2-point correlation function of the lensing term. The lowestmultipoles are considered: (cid:96) = 1 (top-left), (cid:96) = 2 (top-right), (cid:96) = 3 (center-left), (cid:96) = 4 (center-right), (cid:96) = 5 (bottom-left) and (cid:96) = 6 (bottom-right). The monopole is not shown as it is null for the correlation function of the lensing correction. All theselines refers to the same redshift in the correlation function, ranging from z = 0 .
15 (bottom curves) to z = 1 .
55 (top curves),with an interval step of 0.1. As we can notice, higher multipoles exhibit the presence of a second peak which becomes moreimportant at higher redshifts. These plots show scales ranging from k IR = H = 3 × − h Mpc − to k UV = 0 . h Mpc − . Along these plots, we notice several interesting fea-tures: • First of all, higher redshifts exhibit higher peaksthan closer sources. This is somehow expected,since this behavior just reflects the fact that lens-ing is an integrated effect, such that the farer thesource, the bigger the effect. • Secondly, for higher multipoles the peak is driftedtowards the smaller scales. This is also expected since higher multipoles investigates the correlationon smaller angular scales. To be mentioned is thebehavior for the dipole. Here at larger redshift wenotice that a considerable amount of power lies onsuper-Hubble scales. This leads to a degeneracybetween super-Hubble fluctuations and the choiceof the background value. This is due to the factthat we have not taken into account corrections dueto the observer’s peculiar motion. Indeed, the lat-ter contributes to the dipole and this renders the0total 2-point correlation function independent ofthe physics beyond super-Hubbles scales, restoringcausality as shown in [46]. However, as it emergesfrom Figs. 2, the lensing dipole is comparable toother multipoles only for very low redshifts. Hence,truncation of the angular spectra at (cid:96) = 2 intro-duces a negligible error to the lensing contributionto the 2-point correlation function for high redshifts(see also Fig. 4). This method is in line with whatsuggested in [47] and has the practical advantageof getting rid of the dependence on super-Hubblephysics and lower the number of terms involved inthe analysis. • Moreover, we notice the appearance of a secondarypeak at higher multipoles. The relative importanceof the second peak becomes higher for farer sources.This is shown in Fig. 3, where the same plot asin Fig. 2 is done for (cid:96) = 6 and where the redshiftranges from z = 1 (bottom) to z = 7 (top). It is evi-dent how at higher redshift the peaks are such com-petitive and close that it sounds fairer to refer toa range of scales as the dominant ones rather thanjust a single peak. According to our understanding,this behavior reflects the fact that on larger multi-poles and at higher redshifts, at the given angularscale, there is room for ”resonances” in the k -spaceon the transverse plane wrt to the line-of sight. × - × - [ h Mpc - ] k - Ψ ( k ) ℒ ( z , z , k ) FIG. 3. The same as in Fig. 2 for (cid:96) = 6. Now the redshiftranges from z = 1 (bottom curve) to z = 7 (top curve) witha redshift gap of 1. As we can notice, at higher redshift theappearance of a second peak eventually flattens the curves attheir maximum. Figs. 2 also indicate that the (cid:96) -expansion convergesquite slowly to the full angular correlation function. In-deed, until (cid:96) = 6 the order of magnitude of the spectratends to increase or stays constant for the same redshift.This behavior is more evident in the top panel of Fig. 4.Here we plot the first 19 multipoles for the lensing andconsider the cross correlation between the sources at red-shift z = 0 .
75 and other sources with possibly differentredshift, letting z range from 0 .
15 (bluer, bottom) to1 .
55 (redder, top). As we can appreciate from the figure, the order of magnitude for the multipoles at the sameredshift tends to remain the same or increase. For whatconcerns the pair ( z , z ) = (0 . , .
75) (bottom curvein Fig. 4), the multipoles decrease in amplitude after (cid:96) ∼
7. However, this decreasing is not rapid enough toallow convergence to the full result for the chosen first19 multipoles. To make this point clearer, in the bottom × - × - × - × - × - × - ℓ C ℓ ( z , z ) × - × - × - θ ξ l en s ( z , z , θ ) FIG. 4. Top-panel: 19 lowest multipoles of the 2-point angu-lar correlation function of lensing. Here the first redshift isfixed at z = 0 .
75 whereas the second one z varies from 0 . .
55 (top curve), with a redshift step of0 .
1. Bottom-panel: comparison between the numerical evalu-ation of the 2-point angular correlation function for the lens-ing (blue points) with the semi-analytical estimation givenby the sum of the multipoles truncated at (cid:96) = 19 (red curve).Here z = 0 .
75 and z = 0 .
15. The bottom panel correspondsto the re-summation of multipoles of the bottom curve in thetop panel. The comparison between the two figures makesclear that the multipoles for the considered redshift slowlydecreases and the cutoff to be chosen in order to truncatethe multipole expansion is not yet reached at (cid:96) = 19. Thisbehavior is due to the steepness of the peak at ν = 1. Theconvergence becomes even slower at higher redshift, as shownin the other curves in the top panel. panel of Fig. 4 we compare the full numerical estimationof the the angular correlation function at z = 0 .
15 and z = 0 .
75 (blue points) with its truncated (cid:96) -expansionuntil (cid:96) = 19 (with coefficient given in the bottom curveof the top panel). We realize from this plot that inter-mediate angular ranges are in an acceptable agreement(modulo an oscillation due to the truncation in the Leg-endre expansion), whereas the maximum still admits a1 ∼
15% missing correction. Considering that the numeri-cal value of the peak is 3 × − and the order of magni-tude of the C (cid:96) for this bin pair is 10 − , we can estimatethe convergence (cid:96) scale thanks to the argument whichfollows. The value of the angular correlation at its peakcorresponds to the full resummed series of all the C (cid:96) ’s fora given redshift pair. This can be understood by lookingat the Eq. (30) since Legendre’s polynomials all equal1 when evaluated at ν = 1. Thanks to a simple ana-lytical estimation, we can then infer that a 1% error onthe peak’s resolution for the chosen case will be reachedaround at (cid:96) ∼
50. The situation is even worse for theother redshift pairs.Finally, the high value of these convergent angularscales means that the angular correlation function forthe lensing is extremely peaked around its maximumand then tends to correlate between light-signals emittedfrom sources separated by a very narrow angular size.
B. Doppler
The same analysis can be done for the Doppler angu-lar correlation function. The analogous of L L(cid:96) for theDoppler effect can be written exactly for each multi-poles without requiring the polynomial interpolation for g . This is due to the fact that Doppler effect is not in-tegrated along the line-of-sight and then the integrationsof the growth functions factorizes in the ultimate expres-sion in the functions G i in W P V in Eq. (27). This leadsto an important difference when compared with lensing.Indeed, for the case of Doppler, L P V(cid:96) can be written as aproduct of combination of spherical Bessel functions foreach (cid:96) , and this product is symmetric under the exchangeof z i ↔ z j . Moreover, the dependences on the two red-shifts are decoupled. For a matter of clarity, we reportthe simplest cases of the monopole ( (cid:96) = 0) and dipole( (cid:96) = 1) for the Doppler L P V ( z i , z j ) = j ( k ∆ η i ) j ( k ∆ η j ) L P V ( z i , z j ) = 3 (cid:20) j ( k ∆ η i ) k ∆ η i − j ( k ∆ η i ) (cid:21) × (cid:20) j ( k ∆ η j ) k ∆ η j − j ( k ∆ η j ) (cid:21) . (41)It is worth to remark that Eq. (41) can be writtenas L P V(cid:96) = (2 (cid:96) + 1) j (cid:48) (cid:96) ( k ∆ η i ) j (cid:48) (cid:96) ( k ∆ η j ). This is in agree-ment with results of [40]. The same structure occurs forhigher multipoles, but with a more involved combina-tion of higher order j n . Since the product of j n ’s showsa constructive interference only when z ≈ z , L P V(cid:96) ismostly positive only when the two redshift are almost thesame. For larger separation, L P V(cid:96) exhibits an oscillatingbehavior which suppresses the integration in k -space andreduces the amplitude of multipoles. This explain whyDoppler angular correlation is relevant only for sourcesplaced within the same redshift bin. Finally, we com-ment on the fact that this is not the case for lensing, when the L L(cid:96) ’s are positive defined also for large angularseparation. This is due to the fact that lensing is an in-tegrated effect along the line-of-sight and then also twosources with a separation in redshift space contribute tothe amplitude in a non-negligible way.
C. Monopoles for the general 2-point correlationfunction
The above-mentioned features for the spectral decom-position of the 2-point correlation functions can be gen-eralized by evaluating the angular integral over Ω k in Eq.(22). In fact, since the operators ˆ O E do not act on theangular dependences in d Ω k in the kernels W Ei,E (cid:48) j , wecan perform this integration independently by aligningthe z axis of k with the vector r Ei n i − r E (cid:48) j n j . We thenget (cid:90) d Ω k e i k · ( r Ei n i − r E (cid:48) j n j ) =2 π (cid:90) − d (cos θ k ) e ikR cos θ k =4 πj ( kR ) , (42)where R = R ( η E , η E (cid:48) , ν ). In this way, Eq. (22) gets thesimple expression W Ei,E (cid:48) j = ˆ O Ei ˆ O E (cid:48) j (cid:20) g ( η Ei ) g ( η o ) g ( η E (cid:48) j ) g ( η o ) j ( kR ) (cid:21) . (43)This exact expression is as simple as powerful. Indeed,thanks to Eq. (43), we can generate the kernel for anyrelativistic effect E simply as the action of the operatorˆ O E on the function g ( η E ) j ( kR ). Eq. (43) is also usefulsince it provides a generating function for the monopoleof the ξ EE (cid:48) . Indeed, the kernel for the monopole is givenby the action of the angular average on the Eq. (43) W Ei,E (cid:48) j = 12 (cid:90) − dν W Ei,E (cid:48) j = 12 (cid:90) − dν ˆ O Ei ˆ O E (cid:48) j (cid:20) g ( η Ei ) g ( η o ) g ( η E (cid:48) j ) g ( η o ) j ( kR ) (cid:21) = 12 ˆ O Ei ˆ O E (cid:48) j (cid:20) g ( η Ei ) g ( η o ) g ( η E (cid:48) j ) g ( η o ) (cid:90) − dν j ( kR ) (cid:21) = ˆ O Ei (cid:20) g ( η Ei ) g ( η o ) j ( k ∆ η Ei ) (cid:21) × ˆ O E (cid:48) j (cid:20) g ( η E (cid:48) j ) g ( η o ) j ( k ∆ η E (cid:48) j ) (cid:21) , (44)where we used (cid:82) − dνj ( kR ( x, y, ν )) = 2 j ( kx ) j ( ky ).The second-last equality in Eq. (44) is due to the factthat all the ˆ O E ’s commute with the angular average. In-deed, for what concerns the lensing ˆ O L = ∆ , on onehand we recall that the angular average of the angularLaplacian is null. On the other hand, the angular averageof a given function does not depend on the angle, thenits Laplacian is null. Provided that, we can commute ˆ O L O P V (cid:20) g ( η s ) g ( η o ) j ( k ∆ η s ) (cid:21) = k Ξ s G s j ( k ∆ η s )ˆ O SW (cid:20) g ( η s ) g ( η o ) j ( k ∆ η s ) (cid:21) = − (1 + Ξ s ) g ( η s ) g ( η o ) j ( k ∆ η s )ˆ O ISW (cid:20) g ( η ) g ( η o ) j ( k ∆ η ) (cid:21) = − s (cid:90) η o η s dη (cid:18) g (cid:48) ( η ) g ( η o ) j ( k ∆ η ) + g ( η ) g ( η o ) j ( k ∆ η ) (cid:19) ˆ O T D (cid:20) g ( η ) g ( η o ) j ( k ∆ η ) (cid:21) = 2∆ η s (cid:90) η o η s dη g ( η ) g ( η o ) j ( k ∆ η )ˆ O L (cid:20) g ( η ) g ( η o ) j ( k ∆ η ) (cid:21) =0 . (45)These relations are enough to evaluate the monopole forthe whole ξ . We remark that lensing has no monopole.Hence, from Eqs. (45) and the factorization in Eq. (44)we infer that also all the cross-correlations between lens-ing and any other relativistic effect in the luminosity dis-tance have vanishing monopole. This has consequencesfor what concerns the observation of large number ofsources distributed all over the sky. Indeed, in the limitof large number of sources, the spatial average tends tothe angular ones and then σ H is completely given by thesum of the monopoles for each combination of EE (cid:48) . Inthis limit, then, any cross-correlation with lensing can-not contribute to the total effect. This independence onlensing is an interesting point, since this is usually one ofthe most important effect in the analysis of Large ScaleStructure surveys.Total monopole is then expected to be dominated bythe auto correlation of Doppler. This indeed is due tothe fact that Eqs. (45) exhibits a further k in the firstline. The latter then contributes with a k amplitude inthe integration over k -space which then contributes morethan the other pairs.We underline that this features of monopoles followsfrom the fact that P is constant. The same analysis canbe applied to higher multipoles. However, since P (cid:96) ’s ingeneral depend on ν , they do not commute with all theoperators ˆ O E (in particular when E = L ) and this aspectmust be taken into account for higher multipoles. Thecosmological information within LSS contained in highermultipoles has been discussed in recent papers as [48, 49],where the dipolar structure of galaxy number counts hasbeen investigated. We postpone the investigation of theseeffects multipole by multipole within forthcoming surveysto future works. We stress, however, that the analysisperformed in the next sections captures the presence of all the multipoles for the 2-point correlation function.We then conclude this section by remarking that typi-cal surveys have no access to the full sky coverage. Thispractical limitation makes lensing contribution no longervanishing and then dominating the other effects. In thefollowing sections, we will show this point and forecastfor some cases of interest for forthcoming surveys the ex-pected values for σ H . IV. NUMERICAL EVALUATIONS
In this section, we will provide numerical evaluationof the 2-point correlation function ξ . To this aim, wewill focus our analysis only on lensing and doppler terms.Indeed, from Eqs. (27) we have that the kernels regardinglensing and doppler effects, namely W L and W P V , havethe higher number of powers in terms of k , respectively k and k . Moreover, W P V contains the pre-factors Ξand these amplify the effect on small redshift. Motivatedby the same argument, in Eqs. (28) the only term ofinterest for us is W P V L which indeed contains the highestnumber of power in k and is amplified by the prefactor Ξ.To this aim, we consider the linear dimensionless powerspectrum today P ψ = A (cid:18) kk (cid:19) n s − (cid:20) g ( η o ) g ∞ (cid:21) T (cid:18) k . k eq (cid:19) , (46)where T ( k ) is the so-called transfer function which takesinto account the sub-horizon evolution of modes re-entering the horizon during the radiation era. We haveexpressed T ( k ) in the Hu-Eisenstein parametrization [50],3given by T ( q ) = L ( q ) L ( q ) + q C ( q ) ,L ( q ) = log(2 e + 1 . q ) ,C ( q ) =14 . . q . (47)We have integrated over the spectral distribution of fre-quency modes using the following infrared (IR) and ul-traviolet (UV) cutoff values: k IR = 3 × − h Mpc − , k UV = 0 . × h Mpc − . (48) They roughly correspond to the present horizon scale andto the limiting scale of the linear spectral regime, respec-tively. The numerical values of the parameters appearingin Eqs. (17), (46) and (47) have been chosen, accordingto recent cosmological observations [9], as follows A = 2 . × − , n s = 0 . , k = 0 .
05 Mpc − ,k eq =0 . h Ω m , h = 0 . , Ω m = 0 . . (49)With these numerical specifications, we first want to plotthe values for aligned ( ν = 1) and antipodes ( ν = − ξ L , ξ P V and ξ P V,L . To this aim, weunderline that these explicit limits ν = ± W L , W P V and W P V, L in Eqs.(27), (28) and (34). We find indeed that those kernelsbecome W ± P V ij = ±
13 Ξ i Ξ j G i G j k [ j − j ] ( k (∆ η i ∓ ∆ η j )) W ± Lij = 1∆ η i η j (cid:90) ∆ η i dx ∆ η i − xx (cid:90) ∆ η j dy ∆ η j − yy g ( η o − x ) g ( η o − y ) g ( η o ) (cid:104) k x y ( x ∓ y ) j ± k x yx ∓ y j (cid:105) ( k ( x ∓ y )) W ± P V i,Lj = Ξ i ∆ η j G i (cid:90) ∆ η j dx ∆ η j − xx g ( η o − x ) g ( η o ) (cid:40) − k x j + k x (cid:18) x ∆ η i ∓ x ∓ (cid:19) j − k x ∆ η i ∓ x (cid:2) k x (∆ η i ∓ x ) ∓ (cid:3) j (cid:41) ( k (∆ η i ∓ x )) , (50)where ± stands for ν = ± x = η o − η in the in-tegration variable. A further simplification can be donefor both W ± L and W ± P V,L about the growth function g .Just as done for the multipoles in the previous section,we can analytically perform the integrals along the line-of-sight thanks to the polynomial expansion of g ( η o − x ).This decreases the computational time, since it reducesthe evaluation of the function to only 1-dimensional nu-merical integration in k -space. Results are shown in Fig.5. Numerical evaluation in the range of redshift of inter-est for forthcoming surveys (0 . ≤ z ≤ .
85) show thatlensing and doppler are competitive effects for redshiftsmaller than 1. Moreover, cross-correlations betweenthese two effects are always negligible in the explored ranges of z . According to this, one might conclude thatthe contributions due to doppler correction are impor-tant. However, doppler terms are counterbalanced bythe change in sign which occurs around z = 1 . − H ∆ η in the expression of ξ EE (cid:48) . This co-efficient is null precisely when ∆ η = H − and changessign. According to our chosen cosmology in Eq. (49),this switch happens exactly at z = 1 . ξ Lij is highly peaked when ν ≈ ν = −
1. The angular scales ν ijth at which this decay occurs can beestimated by the following approximation method. Webuild an approximated ξ Lijapp as ξ Lijapp ( ν ) = Θ( ν ijth − ν ) (cid:12)(cid:12)(cid:12) ξ Lij + (cid:12)(cid:12)(cid:12) − Θ( ν − ν ijth ) (cid:12)(cid:12)(cid:12) ξ Lij − (cid:12)(cid:12)(cid:12) , (51)4 - × - - × - - × - - × - - × - × - × - × - - × - - × - - × - - × - - × - - × - × - FIG. 5. Aligned-correlations ( ν = 1, right panels) and antipodes-correlations ( ν = −
1, left panels) as function of redshift. Thinsolid lines refers to positive values as specified by the relative label. In the same way, dashed lines refers to negative values.Thick solid lines refers to the redshifts where correlations are null. We notice that the effect involving doppler terms (middleand bottom ones) always exhibit a null value when the redshift for the doppler effect is z ≈ .
6. This value corresponds to thescale where Ξ s = 0, namely ∆ η s = H − s . Figures refer to lensing (top panels), doppler (middle panels) and cross-correlationbetween lensing and doppler (bottom panels). In the regime of interest for us, lensing is always leading with respect to theother terms for the aligned-correlations, whereas doppler turns out to be the leading effect for the antipodes-correlation atredshift smaller than 1. Plots for doppler effect terms are shown in logarithmic scales in order to show better the scaling oftheir values. where Θ( x ) is the Heaviside step function and subscripts in ξ + and ξ − respectively indicate ν = +1 and ν = − ξ Lij shown in Fig.5. Moreover, ν ijth is determined by the further require-ment that the angular average over the full sky of ξ Lijapp isnull, just as what is analytically shown for ξ Lij , namely (cid:90) − dν ξ Lijapp ( ν ) = 0 . (52)Indeed, by combining Eqs. (51) and (52), we have that (cid:16) − ν ijth (cid:17) (cid:12)(cid:12)(cid:12) ξ Lij + (cid:12)(cid:12)(cid:12) = (cid:16) ν ijth + 1 (cid:17) (cid:12)(cid:12)(cid:12) ξ Lij − (cid:12)(cid:12)(cid:12) (53)which then returns ν ijth = (cid:12)(cid:12)(cid:12) ξ Lij + (cid:12)(cid:12)(cid:12) − (cid:12)(cid:12)(cid:12) ξ Lij − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ Lij + (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) ξ Lij − (cid:12)(cid:12)(cid:12) . (54)From top panels of Fig. 5, we notice that (cid:12)(cid:12)(cid:12) ξ Lij + (cid:12)(cid:12)(cid:12) (cid:29) (cid:12)(cid:12)(cid:12) ξ Lij − (cid:12)(cid:12)(cid:12) so ν ijth is expected to be very close but lower than 1.In Fig. 6 we show the numerical values for the range ofredshift of our interest. This gives us a threshold angular FIG. 6. Plot of ν ijth for the lensing 2-point correlation function.As we can see, ξ Lij rapidly peaks at values close to ν = 1. scale. Redshift by redshift, beyond this angular scale thelensing 2-point correlation function for a given pair ofsources is negligible.The previous estimation is very powerful, since ν ijth can be directly compared with the angular sky cover-age of a given survey and tell whether the lensing effectis important or not for the chosen catalog. As an in-stance, for the EDS configuration, the angular sky cov-erage is approximatively 20 deg per observation patch[1]. This solid angle can be associated to a typical cosineas 20 = (cid:0) π (cid:1) (cid:82) π dφ (cid:82) ν EDS dν = 2 × (1 − ν EDS ) /π ,which returns ν EDS = 0 . ν EDS is lower than any angular scale at z > . ξ Lij is at its maxi-mum. On the other hand, when one of the two sourcesstands at redshift z < . ν EDS is always greater than ν ijth . In this regard, the approximation of ξ Lij with ξ Lijapp isexpected to work quite well when one of the source is atredshift closer than 0.5. For fainter sources, the fact thatwe deal with very narrow line-of-sight might show thelimit of our approximation. We remark that for largersky-coverage, just like LSST, our approximation is ex-pected to work well at any redshift. We argue this sincelarger angular openings tend to full-sky coverages and inthis limit case the lensing 2-point correlation function isexpected to be null and the peculiar sources located atthe transition scales become statistically less significant.We will quantify all these aspects in the next section.Finally, we conclude this section by underlining theimportance of the numerical approximation presented inEq. (51). The analytical estimation of ξ Lij ± with thepolynomial expansion of the growth function is numer-ically very easy to implement. Indeed, it requires only1-D numerical integration in k -space and nothing else.Hence, in spite of its simplicity, the fitting function inEq. (51) is powerful since its three parameters can beevaluated quickly for different cosmological parametersand then easily implemented in the Montecarlo analy-sis, where several runs over different cosmologies needto be accessible in short time. With the exact analyti-cal expression of ξ Lij , the multiple line-of-sight integralsare a huge obstacle in this regard. We finally under-line that this fitting method works well only for lensing2-point correlation function. This is due to the partic-ular features that lensing effect is very peaked for twonarrow line-of-sights and its angular average on the fullsky is zero. Neither of these properties simultaneouslyoccurs for any other effect here considered. However,this method remains powerful since lensing is anyway theleading effect among all.In the next section, we will provide some forecast forthe σ H for the case of the forthcoming surveys EDS andLSST. Finally, we will show the goodness of our fittingfunction in Eq. (51) in the estimation of σ H itself. V. A FIRST ESTIMATION FOR NEXTGENERATION SURVEYS
In this section, we want to apply the analytical andnumerical results previously obtained to the case of Su-perluminous Supernovae (SLSNe). In particular, we willconsider technical aspects for the forthcoming surveys ofEDS [1] and LSST [2]. In this regard, we will follow theexpected detection rate of SLSNe claimed in [28] whichwe report in Fig. 7. Having these histograms in mind,we generated two random surveys for the distribution ofSLSNe with the following specifics •
135 sources from EDS are generated in 7 redshift6 o f S L S N e EuclidLSST
FIG. 7. Simulated SLSNe distributions for Euclid Deep Sur-vey and LSST. From [28] bins where 0 . ≤ z ≤ . z = 0 .
5. For this survey, the an-gular distribution covers two line-of-sights at Northand South Poles, with angular opening of 20 deg per line-of-sight, •
929 sources from LSST are generated in 38 redshiftbins where 0 . ≤ z ≤ .
85 and the redshift binwidth is assumed to be ∆ z = 0 .
1. For this survey,the angular distribution spans a broad solid angleof 9000 deg .Results are summarized in Table I. Here we see that σ H /H EDS LSSTLensing 5 . × − . × − Doppler 2 . × − . × − Approximated Lensing 3 × − . × − TABLE I. Forecasts for the variance of H in EDS and LSST.In the first line, exact 2-point correlation function for lensingis considered. In the second line, there are the contributionsfrom 2-point correlation function of peculiar velocities. In thelast line, the expected error from our approximated 2-pointcorrelation function of lensing in Eq. (51) are shown. Thesevalues translate in the following values for the dispersion: forEDS, σ H /H = 0 .
002 for the exact estimation and σ H /H =0 .
002 for the approximated estimation. For LSST, σ H /H =0 . σ H /H = 0 . the dispersion associated to the measure of H , namely σ H ≡ (cid:113) σ H is of ∼ .
2% for EDS but its value dropsof almost 1 order of magnitude for LSST, where it con-tributes with a dispersion of ∼ . σ H = (cid:113) σ H L + σ H P V ≈ σ H L (cid:34) σ H P V σ H L + O (cid:32)(cid:18) σ H P V σ H L (cid:19) (cid:33)(cid:35) . (55)Again from Table I, we then get that the doppler effectcorrects the total cosmic variance associated to H by0 .
2% for LSST and by 0 .
02% for EDS.A further remark is about our approximation schemefor lensing proposed in Eq. (51). From Table I, we seethat the approximated method is in reasonable agree-ment with the exact evaluation done for LSST. The situ-ation gets worse for EDS. We address this behavior to thespecific sky coverage for the chosen surveys. Indeed, forEDS, the sky coverage is quite narrow and comparablefor the redshift of our interest with the threshold scalesestimated by ν ijth . This implies that the particular detailsof the angular dependence in ξ L ( ν ) are quite relevant.On the opposite case, LSST has a very broader angularopening. In this case, the sources are distributed with alarger angular opening and this renders the estimation of σ H L less sensible to the specific values of ν ijth .Finally, we comment on the limit ideal case of full-sky coverage . For a large number of sources, lensingnever contributes to the 2-point correlation function andthen the leading correction is entirely addressed to themonopole of ξ P V . This can be easily evaluated from thekernel in Eq. (41) and results are shown in Fig. 8. In this × - × - × - × - × - FIG. 8. Monopole of the doppler 2-point correlation func-tion. Thick black lines indicate where function is 0. Dashedlines stand for negative values whereas continues lines referto positive values. See [51] for a detailed discussion of this ideal case. σ H is entirely given by the doppler 2-pointcorrelation function. According to the redshift distribu-tion of LSST, we get that σ H /H = 5 . × − . A. Non-linear scales
So far we have taken into account only the linear powerspectrum in our analysis. However, non-linear physicscan affect significantly the amplitude of lensing (see, forinstance [52]). Hence, it is important to understand howmuch our results are robust when non-linear scales in thepower spectrum are taken into account. To this aim,we first adopt the approximation scheme outlined in Eq.(51). In this case, since non-linear scales can enhancelensing up to 1 order of magnitude, we might naivelyexpect that non-linear scales amplify our estimation for σ H by a factor √
10. However, this estimation is toomuch conservative.A more refined investigation about non-linear scalescan be done just by looking at the kernels W ± Lij inEqs. (50). Indeed, these are controlled by the functions j ( z ) /z and j ( z ) /z . Both these functions have theirmaxima in z = 0 and then oscillate around 0. Hence,they mostly contribute to the integrand when their ar-gument is O (1). For the antipodes-correlation, from Eq.(51), we get that the argument of the j n ’s is k ( x + y ).Hence, the most relevant part of the integrand occurswhen k ∼ ( x + y ) − . It follows then that small scalescontribute to ξ Lij − only at very low redshifts, when x + y approaches 0. On top of that, both j n ’s in W − Lij are mul-tiplied by power of xy . This means that a competitivebehavior occurs at low redshifts along the line-of-sightintegrations and small scales tend to be suppressed inthe ultimate evaluation. Thanks to this argument, wecan infer that the value of ξ Lij − is quite insensitive to thesmall scales and then ξ Lij − does not change dramaticallywhen non-linearities are taken into account.On the contrary, the j n ’s in W + Lij depend on k ( x − y ).This means that the integrand of ξ Lij + is heavily sourcedwhen x ≈ y . If we call d = x − y , we get that the j n ’s givesa non-negligible contribution to the total ξ Lij + whenever k ∼ d − . This implies that, along the line-of-sight inte-grations, the smaller the d , the higher the contributionfrom higher k . Hence, ξ Lij + is strongly dependent on thesmaller scales. This is in line with the fact that non-linearities are expected to heavily source lensing. Fromour analysis, it turns out that this is true only for ν = 1,namely when two sources are aligned.The fact that ξ Lij + is highly sensitive to non-linearscales whereas ξ Lij − is almost independent of them hasrelevant implications on the angular correlation scales ν ijth . Indeed, since the fact that lensing vanishes whenintegrated all around the observed sky is a pure geomet-rical effect, hence independent of the investigated scales,Eq. (54) is still a viable approximation. Then we getthat the contribution from non-linear scales constrains ν ijth to be closer to one. In fact, when | ξ Lij − | (cid:28) | ξ Lij + | , Eq.(54) can be expanded as ν ijth ≈ − | ξ Lij − || ξ Lij + | . (56)In this way, if the only effect of non-linearities is to in-crease | ξ Lij + | , for instance, by a factor C , the angularscales conversely tends to 1 with a factor C − . We thenhave two competitive effects: on one hand, the maxi-mum of the 2-point correlation function increases. On theother hand, the angular scales involved in this enhance-ment are less. Moreover, the product (1 − ν th ) | ξ Lij + | ≈ | ξ Lij − | is insensitive to the UV scales. As a consequence,due to non-linear effects in the power spectrum, onlysources located almost along the same line-of sight showa non-negligible correction.To test our arguments, we have then adopted the ap-proximation method outlined in Eqs. (51) and assumedthat ξ Lij − remains the same as for the linear spectrum,whereas the net effect of non-linearities is to amplify ξ Lij + by an overall factor 10. With this approach, we findthat σ H NL = 7 . × − H for EDS and σ H NL =6 . × − H for LSST. Hence, within the specific ofLSST σ H is indeed enhanced by almost one order ofmagnitude, whereas for EDS σ H increases by almost afactor 1.5.Despite of its rudeness, our analytical estimationagrees very well with the exact estimation where wehave chosen k UV = 10 h Mpc − and set a non-linearpower spectrum with the HaloFit model [53, 54]. Inthis case we obtain σ H NL = 1 . × − H for EDSand σ H NL = 7 . × − H for LSST. Moreover, wehave also applied the approximated formula in Eq. (51)with the non-linear values for ξ Lij + and ξ Lij − and we findthat σ H NL = 9 . × − H for EDS and σ H NL =7 . × − H for LSST. We notice that all our estima-tions agree quite well for both surveys.Our results then show that, for EDS, even if non-linearscales are expected to enhance the lensing correction byalmost one order of magnitude, the intrinsic error associ-ated to the measurement of H is almost insensitive to thenon-linear scales, since it becomes σ H NL /H = 0 . H within the specific of LSST,raising σ H to the value σ H NL /H = 0 . VI. SUMMARY AND CONCLUSIONS
In this work, we have studied the impact of cosmolog-ical inhomogeneities on the estimation of H from thehigh redshift Hubble diagram. Our analysis considersthe possibility, discussed in [28], that a statistically rele-vant number of Superluminous Supernovae could be de-tected in the next years by EDS [1] and LSST [2]. Inthis regards, less conservative studies about the Hubblediagram at high redshifts ( z ≤ .
5) have been also in-vestigated by exploiting exact inhomogeneous models ingeneral relativity [55–58].On the contrary, along our study, we have adopteda conservative approach based on linear perturbationswithin the Cosmic Concordance model. In this frame-work, first of all we have derived the 2-point correlationfunction of luminosity distance-redshift relation and pro-vided explicit fully relativistic analytic expressions forits angular spectra. Our derivations agree with thosealready obtained in literature [40], modulo a differentclassification of involved terms. It turns out in our anal-ysis that lensing is the leading effect at the consideredredshift, as one may expect.In particular, as shown in Fig. 4, we get that theangular multipoles expansion of the lensing 2-point cor-relation function at high redshift converges quite slowly.Moreover, lensing dipole rapidly becomes negligible withrespect to the other multipoles. This allows to safely ne-glect its contribution to the 2-point correlation function,since it is expected to be contaminated by the observer’speculiar motion.Furthermore, the role of higher multipoles in the lens-ing spectra plays a relevant role for partial sky-coveragesurveys. Indeed, in the ideal case of large number ofsources distributed all over the sky, lensing does not con-tribute to the total cosmic variance since its monopoleis null. We have shown this also for cross-correlations oflensing with the other relativistic effects. In this idealcase, then, the leading correction to the estimation of H is due to the peculiar velocities of the sources. However,realistic surveys deal with limited sky coverage and thismakes lensing contribution no longer vanishing. In fact, according to the specific of EDS and LSST and to whathas been claimed in [28], we forecast that the intrinsic er-ror from cosmic variance associated to H is of ∼ .
03 %for LSST and 0 . z = 0 .
1) have beenconsidered. Here we extend the analysis of [27] sincethere only the peculiar motion of the sources is takeninto account. In fact, this is the leading correction ex-pected at those redshifts [52]. The interesting result isthat low redshift surveys discussed in [27] admits a cos-mic variance for H of ∼ H which tendsto decrease when higher redshift sources are considered.Finally, we remark that our results are not able toalleviate the tension between local and distant measure-ments of the Hubble constant. However, they indicatethat the analysis of fainter sources does not increase thetheoretical uncertainty on H . The price to pay standsin the fact that the Hubble diagram at higher redshift isno longer model independent. ACKNOWLEDGMENTS
The authors are thankful to Vincenzo Cardone, EneaDi Dio, Ruth Durrer and Kazuya Koyama for use-ful discussions. GF acknowledges support by FCTunder the program “Stimulus” with the grant no.CEECIND/04399/2017/CP1387/CT0026. BF is sup-ported by the PhD program of the University ofPortsmouth. GM is supported in part by INFN underthe program TAsP (
Theoretical Astroparticle Physics ). Appendix A: Kernels derivation
In this appendix, we report some technical aspects of the derivation of Eqs. (27) and (28). We start with thederivation of ξ E in Eqs. (25). There are 5 terms and the geometrical sets for the integration are shown in Fig. 9respectively for local effects (left panel), such as PV and SW, and integrated terms (right panel), like TD, L and ISW. a. Peculiar Velocity As illustrated in Fig. 9, we define x = η o − η i and y = η o − η j . Combining Eqs. (14) and(22), for peculiar velocity operator we obtain W P V ij = Ξ i Ξ j G i G j ∂ x ∂ y (cid:90) d Ω k π e i k · ( x − y ) , (A1)9 O θ ~x~y R ~v i ~v j j i O θ ~x~y Rj i FIG. 9. Geometrical relations between line-of sight variables for local effects (left), i.e. PV and SW, and integrated ones (right),namely L, TD and ISW. where we have used Eqs. (26). Since k has rotational freedom, we align its component along the azimutal axis with x − y , so that the integral over the solid angle becomes (cid:90) d Ω k e i k · ( x − y ) = 2 π (cid:90) − d (cos θ k ) e ikR cos θ k = 4 πj ( kR ) , (A2)where j is the 0-th order spherical Bessel function and R = R ( x, y, ν ) is taken as defined in Eqs. (35). Thus we areleft with the evaluation of the derivatives of the spherical Bessel function ∂ x ∂ y j ( kR ) = k ∂ x [ j (cid:48) ( kR ) ∂ y R ] = k j (cid:48)(cid:48) ( kR ) ∂ x R ∂ y R + k j (cid:48) ( kR ) ∂ x ∂ y R , (A3)where j (cid:48) n ( z ) = ∂ z j n ( z ), j (cid:48)(cid:48) n ( z ) = ∂ z j n ( z ) and so on. Finally, by exploiting the properties of the j n reported in Eq.(C3), we obtain W P V ij = Ξ i Ξ j G i G j k (cid:26) xy (1 − ν ) R j ( kR ) + ν j ( kR ) − j ( kR )] (cid:27) , (A4)in agreement with [27]. b. Lensing From right panel of Fig. 9, we define x = η o − η and y = η o − η (cid:48) . From the definition of the lensingoperator in Eqs. (14), we then obtain W Lij = 1∆ η i ∆ η j (cid:90) η o η i dη η − η i η o − η (cid:90) η o η j dη (cid:48) η (cid:48) − η j η o − η (cid:48) g ( η ) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π ∆ x ∆ y e i k · ( x − y ) . (A5)We recall now that ∆ r = ∂ θ + cot θ ∂ θ + 1 / sin θ ∂ φ is the angular Laplacian which can be either evaluated along the r = x or r = y direction. It can be written in terms of the 3-dimensional one as∆ r = r (cid:16) (cid:126) ∇ · (cid:126) ∇ − ∇ r (cid:17) = r (cid:18) (cid:126) ∇ · (cid:126) ∇ − ∂ r − r ∂ r (cid:19) . (A6)We use Eq. (A6) to replace the Laplacian in the angular coordinates and get easier evaluable derivatives. With sucha substitution, the last integral in Eq. (A5) becomes (cid:90) d Ω k π ∆ ∆ (cid:48) e i k · ( x − y ) = (cid:90) d Ω k π x (cid:18) (cid:126) ∇ x · (cid:126) ∇ x − ∂ x − x ∂ x (cid:19) y (cid:18) (cid:126) ∇ y · (cid:126) ∇ y − ∂ y − y ∂ y (cid:19) e i k · ( x − y ) = (cid:90) d Ω k π x (cid:18) − k − ∂ x − x ∂ x (cid:19) y (cid:18) − k − ∂ y − y ∂ y (cid:19) e i k · ( x − y ) = x (cid:18) − k − ∂ x − x ∂ x (cid:19) y (cid:18) − k − ∂ y − y ∂ y (cid:19) j ( kR ) , (A7)0where last equality has been obtained thanks to Eq. (A2). Hence, by using the recursive relations of the j n ’s in Eq.(C3), after a bit of algebra, W Lij can be rewritten as W Lij = 1∆ η i η j (cid:90) η o η i dη η − η i η o − η (cid:90) η o η j dη (cid:48) η (cid:48) − η j η o − η (cid:48) g ( η ) g ( η (cid:48) ) g ( η o ) (cid:2) k H j ( kR ) − k H Lj ( kR ) + k (cid:0) L − H (cid:1) j ( kR ) + 4 kLj ( kR ) (cid:3) (A8)where H , L and R are defined in Eqs. (35). c. Sachs-Wolfe In analogy with what has been done for the peculiar velocity, we refer to left panel of Fig. 9,we define x = η o − η i and y = η o − η j and consider the operator ˆ O SW in Eqs. (14). Since no spatial derivatives areconsidered, we simply get, though Eq. (A2) W SW ij =(1 + Ξ i )(1 + Ξ j ) g ( η i ) g ( η o ) g ( η j ) g ( η o ) (cid:90) d Ω k π e i k · ( x − y ) =(1 + Ξ i )(1 + Ξ j ) g ( η i ) g ( η o ) g ( η j ) g ( η o ) j ( kR ) , (A9)where R = R ( η i , η j , ν ) is taken from Eqs. (26). d. Integrated Sachs-Wolfe In analogy with the geometrical set used for lensing (right panel in Fig. 9), we define x = η o − η and y = χ (cid:48) = η o − η (cid:48) . Given that we simply obtain W ISW ij =4 Ξ i Ξ j (cid:90) η o η i dη (cid:90) η o η j dη (cid:48) ∂ η g ( η ) g ( η o ) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π e i k · ( x − y ) =4Ξ i Ξ j (cid:90) η o η i dη (cid:90) η o η j dη (cid:48) ∂ η g ( η ) g ( η o ) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) j ( kR ) , (A10)where we have used again Eq. (A2). We stress that now R = R ( η, η (cid:48) , ν ) from Eqs. (26) is integrated along the twoline-of-sights, differently from the case of PV and SW. e. Time Delay TD term can be evaluated in complete analogy with ISW. The only difference is that the growthfunctions are not derived. Following then the evaluation for ISW, we get W T Dij = 4∆ η i ∆ η j (cid:90) η o η i dη (cid:90) η o η j dη (cid:48) g ( η ) g ( η o ) g ( η (cid:48) ) g ( η o ) j ( kR ) . (A11)There are then 10 mixed terms that come out of the sum in Eq. (20) when E (cid:54) = E (cid:48) . As done for the pure terms,we will use the auxiliary variables x, y , whose definitions will be given for each term and differ whether the effect islocal or integrated along the line-of-sight. To our aim, we first report the following useful expressions (cid:90) d Ω k π ∂ x e i k · ( x − y ) = k νy − xR j ( kR ) , (cid:90) d Ω k π ∆ x e i k · ( x − y ) = k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) , (cid:90) d Ω k π ∂ x ∆ y e i k · ( x − y ) = − k y ( x − νy )( y − νx ) R j ( kR ) + k y (cid:20) y ( x − νy ) R − ν (cid:21) j ( kR ) − kyR (cid:2) k y ( x − νy ) − ν (cid:3) j ( kR ) . (A12)where R , L and H are taken from Eqs. (35).1 f. Peculiar Velocity - Lensing From the definition of the operators in Eqs. (14), we define x = η o − η i and y = η o − η (cid:48) , since PV is a local effect whereas L is integrated along the line-of-sight. With this, we get W P V i,Lj = Ξ i ∆ η j (cid:90) η (0) i η in dη a ( η ) a ( η o ) g ( η ) g ( η o ) (cid:90) η o η (0) j dη (cid:48) η (cid:48) − η (0) j η o − η (cid:48) (cid:90) d Ω k π g ( η (cid:48) ) g ( η o ) ∂ x ∆ y e i k · ( x − y ) . (A13)Hence, by using Eqs. (A12), we obtain W P V i,Lj = Ξ i ∆ η j G i (cid:90) η o η (0) j dη (cid:48) η (cid:48) − η (0) j η o − η (cid:48) (cid:26) − k y ( x − νy )( y − νx ) R j + k y (cid:20) y ( x − νy ) R − ν (cid:21) j − kyR (cid:2) k y ( x − νy ) − ν (cid:3) j (cid:27) , (A14)where we have omitted the dependence on kR ( x, y, ν ) in the j n ’s. g. Peculiar Velocity - Sachs-Wolfe Since both effects are local effects, we define x = η o − η i and y = η o − η j .From the definitions of the respective operators in Eqs. (14), we find W P V i,SW j =Ξ i (1 + Ξ j ) (cid:90) η i η in dη a ( η ) a ( η o ) g ( η ) g ( η o ) g ( η j ) g ( η o ) (cid:90) d Ω k π ∂ x e i k · ( x − y ) =Ξ i (1 + Ξ j ) G i g ( η j ) g ( η o ) k νy − xR j ( kR ) , (A15)where we made use again of Eqs. (A12). h. Peculiar Velocity - Integrated Sachs-Wolfe Since PV is a local effect whereas ISW is integrated along theline-of-sight, we define x = η o − η i and y = η o − η (cid:48) . Hence, from the definitions of the operators in Eqs. (14), we get W P V i,ISW j =2 Ξ i Ξ j (cid:90) η i η in dη a ( η ) a ( η o ) g ( η ) g ( η o ) (cid:90) η o η (0) j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π ∂ x e i k · ( x − y ) =2 Ξ i Ξ j G i (cid:90) η o η (0) j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) k νy − xR j ( kR ) , (A16)where we have used as before Eqs. (A12). i. Peculiar Velocity - Time Delay Since PV is a local effect whereas TD is integrated along the line-of-sight, wedefine x = η o − η i and y = η o − η (cid:48) . Hence, from the definitions of the operators in Eqs. (14), we obtain W P V i,T Dj = − i ∆ η j (cid:90) η i η in dη a ( η ) a ( η o ) g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π ∂ x e i k · ( x − y ) = − i ∆ η j G i (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) k νy − xR j ( kR ) (A17)where we have used as before Eqs. (A12). j. Lensing - Sachs-Wolfe Since L is an integrated effect and SW is a local one, we set the geometrical variablesas x = η o − η and y = η o − η j . Hence, the definition of the operators in Eqs. (14) leads to W Li,SW j = 1 + Ξ j ∆ η i g ( η j ) g ( η o ) (cid:90) η o η i dη η − η j η o − η g ( η ) g ( η o ) (cid:90) d Ω k π ∆ y e i k · ( x − y ) = 1 + Ξ j ∆ η i g ( η j ) g ( η o ) (cid:90) η o η i dη η − η j η o − η g ( η ) g ( η o ) (cid:20) k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) (cid:21) , (A18)2thanks to the Eqs. (A12). k. Lensing - Integral Sachs Wolfe Since both effects are integrated along the line-of-sights, we define x = η o − η and y = η o − η (cid:48) . Hence, the action of the operators (14), combined with Eqs. (A12), gives W Li,ISW j =2 Ξ j ∆ η i (cid:90) η o η i dη η − η j η o − η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π ∆ y e i k · ( x − y ) =2 Ξ j ∆ η i (cid:90) η o η i dη η − η j η o − η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) (cid:20) k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) (cid:21) . (A19) l. Lensing - Time Delay Since L and TD are both integrated along the line-of-sights, we define x = η o − η and y = η o − η (cid:48) , so, with the definition (14), we simply evaluate W Li,T Dj = − η i ∆ η j (cid:90) η o η i dη η − η j η o − η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π ∆ y e i k · ( x − y ) = − η i ∆ η i (cid:90) η o η i dη η − η j η o − η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) (cid:20) k H j ( kR ) − k (cid:18) H R − L (cid:19) j ( kR ) (cid:21) , (A20)thanks to the relations (A12). m. Time Delay - Integrated Sachs-Wolfe Just as for the previous term, we define x = η o − η and y = η o − η (cid:48) since both TD and ISW are integrated effects. Hence, from the operators defined in Eqs. (14), we obtain W T Di,ISW j =2 (1 + Ξ i )Ξ j g ( η i ) g ( η o ) (cid:90) η o η j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π e i k · ( x − y ) =2 (1 + Ξ i )Ξ j g ( η i ) g ( η o ) (cid:90) η o η j dη (cid:48) ∂ η (cid:48) g ( η (cid:48) ) g ( η o ) j ( kR ) , (A21)where we have applied Eq. (A2). n. Sachs-Wolfe - Time Delay SW and TD are respectively local and integrated effects. Because of that, wedefine x = η o − η i and y = η o − η (cid:48) and then, from the operators (14), we obtain W SW i,T Dj = − i ∆ η j g ( η i ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π e i k · ( x − y ) = − i ∆ η j g ( η i ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) j ( kR ) , (A22)where we have applied Eq. (A2). o. Integrated Sachs-Wolfe - Time Delay ISW and TD are both integrated effects. Hence, we define x = η o − η and y = η o − η (cid:48) and then, from the operators (14), we obtain W ISW i,T Dj = − i ∆ η j (cid:90) η o η (0) i dη ∂ η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) (cid:90) d Ω k π e i k · ( x − y ) = − i ∆ η j (cid:90) η o η i dη ∂ η g ( η ) g ( η o ) (cid:90) η o η j dη (cid:48) g ( η (cid:48) ) g ( η o ) j ( kR ) , (A23)where we have applied Eq. (A2).3 Appendix B: Useful Expressions
Here we report the explicit value for the Q (cid:96)n in Eq. (40) for (cid:96) = 1 , ,
3, which allows to write analytically the kernelfor the dipole of the lensing 2-point correlation function in Eq. (39). For (cid:96) = 1, we obtain Q ( z ) = − z SinInt( zk ) − sin( zk )2 k z − k cos( zk ) + 1 k ,Q ( z ) = 1 k SinInt( zk ) − zk ,Q ( z ) = − zk SinInt( zk ) − zk ) k + 2 k ,Q ( z ) = 3 sin( zk ) k − zk − z cos( zk ) k ,Q ( z ) = − z cos( zk ) k + 5 z sin( zk ) k + 8 cos( zk ) k − k ,Q ( z ) = 7 z sin( zk ) k − z cos( zk ) k + 8 zk −
30 sin( zk ) k + 22 z cos( zk ) k . (B1)For (cid:96) = 2, we have Q ( z ) = − sin( kz )2 k z + cos( kz )2 k z + SinInt( kz )2 k − z ,Q ( z ) = − kz )2 k z − cos( kz )2 k + 2 k − z SinInt( kz )2 k ,Q ( z ) = 3SinInt( kz ) k − zk − sin( kz ) k ,Q ( z ) = − z SinInt( kz ) k − z sin( kz ) k − kz ) k + 8 k ,Q ( z ) = − z sin( kz ) k − zk + 15 sin( kz ) k − z cos( kz ) k ,Q ( z ) = − z sin( kz ) k − z cos( kz ) k + 33 z sin( kz ) k + 48 cos( kz ) k − k . (B2)Finally, for (cid:96) = 3 we get Q ( z ) = − kz )4 k z + 5 cos( kz )4 k z − sin( kz )8 k z − z SinInt( kz )8 − cos( kz )8 k + 23 k ,Q ( z ) = 3SinInt( kz )2 k − kz )2 k z + 5 cos( kz )2 k z − z k ,Q ( z ) = − z SinInt( kz )2 k −
15 sin( kz )2 k z − cos( kz )2 k + 8 k ,Q ( z ) = 15SinInt( kz ) k − zk − kz ) k + z cos( kz ) k ,Q ( z ) = − z SinInt( kz ) k + z cos( kz ) k − z sin( kz ) k −
48 cos( kz ) k + 48 k ,Q ( z ) = − z sin( kz ) k + z cos( kz ) k − zk + 105 sin( kz ) k − z cos( kz ) k . (B3) Appendix C: Spherical Bessel functions
Spherical Bessel functions j n ( z ) are generated for any integer n through the recursive formula j n ( z ) = ( − n z n (cid:18) z ddz (cid:19) n (cid:18) sin zz (cid:19) , (C1)and are the fundamental solutions of the ODE z d dz j n ( z ) + 2 z ddz j n ( z ) + (cid:2) z − n ( n + 1) (cid:3) j n ( z ) = 0 . (C2)4The useful relations that we report for our purposes are the following j n − ( z ) = 2 n + 1 z j n ( z ) − j n +1 ( z ) ,j (cid:48) n ( z ) = n n + 1 j n − ( z ) − n + 12 n + 1 j n +1 ( z ) ,j (cid:48) n ( z ) = j n − ( z ) − n + 1 z j n ( z , ) j (cid:48) n ( z ) = − j n +1 ( z ) + nz j n ( z ) , (C3)where, in particular, the first equation allows to extend the definition of the j n ( z ) also when n is negative, whereasthe remaining three equations allows to express also their derivatives in terms of the j n ’s themselves. Appendix D: Legendre polynomials and multipole expansion
Let be f (cos θ ) a function which does not depend on the azimutal angle φ . In terms of ν ≡ n x · n y = cos θ , we canwrite the angular laplacian ∆ as ∆ f ( ν ) = ∂ ν (cid:2)(cid:0) − ν (cid:1) ∂ ν f ( ν ) (cid:3) . (D1)In this way, the multipole expansion of ∆ f in terms of the Legendre polynomials P (cid:96) ( ν ) can be integrated by parttwice. This leads to (cid:90) − dν P (cid:96) ( ν )∆ f ( ν ) = (cid:90) − dν P (cid:96) ( ν ) ∂ ν (cid:2)(cid:0) − ν (cid:1) ∂ ν f ( ν ) (cid:3) = − (cid:90) − dν ∂ ν P (cid:96) ( ν ) (cid:2)(cid:0) − ν (cid:1) ∂ ν f ( ν ) (cid:3) = (cid:90) − dν ∆ P (cid:96) ( ν ) f ( ν ) = − (cid:96) ( (cid:96) + 1) (cid:90) − dν P (cid:96) ( ν ) f ( ν ) , (D2)where in the last line we have used the fact that P (cid:96) ( ν ) are eigenfunctions of ∆ . This simple result enlightens thecase of lensing multipoles in Eq. (39). Indeed, those coefficients can be formally written in the following form (cid:90) − dν P (cid:96) ( ν )∆ x ∆ y f ( R ( ν )) . (D3)In fact, since (D3) is invariant under rotation of the coordinate system in the plane spanned by x and y , we can applyEq. (D2) as follows: we first rotate the coordinate system in the integral in order to align them with the direction x ,hence we apply Eq. (D2). After that, we align again the coordinate system with y and finally apply Eq. (D2) again.Since Eq. (D2) is applied twice, we get the prefactor (cid:96) ( (cid:96) + 1) . [1] R. Laureijs et al. (EUCLID), (2011), arXiv:1110.3193[astro-ph.CO].[2] P. A. Abell et al. (LSST Science, LSST Project), (2009),arXiv:0912.0201 [astro-ph.IM].[3] G. Mellema et al., Exper. Astron. , 235 (2013),arXiv:1210.0197 [astro-ph.CO].[4] D. R. DeBoer et al., Publ. Astron. Soc. Pac. , 045001(2017), arXiv:1606.07473 [astro-ph.IM].[5] G. Marozzi, G. Fanizza, E. Di Dio, and R. Durrer, JCAP , 028 (2016), arXiv:1605.08761 [astro-ph.CO].[6] G. Marozzi, G. Fanizza, E. Di Dio, and R. Durrer, Phys. Rev. Lett. , 211301 (2017), arXiv:1612.07650 [astro-ph.CO].[7] G. Marozzi, G. Fanizza, E. Di Dio, and R. Durrer,Phys. Rev. D , 023535 (2018), arXiv:1612.07263 [astro-ph.CO].[8] E. Di Dio, R. Durrer, G. Fanizza, and G. Marozzi, Phys.Rev. D , 043508 (2019), arXiv:1905.12573 [astro-ph.CO].[9] N. Aghanim et al. (Planck), Astron. Astrophys. , A6(2020), arXiv:1807.06209 [astro-ph.CO].[10] A. G. Riess, S. Casertano, W. Yuan, J. B. Bow- ers, L. Macri, J. C. Zinn, and D. Scolnic, (2020),arXiv:2012.08534 [astro-ph.CO].[11] W. L. Freedman, B. F. Madore, D. Hatt, T. J. Hoyt,I. S. Jang, R. L. Beaton, C. R. Burns, M. G. Lee, A. J.Monson, J. R. Neeley, and et al., The AstrophysicalJournal , 34 (2019).[12] W. L. Freedman, B. F. Madore, T. Hoyt, I. S. Jang,R. Beaton, M. G. Lee, A. Monson, J. Neeley, and J. Rich,The Astrophysical Journal , 57 (2020).[13] K. C. Wong, S. H. Suyu, G. C.-F. Chen, C. E. Rusu,M. Millon, D. Sluse, V. Bonvin, C. D. Fassnacht,S. Taubenberger, M. W. Auger, and et al., Monthly No-tices of the Royal Astronomical Society , 1420?1439(2019).[14] W. Yuan, A. G. Riess, L. M. Macri, S. Casertano, andD. M. Scolnic, The Astrophysical Journal , 61 (2019).[15] D. W. Pesce, J. A. Braatz, M. J. Reid, A. G. Riess,D. Scolnic, J. J. Condon, F. Gao, C. Henkel, C. M. V.Impellizzeri, C. Y. Kuo, and et al., The AstrophysicalJournal , L1 (2020).[16] L. Verde, T. Treu, and A. G. Riess, Nature Astronomy , 891 (2019).[17] J. L. Bernal, L. Verde, and A. G. Riess, JCAP ,019 (2016), arXiv:1607.05617 [astro-ph.CO].[18] K. Jedamzik, L. Pogosian, and G.-B. Zhao, (2020),arXiv:2010.04158 [astro-ph.CO].[19] W. Beenakker and D. Venhoek, (2021), arXiv:2101.01372[astro-ph.CO].[20] L. Knox and M. Millea, Phys. Rev. D , 043533 (2020),arXiv:1908.03663 [astro-ph.CO].[21] E. Di Valentino et al., (2020), arXiv:2008.11284 [astro-ph.CO].[22] E. Di Valentino, A. Melchiorri, and J. Silk, Phys. Lett.B , 242 (2016), arXiv:1606.00634 [astro-ph.CO].[23] E. Di Valentino, A. Melchiorri, and O. Mena, Phys. Rev.D , 043503 (2017), arXiv:1704.08342 [astro-ph.CO].[24] S. Vagnozzi, Phys. Rev. D , 023518 (2020),arXiv:1907.07569 [astro-ph.CO].[25] M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A. A.Starobinsky, and C. Umilt`a, JCAP , 044 (2020),arXiv:2004.14349 [astro-ph.CO].[26] E. J. Baxter and B. D. Sherwin, Monthly Notices of theRoyal Astronomical Society , 1823?1835 (2020).[27] I. Ben-Dayan, R. Durrer, G. Marozzi, and D. J. Schwarz,Phys. Rev. Lett. , 221301 (2014), arXiv:1401.7973[astro-ph.CO].[28] C. Inserra et al. (DES), (2020), arXiv:2004.12218 [astro-ph.CO].[29] K. Abazajian et al., (2019), arXiv:1907.04473 [astro-ph.IM].[30] M. Gasperini, G. Marozzi, F. Nugier, and G. Veneziano,JCAP , 008 (2011), arXiv:1104.1167 [astro-ph.CO].[31] C. Bonvin, C. Clarkson, R. Durrer, R. Maartens, andO. Umeh, JCAP , 040 (2015), arXiv:1504.01676 [astro-ph.CO].[32] A. Heinesen, P. Mourier, and T. Buchert, Class. Quant.Grav. , 075001 (2019), arXiv:1811.01374 [gr-qc].[33] G. Fanizza, M. Gasperini, G. Marozzi, and G. Veneziano,JCAP , 017 (2020), arXiv:1911.09469 [gr-qc]. [34] I. Ben-Dayan, M. Gasperini, G. Marozzi, F. Nugier, andG. Veneziano, JCAP , 036 (2012), arXiv:1202.1247[astro-ph.CO].[35] P. Fleury, C. Clarkson, and R. Maartens, JCAP , 062(2017), arXiv:1612.03726 [astro-ph.CO].[36] J. Yoo and R. Durrer, JCAP , 016 (2017),arXiv:1705.05839 [astro-ph.CO].[37] G. Fanizza, M. Gasperini, G. Marozzi, and G. Veneziano,Phys. Lett. B , 505 (2016), arXiv:1512.08489 [astro-ph.CO].[38] G. Marozzi, Class. Quant. Grav. , 045004 (2015),[Erratum: Class.Quant.Grav. 32, 179501 (2015)],arXiv:1406.1135 [astro-ph.CO].[39] G. Fanizza, J. Yoo, and S. G. Biern, JCAP , 037(2018), arXiv:1805.05959 [gr-qc].[40] C. Bonvin, R. Durrer, and M. Gasparini, Phys. Rev.D , 023523 (2006), [Erratum: Phys.Rev.D 85, 029901(2012)], arXiv:astro-ph/0511183.[41] I. Ben-Dayan, G. Marozzi, F. Nugier, and G. Veneziano,JCAP , 045 (2012), arXiv:1209.4326 [astro-ph.CO].[42] G. Fanizza, M. Gasperini, G. Marozzi, and G. Veneziano,JCAP , 020 (2015), arXiv:1506.02003 [astro-ph.CO].[43] O. Umeh, C. Clarkson, and R. Maartens, Class. Quant.Grav. , 205001 (2014), arXiv:1402.1933 [astro-ph.CO].[44] P. Peter and J.-P. Uzan, Primordial Cosmology, OxfordGraduate Texts (Oxford University Press, 2013).[45] B. Fiorini, Master Thesis (2018).[46] F. Scaccabarozzi, J. Yoo, and S. G. Biern, JCAP ,024 (2018), arXiv:1807.09796 [astro-ph.CO].[47] V. Tansella, G. Jelic-Cizmek, C. Bonvin, and R. Durrer,JCAP , 032 (2018), arXiv:1806.11090 [astro-ph.CO].[48] E. Di Dio and U. Seljak, JCAP , 050 (2019),arXiv:1811.03054 [astro-ph.CO].[49] F. Beutler and E. Di Dio, JCAP , 048 (2020),arXiv:2004.08014 [astro-ph.CO].[50] D. J. Eisenstein and W. Hu, Astrophys. J. , 605(1998), arXiv:astro-ph/9709112.[51] J. Yoo, Phys. Rev. D , 043507 (2020),arXiv:1911.07869 [astro-ph.CO].[52] I. Ben-Dayan, M. Gasperini, G. Marozzi, F. Nugier, andG. Veneziano, JCAP , 002 (2013), arXiv:1302.0740[astro-ph.CO].[53] R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White,C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou,and H. M. P. Couchman, Monthly Notices of the RoyalAstronomical Society , 1311?1332 (2003).[54] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, andM. Oguri, The Astrophysical Journal , 152 (2012).[55] L. Cosmai, G. Fanizza, M. Gasperini, and L. Tedesco,Class. Quant. Grav. , 095011 (2013), arXiv:1303.5484[gr-qc].[56] A. E. Romano, Int. J. Mod. Phys. D , 1850102 (2018),arXiv:1609.04081 [astro-ph.CO].[57] L. Cosmai, G. Fanizza, F. Sylos Labini, L. Pietronero,and L. Tedesco, Class. Quant. Grav. , 045007 (2019),arXiv:1810.06318 [astro-ph.CO].[58] S. A. Vallejo-Pe˜na and A. E. Romano, JCAP03