Cross-correlating radial peculiar velocities and CMB lensing convergence

We study, for the first time, the cross correlation between the angular distribution of radial peculiar velocities (PV) and the lensing convergence of cosmic microwave background (CMB) photons. We derive theoretical expectations for the signal and its covariance and assess its detectability with existing and forthcoming surveys. We find that such cross-correlations are expected to improve constraints on different gravitational models by partially breaking degeneracies with the matter density. We identify in the distance-scaling dispersion of the peculiar velocities the most relevant source of noise in the cross correlation. For this reason, we also study how the above picture changes assuming a redshift-independent scatter for the PV, obtained for example using a reconstruction technique. Our results show that the cross correlation might be detected in the near future combining PV measurements from DESI and the convergence map from CMB-S4. Using realistic direct PV measurements we predict a cumulative signal-to-noise ratio of approximately $3.8 \sigma$ using data on angular scales $3 \leq \ell \leq 200$. For an idealized reconstructed peculiar velocity map extending up to redshift $z=0.15$ and a smoothing scale of $4$ Mpc $h^{-1}$ we predict a cumulative signal-to-noise ratio of approximately $ 27 \sigma$ from angular scales $3 \leq \ell \leq200 $. We conclude that currently reconstructed peculiar velocities have more constraining power than directly observed ones, even though they are more cosmological-model dependent.


Introduction
Nowadays, at least on cosmological scales, the most commonly accepted description for the evolution of the Universe is the ΛCDM model [1][2][3][4], which consistently accounts for a smooth transition between three distinct stages of evolution. More precisely, it describes an expanding Universe which is initially filled mostly by radiation, i.e. photons and neutrinos. The latter are diluted to the point that their density becomes comparable to or smaller than the Cold Dark Matter (CDM) one. At this stage the photons are initially confined by the presence of baryons in a dense plasma, but eventually become free once that the Universe expanded enough to make the Thomson scattering inefficient. The last photons that scattered off the electrons constitute the Cosmic Microwave Background (CMB). The initially small density and temperature fluctuations of matter and radiation are correlated, and trigger the gravitational collapse of matter. Hotter regions become denser and and colder ones emptier. Eventually, this led to the formation of a complex web of CDM structures where baryons are trapped, and undergo a similar process resulting in the formation of stars and galaxies. When the cosmic expansion dilutes the CDM density enough, the most abundant component of the Universe becomes the Cosmological Constant Λ. Maybe coincidentally [5,6], this happens about half the age of the universe ago. As the name itself suggests, Λ does not dilute and we enter the final stage of the evolution of the Universe.
From an observational point of view, it is relatively simple to observe the leading character of the first act in the ΛCDM expansion history by looking at the CMB. These photons travelled from their last Thomson scattering until today, and their distribution has been probed by space-based experiments like WMAP and Planck [1,7], and ground-based like the South Pole Telescope (SPT) [8] and the Atacama Cosmology Telescope (ACT) [9]. On the other hand, it is impossible to probe with direct observations the distribution of CDM and Λ since they only interact through gravity. However, baryonic structures are formed within larger CDM haloes [10], and hence trace the underlying CDM distribution. The CMB photons are also sensitive to the matter distribution because of the gravitational lensing induced on them by the Large Scale Structure (LSS), hence the CMB itself is also a CDM tracer. This, in turn, suggests that the signals of the CMB and of any tracer of the CDM distribution should be correlated. This hypothesis was tested (and confirmed) using a variety of CDM tracers such as the spatial distributions of galaxies, Type Ia supernovae, and quasars, see for example Refs. [11][12][13][14][15][16].
It is important to stress that these tracers are non trivially related to the underlying total matter distribution δ = (ρ −ρ)/ρ, where ρ is the total matter density and ρ its mean. Indeed, one needs to introduce auxiliary parameters not predicted by the ΛCDM model itself, like the galaxy bias b relating the CDM density to the number density of galaxies. These auxiliary parameters and the cosmological ones may be (and often are) degenerate in the actual measurement, leading to potential model dependent biases in the cosmological inference. Luckily, the self-correlation and the cross correlations of different tracers may break these degeneracies, making this subject an active field of research [17,18]. Furthermore, over the last decade, high precision cosmology measurements have challenged our understanding of the growth and distribution of LSS due to an established ∼ 3σ tension on the observed clustering rate at the pivot scale of 8 Mpc between early and late Universe cosmological probes, usually referred to as S 8 tension [19][20][21]. This tension, interestingly, is closely related to the one on the value of the Hubble factor today H 0 [22][23][24], as it is difficult to alleviate either of the two without worsening the other [25][26][27]. These challenges to the standard model paradigm advocate for new independent probes of the LSS distribution.
With the above motivations in mind, in this work we investigate for the first time the cross-correlation between the CMB lensing signal and the distribution of radial peculiar velocities (PV) of galaxies along the line of sight. Indeed, due to the presence of LSS, galaxies are not at rest with respect to the cosmological rest frame and move from underdense regions to overdense ones, following the strength of the local gravitational potential. We define the radial PV as the projection along the line of sight of the velocity induced by these inhomogeneities. In recent years, several works have shown [28][29][30][31][32][33][34][35][36] the positive impact of PV measurements in combination with redshift survey data in constraining cosmological parameters. Ongoing and upcoming surveys like DESI [37] and 4HS [38] will provide an order of magnitude more PV measurements than existing catalogs like SDSS [39], with increased precision and redshift depth. For these reasons, cosmological inference with PV is an active field of research and a promising avenue to test our understanding of gravity.
At linear order in perturbation theory, the gradient of the PV is proportional to the matter density contrast δ, making them a powerful tracer of the CDM distribution with no dependence on the galaxy bias b(z). State of art PV surveys usually target objects located within redshift z ≤ 0.1. On the other hand, most of the gravitational lensing is induced by the LSS at higher redshift, z ∼ 1. Because of this huge difference, one would naively expect the signal from the cross correlation of the two to be essentially vanishing. Surprisingly, as we are going to show, this is not the case. Even if intrinsically small, our results indicate that current and forthcoming experiments may be able to detect the aforementioned cross-correlation. Our findings show that the main obstacle towards the detection of the signal is the scale-dependent dispersion associated with direct PV measurements. We also find that reconstructed PV can alleviate this problem, but require model dependent assumptions on cosmological parameters such as the linear galaxy bias b(z).
The structure of the paper is the following: in Sec. 2 we briefly review the formalism used to define angular cross-correlations between cosmological observables. In Sec. 3 we compute the theoretical expectation for the cross-correlation between the CMB convergence κ, a measurement of the CMB lensing, and radial PV as function of the angular scale . We also discuss the physical properties of the signal, and derive the covariance matrix. In Sec. 4 we introduce the datasets used in this work, and in Sec. 5 we present our results for the expected cumulative signal-to-noise ratio (S/N) for existing and forthcoming surveys. Finally, in Sec. 6 we discuss our findings and address future shortcomings.

Theory of angular correlations for Cosmological observables
In this section, following closely Ref. [40], we briefly review the theory of angular statistics for cosmological observables. To begin with, we note that any observable O α (n, z) can be expanded in spherical harmonics with respect to its angular position on the skyn: The correlation function between the spherical harmonic coefficients a α lm 's defines the angular cross spectrum, On the other hand, our theoretical predictions are usually derived from some complicated system of differential equations which are often more conveniently handled in Fourier space rather than in real space. The Fourier transform of O reads, 3) The plane waves may be rewritten using the spherical harmonics addition theorem, where the j are the spherical Bessel functions of the first kind. After substitution in Eq. (2.1), exploiting the orthogonality of the spherical harmonics, we finally obtain a α m (z) = 4πi For cosmological purposes, all the observables of interest can be split into a time and scale (and observable) dependent transfer function T α (k, z) and a primordial perturbation R(k) as O α (k, z) = T α (k, z)R(k). We assume that the primordial perturbation is Gaussian [41][42][43] and satisfies where P R (k) is the power spectrum of the primordial perturbation. The above expression implicitly assumes that the primordial perturbation is adiabatic, i.e. induces inhomogeneities in the spatial curvature. This assumption is reasonable in light of the Planck constraints on entropic perturbations [44], even though certain forms of isocurvature perturbations, for example compensated ones, are not well constrained by Planck [45]. Note that by construction we are assuming that different modes k, k are independent of each other. Substituting expressions 2.5 and 2.6 into 2.2 we arrive at In actual measurements, rather than O, one usually observes its average over a redshift bin weighted by some kernel W (z). Accordingly, this can be absorbed in to the definition of the angular power spectrum where the kernel ∆ α is defined as

Cross correlation between CMB lensing and Peculiar Velocities
The general derivation in the previous section worked primarily with redshifts, as this is the observable property of galaxy surveys. However, for the purposes of both CMB lensing and peculiar velocity theory, it is useful to consider distances. As such, for the remainder of this work we provide expressions in terms of the comoving distance where Ω m and Ω Λ are the matter and cosmological constant energy densities and where H 0 = H(z = 0) is the value of the Hubble factor today. The necessary factors for converting between integration over z and integration over χ have generally been absorbed into the definition of the window function W (χ) below.

Theoretical prediction
The Einstein field equations for scalar perturbations of a flat FLRW background in the Newtonian Gauge give, at linear order, the following relation between the peculiar velocity field v(r) and the dark matter density contrast δ(r) at position r: where a is the scale factor, H the Hubble factor and f = d ln D/d ln a the growth rate, which is itself a function of the growth function D, the growing mode of the linear density contrast. These are in general functions of the redshift, but in the following we will omit their explicit dependence. In Fourier space, the relation between velocity and density becomes In current practical experiments, however, one can usually measure only the radial projection of the velocity field, where k ≡ |k| andr ≡ r/|r|. The Fourier transform of this reads The above expression may be simplified, noticing that from which, using again the spherical harmonics addition theorem, the spherical harmonic coefficients for the radial velocity field may be written as a u m (χ) = (4πi )aHf Projecting the above along the line of sight we may define the radial velocity kernel ∆ u (k) as where we have introduced the window function W u (χ), W u (χ) = Hf a dn dχ , (3.9) and the distribution of sources dn/dχ. For comparison, the kernel for the galaxy distribution, neglecting the magnification bias, can be written where b(χ) is the linear galaxy bias. The CMB lensing signal, following Ref. [12], can be expressed in terms of the convergence field κ(n) in the directionn: where c is the speed of light and χ CMB the comoving distance to the surface of last scattering. From the above we can define the convergence kernel ∆ κ (k) and the convergence projection kernel W κ (χ) as: Combining expressions (3.8) and (3.12), we finally obtain an expression for their angular cross correlation which in terms of the window functions W i becomes, Note that since we are restricting our analysis to linear scales where Eq. (3.2) is a valid approximation, we have defined the matter power spectrum P m (k, Looking at the cosmological parameters entering in the kernels for the galaxy density, velocity and CMB convergence, we can see that the benefit of measuring the cross-correlation C uk is that it provides additional constraining power on the growth rate of structure compared to measuring the auto-correlation of the PVs, C uu , alone. It is also independent of galaxy bias. 2 Furthermore, as the CMB-lensing kernel contains information on the matter density, this may allow one to partially break the degeneracies between different models of gravity, with different growth rates of structure, and the matter density. This is simplest to see if we consider the common 'γ' parameterisation of the growth rate, f (z) = Ω γ m (z), where γ = const ≈ 0.55 for General Relativity, but differs in other gravitational theories [47,48]. Including measurements of C uk alongside the auto-correlations of the velocity or CMB-convergence clearly provides a route to further break the degeneracy between Ω m and γ(z) in the above expression.

Peeking "beyond" the survey
An interesting feature of the radial PV field is that it is effectively sensitive to the CDM distribution outside the scope of the PV survey itself. The same is true for its cross correlation with the convergence field κ, as reflected by the derivatives of the spherical Bessel function appearing in Eq. (3.14). This and other features of the C uκ are more easily understood, at qualitative level, estimating Eq.(3.15) using the Limber approximation to compute the integrals over χ 1 and χ 2 (see Appendix A for the details of the derivation): (3.16) As we are going to show in Sec. 5, most of the signal in C uκ comes from large angular scales, i.e. low-modes, for which the Limber approximation is inaccurate. With the exquisite precision expected from forthcoming cosmological surveys, alternative approximation schemes are currently developed in order to evaluate this type of integrals with greater accuracy, see for example Refs. [49,50]. However, all the result presented in this paper were obtained integrating numerically Eq.(3.15), without the use of any approximation. On the other hand, a number of physical properties of the signal are straightforwardly drawn from Eq. (3.16). For example, one can appreciate that the first term within the square brackets combines the κ and u fields evaluated at χ ± = ( ± 1/2)/k, which for a given are separated by a comoving distance proportional to k −1 . As an illustrative example, consider a galaxy in the north polar direction at a comoving distance χ − ≈ 300h −1 Mpc, i.e. at the edge of current PV surveys. According to Eq. (3.16), the radial peculiar velocity of this galaxy is correlated with the lensing induced on the trajectory of a CMB photon detected at the equatorial plane, i.e. for = 2, due to the CDM distribution at a comoving distance of almost χ + ≈ 500h −1 Mpc. 3 Put another way, the same distribution of matter located beyond the edge of the survey is responsible both for the large-scale motions of galaxies within the survey and the lensing of CMB photons, giving rise to a non-zero correlation. A schematic representation of the above effect is given in Fig. 1. Figure 1. Schematic representation of a lensed CMB photon γ (in yellow) and of the radial peculiar velocity u (in red) of a galaxy at the edge of the PV survey, i.e. at a distance χ survey , separated by an angular scale . According to Eq. (3.16), their correlation includes a contribution from the matter density distribution (in blue) at a distance

Overlap of the tracers
We easily realize from Eq. (3.14) that, for a given , the amount of correlation between the u and κ fields is essentially given by the overlap of the two kernels ∆ u ∆ κ weighted by the power spectrum P (k). In Fig. 2 we compare the product of these kernels with the galaxy-lensing and velocity-velocity ones, given as a function of the scale k at fixed angular mode for an idealized uniform distribution of sources which extends up to a distance χ(z max ), with χ(z max ) being the maximum distance probed by the survey. Note that these curves peak around k = /χ(z max ) and quickly decay for smaller values of k, as expected in virtue of the Limber approximation. Indeed, using the latter approximation the galaxy density kernel may be written: which for a uniform distribution of sources is roughly constant and non-vanishing only for z ≤ z max , whereas the radial PV kernel may be written: It is straightforward to realize that, since the lensing convergence kernel W κ at small redshift grows roughly linearly (see the dashed line in Fig. 3), the product ∆ g ∆ κ is maximum at k = ( + 1 2 )/χ(z max ). A similar statement holds true for the product between the velocity and lensing kernels. Indeed, the difference between the two terms in Eq. (3.18) for large is always positive and close to zero unless ( − 1 2 )/χ(z max ) < k < ( + 1 2 )/χ(z max ), for which the second term is vanishing due to the lack of observed sources beyond χ(z max ).
Notice that the decaying of the peaks at large scales, i.e. towards smaller k, for the red curves in Fig. 2 is slower compared to the green and blue curves, and reflects the lack of precision of the Limber approximation used in Eqs. (3.17),(3.18) for small . Another interesting feature which characterizes the product ∆ κ ∆ u is that, for a given , it oscillates between positive and negative values. As one may deduce from Fig. 2, its integral over k is essentially dominated by the value of ∆ κ ∆ u at k = /χ(z max ). We also notice that, since the matter power spectrum is peaked around the scale of the horizon at matter-radiation equality, k = k eq and the amplitude of ∆ u for a uniform and constant distribution of sources 4 is a decreasing function of , the correlation between the two fields is stronger at angular scales larger than ≤ k eq χ(z max ). On the other hand, at distances smaller than χ(z max ) the fluctuations of the radial velocity field average to almost zero. We conclude that most of the cosmological information from this correlation function comes from sources at the edge of the peculiar velocity survey, because of the lack of observed sources beyond it. For this reason, the expected signal is significantly smaller than, for example, the one from the galaxy density and CMB 4 i.e for a distribution of sources such that, given the small redshift window considered here, lensing cross correlation which contains contributions from all the sources within the galaxy survey.

Covariance
To assess quantitatively the statistical significance of the cross correlation detection (3.14) with future realistic datasets we need to associate to it a noise, i.e. we need to build a covariance matrix. This is defined as the non-connected four-point correlation function Cov(C uκ , C uκ ) ≡ Σ uκ = aũ , * m , aκ m , aũ m , aκ , * m , where we have defined our measured radial velocity and convergenceũ,κ as with ξ k,u being the random (Gaussian) noises. Expanding in spherical harmonics and making use of the Wick theorem (hence assuming that both κ and u are Gaussian fields) we can write the covariance as (see Eq. 15 of [12]) where f sky is the sky fraction where PV and CMB observations overlap, and N κκ , N uu are the noise spectra of the two fields. We assume that the shot noise relative to the radial velocity field, being induced by the underlying galaxy distribution, is isotropic (does not depend on l) and can be estimated (see Appendix B) as wheren is the angular number density of galaxies, and α is a scaling factor that depends of the type of tracer used for the peculiar velocity measurements. More details on these scalings are given in Section 4. For comparison, the usually assumed noise spectrum for the galaxy angular power spectrum is N gg = 1/n; the peculiar velocity field contains the same contribution from shot-noise as the galaxy density field, but has an additional term arising from the typical uncertainties in each PV measurement (αH 0 χ) weighted and integrated along the line-of-sight. The noise power spectrum of the reconstructed CMB convergence field κ depends non trivially on the specifics of the CMB experiment and on the particular reconstruction technique. The most commonly adopted method, the quadratic estimator, is based on the idea that lensing modifies the CMB statistics and induces a correlation between previously independent CMB harmonic modes [51,52]. By examining the correlation between modes of two CMB maps X, Y ∈ T, E, B, we can obtain a noisy estimate of φ XY m = 2 ( +1) κ XY m with associated noise power spectrum 5 where g XY L and f XY L are weight functions that depend on the power spectrum of the CMB fields C XY L and their overall noise levels N XY (see [53] for the exact expressions). The different temperature and polarization based estimators can be combined into a 5 The expression is strictly valid for the auto-spectrum of φ XY , i.e. φ XY φ XY, * , for the general form see [53]. minimum-variance estimate, which is the nominal lensing reconstruction we assume for the forecasts presented in this work. The CMB lensing reconstruction noise curves for the surveys considered in this work, together with the signal power spectrum, are shown in Fig. 4. Quadratic estimators will become sub-optimal at the instrumental noise levels soon-to-be reached by the most sensitive experiments, but a variety of methods based on maximum-likelihood and Bayesian techniques are being developed to restore optimality [see, e.g., 54-56].

Peculiar velocity and CMB datasets
So far we derived a theoretical prediction for the cross correlation signal and discussed its physical properties assuming idealized distributions of PV sources and lensed photons. To assess its potential as a cosmological probe, however, we should consider realistic datasets from existing and forthcoming experiments. In this section, we introduce the various PV and CMB catalogs we use to derive the results presented in Sec. 5. The number and the distribution of direct and reconstructed PV sources are summarized in Fig. 3 and Table 1. The signals and the noises for the radial velocity auto-correlation function and those of the various CMB lensing experiments are reported in Fig. 4.

PV Survey
Area (deg 2 ) Depth (z max )n (sr −1 ) α  observatory, in Chile. Whilst having similar uncertainty to DESI (α ≈ 0.2 in Eq. (3.23)), 4MOST will improve on the sky coverage, which will reach roughly 17000 deg 2 , and on the number of targets, expected to be of order ≈ 450 × 10 3 with redshift depth z ∼ 0.15. This survey is highly complementary to DESI in that it covers mostly non-overlapping parts of the sky, and so these to datasets could be combined to prove an area of order of 30000 deg 2 .
LSST The Legacy Survey of Space and Time (LSST) 8 is expected to be operational and running by 2024 at the Vera C. Rubin Observatory [57], currently under construction in Chile. Over it's 10-year program, it is expected to detect millions of Type Ia supernovae, of which some small fraction will have sufficient lightcurve measurements (from LSST or other follow-up campaigns) and host galaxy redshifts/properties to be cosmologically useful. Following [32], we consider a hypothetical 10-year survey consisting of ≈ 160 × 10 3 galaxies hosting Type Ia Supernovae with threshold J-band magnitude J < 19, over a an area of ∼ 18000 deg 2 and a redshift depth of order z ∼ 0.5. Such a number of objects may be feasible to obtain given already planned or ongoing spectroscopic surveys within the LSST footprint (for instance 4MOST mentioned previously) Although  relatively rare, the benefit of using Type Ia supernoave peculiar velocities is that the uncertainty is much better, scaling as 5-10% of the distance to the target. In this work we use an optimistic (but not overly so) value of α = 0.05 in Eq. (3.23).

Reconstructed PV
As discussed above, ongoing and future peculiar velocity surveys such as DESI, 4HS, and LSST will be an order of magnitude larger than previous surveys. However, with the current available methods for direct PV measurements (e.g. Tully-Fisher, Fundamental Plane, and SN Ia) the typical uncertainty scales as a considerable fraction of the distance. Thus, beyond redshift z ≈ 0.1, the uncertainty can become an order of magnitude bigger than the peculiar velocity itself, the latter being typically of order a few hundred km s −1 . This stands as a stumbling block in extending the PV surveys to higher redshifts, where the overlap with the CMB lensing kernel, and hence the signal in C uk is largest. Luckily it is possible to to infer indirectly the PV field from the underlying density field. Indeed, in the linear regime the peculiar velocity is attributed to the growing clustering of the total matter (both dark and luminous matter) via Eq. (3.2), which can be integrated to reconstruct the PV field. This equation however suffers a few limitations: first, we can only measure the galaxy density δ g and not the total density δ; second, PV are sensitive to scales beyond the redshift survey limit where we have no direct information. We can overcome these obstacles by assuming that linear biasing holds, δ = δ g /b where b is the linear biasing parameter, and adding an extra free parameter V ext to Eq. 3.2, accounting for the structures outside the survey limits. We can then readily integrate Eq.(3.2) to get the PV field as a function of V ext and b. These parameters can be calibrated using a small number of directly measured peculiar velocities, providing a reconstructed PV map. This method has a long history, see for example Refs. [58][59][60], and was recently used by [61] to reconstruct the velocity field in the local universe using positions of ≈ 70000 galaxies calibrated with ≈ 3000 directly observed PVs. The final map depends on the smoothing scale chosen for the reconstruction, and Ref. [61] concludes that a smoothing length of 4 Mpc h −1 gives the best peculiar velocity uncertainty of 250 km s −1 for the field galaxies and 150 km s −1 for groups. The assumption of a constant full-sky noise for reconstruction is clearly unrealistic, given the local environmental dependence of the noise in the density map from which the reconstruction is derived. On the other hand, the choice of a constant dispersion of 250 km s −1 is a conservative one, if compared for example with the average dispersion obtained from the space-dependent noise map of Ref. [62] (which includes the volume covered by [61]).Thus, even if slightly unrealistic, we expect our constant noise map to overestimate the typical errors in the reconstruction, making the conclusions we draw from it robust.
In this work we assume that their estimation of 250 km s −1 with a smoothing scale of 4 Mpc h −1 can be extrapolated up to redshift z < 0.15 covering the full sky. Redshift surveys covering this volume are expected to be available in the near future from surveys such as DESI, 4HS and others not mentioned above. Thus our idealized reconstructed PV catalog contains ≈ 1.6 × 10 7 objects (implyingn = 1.2 × 10 6 ) with a constant dispersion σ rec = 250 km s −1 replacing the factor αH 0 χ → σ rec in Eq. (3.23).

CMB lensing catalogs
PL18 The European Space Agency's Planck spacecraft 9 operated from 2009 to 2013 and performed full-sky measurements of the CMB intensity and polarization anisotropies across nine frequencies between 30 and 857 GHz with the High Frequency Instrument (HFI) and the Low Frequency Instrument (LFI). The final data were released in 2018, and contain the lensing potential map reconstructed from foreground-cleaned CMB temperature and polarization maps [63]. The CMB lensing convergence map covers roughly 70% of the sky (we assume an area of 27500 deg 2 in this work). The S/N forecasts below are obtained using the true CMB lensing noise curve measured from the data and available on the Planck Legacy Archive. 10 SO The Simons Observatory is a new-generation ground-based array of millimeterwave telescopes currently under construction in the high Atacama Desert of Chile [64] with first light planned for 2023. It will be equipped with a 6-meter Large-Aperture Telescope (LAT) and three smaller telescopes (SATs) and will image the CMB anisotropies over about 40% of the sky (16500 deg 2 ) in six frequencies between 27 and 280 GHz. To predict the feasibility of a possible detection of the cross-correlation signal, we use the official CMB lensing noise curve 11 released by the SO collaboration calculated assuming foreground cleaned CMB maps and using CMB multipoles between 30 ≤ ≤ 3000 in temperature and 30 ≤ ≤ 5000 in polarization.
CMB-S4 The fourth-generation ground-based CMB experiment, CMB-S4, is currently in its design stages and expected to start operations later this decade [65]. In its current configuration, 21 telescopes between the South Pole and the Atacama Desert in Chile will survey roughly 70% of the sky (again, in this work we assume 27500 deg 2 ). Relevant to this work, the LAT survey will use 6-metre class telescopes in six frequency bands from 30 to 270 GHz. Adopting the individual frequency noise levels from [56], we calculate the CMB lensing reconstruction noise curve from component-separated CMB using temperature and polarization information (in the 30 ≤ ≤ 3000 and 30 ≤ ≤ 5000 ranges respectively).

Results
The main results of this work, i.e. the cumulative Signal-to-Noise ratio (S/N) at angular scales 3 ≤ ≤ 200 for different PV and CMB surveys are reported in Fig. 5 and in Table 2. In all cases, we assume complete overlap between the PV and CMB lensing surveys, and take the smallest area of the two when computing the sky fraction f sky . The biggest angular scale that can be be probed by all the surveys considered in this work corresponds to a lower bound on the the angular mode = 3, which allows for a fairer comparison of the cumulative S/N from different surveys. It is clear from Fig. 5 that the cumulative S/N starts to plateau at small scales, i.e. large multipole . One cause for this is due to the lack of intrinsic signal at these scales, which can be seen by looking at Fig. 4, and is as expected by virtue of Eq. (3.16) which explicitly shows that C uκ scales with −1 . This in turns implies that a more conservative choice of angular scales without the first low-modes results in a significant degradation of the S/N, as one can appreciate comparing the upper and lower panels of table 2. We dedicate the rest of this section to a more thorough discussion of our results, particularly focused on the impact of the PV errors. CosmicVariance Planck CMB S4 Figure 5. The cumulative signal-to-noise ratio, with ≥ 3, as a function of maximum -mode for different representative combinations of PV and CMB dataset. The solid lines show the sample variance limit -the case where the shot/instrumental noise in both the CMB lensing convergence and peculiar maps is assumed to be zero.

Cross Correlation using direct PV measurements
We see from Table 2 that in the best case scenario, obtained by combining the direct PV measurements with the convergence map forecast for CMB S4, the cumulative S/N is above the 3σ threshold and therefore we could potentially marginally detect the cross correlation signal (see the dashed red line in Fig. 5). On the other hand, if we consider idealized PV surveys with no error (so that the only source of uncertainty is the sample variance), the cumulative S/N becomes greater than ≥ 10 already at relatively small ≈ 20 (see the solid red line in Fig. 5). One can identify the causes of such a degradation of the S/N ratio by looking at the terms entering the covariance matrix of Eq. (3.22) in Fig. 4. We see, comparing the solid and dotted lines in the upper right panel, that the Planck noise is always comparable or bigger than the intrinsic signal, whereas the expected noise from SO and CMB S4 is always smaller. This explains the difference between the first and the other columns in Table 2. On the other hand, by looking at the upper left panel of Fig. 4, we easily realize that the largest source of uncertainty comes from the intrinsic dispersion of PV measurements, for which the noise is bigger than the signal's theoretical prediction already for ≥ 5. This noise is largely dominated by the uncertainty on the position of the source, usually derived from scaling relations such as the Fundamental Plane (FP) or the Tully-Fisher (TF) ones. As a result, the dispersion on direct PV measurements usually scales with the distance to the source, leading to a quick degradation of the S/N. This explains the earlier flattening of the dashed and dotted red curves in Fig. 5 compared to the blue ones. In Fig.6 one can appreciate how the S/N changes as a function of α for a DESI-like survey combined with Planck or CMB-S4.

Cross correlation using reconstructed PV
We know that at redshift 0 ≤ z ≤ 1, the lensing kernel grows monotonically (see for example the dashed line of Fig. 3), hence the cross correlation signal should increase with PV measurements at higher redshifts. On the other hand, as previously discussed, the noise associated to direct PV measurements scales with the distance and causes a worsening of the overall S/N. To understand how this picture would change with PV dispersions without such scaling, we studied the cross correlation with a PV map obtained through a reconstruction technique, for which the dispersion is a constant determined by the choice of the reconstruction smoothing scale. For an idealized reconstructed PV field which extends up to redshift z = 0.15, with a smoothing scale of 4 Mpc h −1 and a constant dispersion σ rec = 250 km s −1 we see from Fig. 5 that the cumulative S/N improves greatly and approaches the sample variance limit. To estimate the dependence of the signal from the constant dispersion assumed for the reconstructed map we plot in Fig. 7 the cumulative S/N as a function of the value of the constant dispersion σ rec , from which we see that the cumulative S/N is still ≥ 10 even if the constant noise is of order σ rec ∼ 10 3 . It is important to stress that the reconstructed PV map depends somehow on the galaxy density field and the galaxy bias b(z), and therefore cannot be considered as a completely independent tracer of the underlying DM distribution. However, this exercise also shows the potential of C uκ as an unbiased cosmological probe if direct PV measurements with the precision of the reconstructed ones will be available in the future. Some examples could be through the construction of direct PV surveys with higher number density than we have considered here, the use of ancillary properties such as metallicity or age in the standard Tully-Fisher or Fundamental Plane relations to reduce their intrinsic scatter [66][67][68], the use of 'twinned' supernovae to obtain more accurate distance moduli [69], or alternate distance indicators such as Brightest Cluster Galaxies [70], Type-II supernovae [71], or even gravitational waves [72]. However, while promising, more work is needed in this area to realise the potential of, and implement, these improvements in future surveys. rate f , whereas the lensing convergence to the matter distribution Ω m , making them a promising tool to explore different gravitational models and current cosmological parameters tensions. In this work, for the first time, we assess the significance of their statistical cross correlation from current and forthcoming experiments. Because of the geometry of the problem, one would expect the signal to be intrinsically small. Indeed, the CMB lensing peaks at relatively high redshift z ∼ O(1), whereas peculiar motions observations are available only at closer redshifts z ∼ O(0.1). However, contrary to these expectations, our theoretical prediction shows that with perfect experiments, i.e. at the sample variance limit, the signal would be significantly detectable with a cumulative S/N ≥ 30σ (standard deviations) at angular scales 3 ≤ ≤ 100. If realistic noises are taken into account, however, the S/N drops for direct PV measurements to a few σ's, making their detection marginal. Looking at Eq. (3.22) and Fig. 4 it is straightforward to realize that the observed PV uncertainties are responsible for such a degradation of the cumulative S/N, since the noise starts to dominate over the signal already at relatively large scales (small ). This is due to the fact that the error on the source's PV is usually proportional to its comoving distance to the observer, and thus increases with the redshift. Another option is to use a reconstructed peculiar velocity map, for which the uncertainty is approximately constant with the redshift, depending on the choice of the reconstruction smoothing scale. For a reconstruction that extends up to redshift z ≤ 0.15, a smoothing scale of 4h −1 Mpc (i.e. one likely to be feasible with DESI and 4HS data) and CMB S4 observations, the dashed blue line in Fig. 5 shows that the cumulative S/N approaches the sample variance limit. Whilst a reconstructed PV map would not be completely independent of the galaxy bias b(z), the above exercise shows the expected significance of such a cross-correlation if the errors on PV observations do not grow with redshift.
In summary, this paper shows that the radial peculiar velocities of galaxies at low redshift and the convergence map of CMB photons on angular scales O(1) ≤ ≤ O(100) are correlated, and such correlation might be detectable with forthcoming surveys in the near future. It would be interesting to assess the constraining power of such a probe for cosmological inference of ΛCDM derived parameters like the growth rate f or the Hubble factor today H 0 , as well as for testing different theories of gravity and cosmology beyond the standard model. We will address these compelling questions in future works.

Acknowledgement
We are grateful to Hayley Valiantis for their support during the preparation of this paper. Part of the work was financially supported by the Australian Government through the Australian Research Council Laureate Fellowship grant FL180100168.

A Limber approximation for C uκ
The spherical Bessel functions of order n may be written using which allows us to rewrite Eq. (3.14) as the sum of two terms: C uκ ≡ 2 π dk dχ 1 dχ 2 W u (χ 1 )W κ (χ 2 ) j −1 (kχ 1 ) − + 1 kχ 1 j (kχ 1 ) j (kχ 2 ) P m (k, χ 1 , χ 2 ) . B Shot noise for the projected radial peculiar velocity field The effect of noise in the radial peculiar velocity field can be computed following the same derivation as for any observable in Section 2. If we define the observed peculiar velocity field asũ(n, χ) = u(n, χ)+ξ u (n, χ) then the observed angular velocity-velocity power spectrum can be writtenC uu = C uu + N uu where N uu denotes the contribution from the noise in the velocity field N uu = dΩndΩn Y * m (n)Y m (n ) ξ u (n, χ 1 )ξ u (n , χ 2 ) . (B.1) In practice, the angular power spectrum is computed by projecting the observable over the redshift bin. For the radial peculiar velocity this projection kernel is dn/dχ. This act of projection also applies to the noise spectrum, such that N uu ⇒ N uu = dn dχ 1 dχ 1 dn dχ 2 dχ 2 dΩndΩn Y * m (n)Y m (n ) ξ(n, χ 1 )ξ(n , χ 2 ) (B.2) As discussed in [28,32,74,75] the error in the radial peculiar velocity field generally scales as ξ u (n, χ) = αH 0 χ/ √ N where the 1/ √ N term arises from the averaging of N independent radial peculiar velocity measurements and typically α = 0.2 for Tully-Fisher or Fundamental Plane distance measurements and α = 0.05 for Type Ia Supernovae. Crucially, the errors typically depend on the distance to the object, but not on the direction on the sky. Hence, if we substitute this expression into Eq. B.2 we can integrate out the spherical harmonics where the resulting 4π has been used to convert the number of galaxies to an angular number densityn. This expression matches that in Eq. 3.23 from the main text.