Cross-correlation between the CMB lensing potential measured by Planck and high-z sub-mm galaxies detected by the Herschel-ATLAS survey

We present the first measurement of the correlation between the map of the cosmic microwave background (CMB) lensing potential derived from the \emph{Planck} nominal mission data and $z\gtrsim 1.5$ galaxies detected by the \emph{Herschel}-ATLAS (H-ATLAS) survey covering about $600\,\hbox{deg}^2$, i.e. about 1.4\% of the sky. We reject the hypothesis that there is no correlation between CMB lensing and galaxy detection at a $20\,\sigma$ significance, checking the result by performing a number of null tests. The significance of the detection of the theoretically expected cross-correlation signal is found to be $10\,\sigma$. The galaxy bias parameter, $b$, derived from a joint analysis of the cross-power spectrum and of the auto-power spectrum of the galaxy density contrast is found to be $b=2.80^{+0.12}_{-0.11}$, consistent with earlier estimates for H-ATLAS galaxies at similar redshifts. On the other hand, the amplitude of the cross-correlation is found to be a factor $1.62 \pm 0.16$ higher than expected from the standard model and also found by cross-correlation analyses with other tracers of the large-scale structure. The enhancement due to lensing magnification can account for only a fraction of the excess cross-correlation signal. We suggest that part of it may be due to an incomplete removal of the contamination of the CIB, that includes the H-ATLAS sources we are cross-correlating with. In any case, the highly significant detection reported here using a catalog covering only 1.4\% of the sky demonstrates the potential of CMB lensing correlations with submillimeter surveys.


INTRODUCTION
Cosmological observations carried out in the last two decades have enabled the establishment of the standard cosmological model. In this picture, observed galaxies form in matter over-densities that are the result of the growth, driven by gravitational instabilities in an expanding Universe, of primordial inhomogeneities generated during an inflationary epoch. A picture of primordial inhomogeneities at an early stage of their evolution is provided by observations of the cosmic microwave background (CMB) anisotropy.
However, this picture is to some extent distorted by interactions of the CMB photons with matter inhomogeneities encountered during their travel from the lastscattering surface to the observer. On the other hand, these effects are a useful source of information on the large-scale structure of the Universe. One of these effects is gravitational lensing, causing small but coherent deflections of the observed CMB temperature and polarisation anisotropies, with a typical amplitude of 2 . Specific statistical signatures of lensing enable the reconstruction of the gravitational potential integrated along the line of sight from observed CMB maps (Hu & Okamoto 2002;Hirata & Seljak 2003).
In recent years CMB lensing has been measured in a number of CMB experiments. The first detections were made via cross-correlations with large-scale structure probed by galaxy surveys (Smith et al. 2007;Hirata et al. 2008;Feng et al. 2012;Bleem et al. 2012;Sherwin et al. 2012;Geach et al. 2013;Holder et al. 2013). The higher sensitivity and resolution of recent CMB experiments, such as ACT, SPT and Planck, have enabled an internal detection of lensing using CMB data alone (Das et al. 2011;Keisler et al. 2011;Das et al. 2014;van Engelen et al. 2012); the measurement with the highest signal-to-noise ratio, around 25σ, was reported last year by the Planck team (Planck Collaboration XVII 2013).
As already mentioned, the CMB lensing potential is an integrated measure of the matter distribution in the Universe, up to the last-scattering surface. As illustrated by Fig. 1, it has a broad kernel, peaking at z 2 but slowly varying from z 1 to z > ∼ 4. A powerful way to reconstruct the evolution of the gravitational potential with cosmic time is to turn to cross-correlations with source samples covering narrow redshift ranges. Several catalogs have been exploited for this purpose, such as those from the NRAO VLA Sky Survey (NVSS), the Sloan Digital Sky Survey (SDSS), the Wide Field Survey Infrared Explorer (WISE). These surveys cover large areas of the sky but detected sources are mostly at z < ∼ 1. The Herschel Terahertz Large Area survey (H-ATLAS; Eales et al. 2010) has provided large samples of galaxies up to substantial redshifts (Lapi et al. 2011;González-Nuevo et al. 2012). In this paper we present the first investigation of the cross-correlation between the CMB lensing potential measured by Planck and Herschel -selected galaxies with estimated redshifts z > ∼ 1.5, i.e. at redshifts higher and closer to the peak of the lensing potential kernel than those of source samples considered so far. Our choice of restricting the analysis to z > ∼ 1.5 has a twofold motivation. First, since we aim at reconstructing the evolution of the lensing potential at higher redshifts than done with other galaxy samples, it is expedient to remove the dilution of the signal by low-z sources. Second, as shown by Lapi et al. (2011) and González-Nuevo et al. (2012), the adopted approach for estimating photometric redshifts becomes unreliable at z < ∼ 1. The outline of this paper is as follows. In Section 2 we describe the theoretical background while the data are introduced in Section 3. The estimator of the crosscorrelation power spectrum and the simulations used for validation of the algorithm and the error estimation are presented in Section 4. The measured auto-and crosspower spectra, as well as the null tests, are reported in Section 5. In Section 6 we analyze the constraints on the galaxy bias and in Section 7 we discuss the potential systematic effects affecting the cross-correlation. Finally in Section 8 we summarize our results.
Throughout this paper we adopt the fiducial flat ΛCDM cosmology with best-fit Planck + WP + highL + lensing cosmological parameters as provided by the Planck team in Planck Collaboration XVI (2013). Here WP refers to WMAP polarization data at low multipoles, highL refers to the inclusion of high-resolution CMB data of the Atacama Cosmology Telescope (ACT) and South Pole Telescope (SPT) experiments and lensing refers to the inclusion of Planck CMB lensing data in the parameter likelihood.

THEORETICAL BACKGROUND
The effect of gravitational lensing on CMB photons can be described as a re-mapping of the unlensed temperature anisotropies Θ(n) by a two-dimensional vector field in the sky, namely the deflection field d(n) (Lewis & Challinor 2006): whereΘ(n) are the lensed temperature anisotropies and φ(n) is the CMB lensing potential: In this equation, χ(z) is the comoving distance to redshift z, χ * is the comoving distance to the last scattering surface at z * 1090, H(z) is the Hubble factor at redshift z, c is the speed of light, and Ψ(χ(z)n, z) is the 3D gravitational potential at a point on the photon path given by χ(z)n. Note that the deflection angle is given by d(n) = ∇φ(n), where ∇ is the the twodimensional gradient on the sphere. Since the lensing potential is an integrated measure of the projected gravitational potential, taking the two-dimensional Laplacian of the lensing potential we can define the lensing convergence κ(n) = − 1 2 ∇ 2 φ(n) which depends on the projected matter overdensity δ (Bartelmann & Schneider 2001): The lensing kernel W κ is where Ω m and H 0 are the present-day values of the Hubble and matter density parameters, respectively. The galaxy overdensity g(n) in a given direction on the sky is also expressed as a line-of-sight integral of the matter overdensity, where the kernel is The galaxy overdensity kernel is the sum of two terms: the first one is given by the product of the linear bias b(z) and the redshift distribution dN/dz; and the second one takes into account the effect of gravitational magnification on the observed density of foreground sources (magnification bias; Ho et al. 2008;Xia et al. 2009). This effect depends on the slope, α(z), of their integral counts (N (> S) ∝ S −α ) below the adopted flux density limit. and cumulative signal to noise ratio (red lines; right axis) evaluated from min = 100 for fiducial models with b = 3 and α = 1 (no magnification, dashed lines) and α = 3 (solid lines).
Given the sharply peaked redshift distribution of our sources (see Fig. 1) we can safely assume a redshift-and scale-independent linear bias (b(z) = constant). Previous analyses of the clustering properties of sub-mm galaxies (Xia et al. 2012;Cai et al. 2013) indicate b 3 at the redshifts of interest here and we adopt this as our reference value.
Recent work by González-Nuevo et al. (2014) has shown that the magnification bias by weak lensing is substantial for high-z H-ATLAS sources selected with the same criteria as the present sample (see the sub-section 3.2). This is because the source counts are steep, although their slope below the adopted flux density limit (S 250µm = 35 mJy) is uncertain. The data (Béthermin et al. 2012) indicate, at this limit, α 2 while for the high-z galaxy sub-sample considered in this work we find α 3. In the following we adopt the latter as our fiducial value. The effect of different choices for this parameter value is examined in Section 7.
Since the relevant angular scales are much smaller than 1 radian (multipoles 100) the theoretical angular cross-correlation can be computed using the Limber approximation (Limber 1953) as where P (k, z) is the matter power spectrum, which we computed using the CAMB 14 code (Lewis et al. 2000). The non-linear evolution of the matter power spectrum was taken into account using the HALOFIT prescription (Smith et al. 2003); more extended discussion on the effect of the non-linear evolution in CMB lensing maps based on N-body simulations is carried out by Antolini et al. (2014). The CMB convergence, W κ (z), and the galaxy redshift distribution dN/dz of the sample analyzed in this work (see sub-section 3.2) are shown in Figure 1. Again under the Limber approximation the CMB convergence, C κκ , and the galaxy, C gg , auto-spectra can be 14 available at http://camb.info evaluated as: The mean redshift probed by the cross-correlation between CMB lensing and our sample is We can predict the signal to noise ratio (S/N) of the convergence-density correlation assuming that both the galaxy overdensity and the lensing fields behave as Gaussian random fields, so that the variance of C κg is where f sky is to the sky fraction covered by both the galaxy and the lensing surveys, N κκ is the noise of the lensing field, and N gg = 1/n is the shot-noise associated with the galaxy field, inversely proportional to the mean number of sources per steradian,n. The signal to noise ratio at multipole is then and the cumulative signal to noise ratio for multipoles up to max is In Figure 2 we show both the signal to noise per multipole and the cumulative one computed using the specifications for the Planck lensing noise (see sub-section 3.1) and the mean surface density of our source sample. It must be noted that, due to the limited area covered by the H-ATLAS survey (and split into 5 fields) the cross-correlation is only meaningful on scales below a few degrees. We have therefore limited our analysis to ≥ min = 100. This restriction prevents us from exploiting the peak at ∼ 100 of the signal to noise per multipole. The cumulative signal to noise ratio saturates at ∼ 1000. If b = 3 and α = 3 we expect S/N 6.

Planck data
We used the publicly released Planck CMB lensing potential map derived from the first 15.5 months of observations (Planck Collaboration XVII 2013). The Planck satellite observed the sky with high angular resolution in nine frequency bands, from 30 to 857 GHz (Planck Collaboration I 2013). The angular resolution (10 , 7 and 5 ) and the noise level (105, 45 and 60 µK arcmin) of the 100, 143 and 217 GHz frequency channels, respectively, make them the most suitable for estimation of the gravitational lensing potential. Nevertheless, the released map is based on a minimum variance combination  of the 143 and 217 GHz temperature anisotropy maps only, as adding the 100 GHz map yields a negligible improvement (Planck Collaboration XVII 2013). The maps are in the HEALPix 15 (Górski et al. 2005) format with a resolution parameter N side = 2048, corresponding to 50 331 648 pixels over the sky, with a pixel size of ∼ 1.7 .
The power spectrum of the lensing potential is very red and this may introduce a bias when we estimate it within multipole bins. To avoid this problem we decided to convert the lensing potential map, φ, into the convergence map, κ, which has a much less red power spectrum. This was done using the relation between the spherical harmonic coefficients of these quantities estimated on the full sky (Hu 2000) The convergence spherical harmonic coefficients were transformed to a map with resolution parameter N side = 512 corresponding to a pixel size of ∼ 7 . This resolution is sufficient for our analysis because the data noise level enables us to detect cross-correlations between the convergence and the galaxy density field only for angular scales larger than ∼ 20 ( 540). The convergence auto-power spectrum recovered on approximately 60% of the sky using a modified version of the mask provided by the Planck collaboration is shown in Fig. 3. There is good agreement with theoretical expectations up to 650. The discrepancies at higher multipoles are out of the range of interest for the present paper.

Herschel fields
We exploited the data collected by the Herschel Space Observatory (Pilbratt et al. 2010) in the context of the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010), an open-time key program that has surveyed 600 deg 2 with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) in five bands, from -Redshift distribution of H-ATLAS galaxies for the combined set of patches used in the analysis. The (blue) histogram is the empirical redshift distributions, the dashed (orange) line is the half-normal fit to dN/dz as described in text, while the solid (green) line represents the convolved dN/dz that takes into account errors on photo-z estimation and is used as the fiducial distribution in our analysis. The values of the parameters µ and σ given in the box are the best-fit values and are used in the analytic expression for dN/dz adopted in calculations.
100 to 500 µm. The H-ATLAS map-making is described by Pascale et al. (2011) for SPIRE and by Ibar et al. (2010) for PACS. The procedures for source extraction and catalogue generation can be found in Rigby et al. (2011), Maddox et al. (in preparation) and Valiante et al. (in preparation).
The z < ∼ 1 galaxies detected by the H-ATLAS survey are mostly late-type and starburst galaxies with moderate star formation rates and relatively weak clustering Guo et al. 2011). High-z galaxies are forming stars at high rates (≥ few hundred M yr −1 ) and are much more strongly clustered (Maddox et al. 2010;Xia et al. 2012), implying that they are tracers of large scale over-densities. Their properties are consistent with them being the progenitors of local massive ellipti- cal galaxies (Lapi et al. 2011). We aim at correlating high-z H-ATLAS galaxies with the Planck CMB lensing map.
Our final sample comprises a total of 99 823 sources, of which 9099 are in G09, 8751 in G12, 9279 in G15, 28245 in NGP and 44449 in SGP. The specifics of each patch are summarized in Table 1. The redshift distribution of the population is needed in order to predict the amplitude of the cross-correlation. Estimating the uncertainties in the redshift distribution due to photometric redshift errors is not a trivial task. In order to get a feeling of the impact of these errors on the estimated dN/dz, we convolved the Gaussian fit to the redshift distribution for the full H-ATLAS sample with a Gaussian kernel whose dispersion reflected the rms uncertainties on estimated redshifts as estimated by González-Nuevo et al. (2012). Then we cut the convolved redshift distribution at z = 1.5 and fitted the portion at higher z with a half-normal distribution normalized to unity: The redshift distribution of the galaxies in our map is shown in Fig. 4. We built an overdensity map at a resolution N side = 512 defined by where n(n) is the number of objects in a given pixel, and n is the mean number of objects per pixel.
The CMB convergence and galaxy overdensity fields in the different patches are shown in Fig. 5. We filtered out from these fields multipoles 400 dominated by noise.

Estimator
We computed the angular power spectra within the regions covered by the H-ATLAS survey using a pseudo-C estimator based on the MASTER algorithm (Hivon et al. 2002). These regions are inside the area used in the estimation of the CMB lensing map. For a survey that covers only a fraction of the sky different modes of the true cross power spectrum C κg are coupled (Hauser & Peebles 1973). The coupling can be described by the mode-mode coupling matrix M which relates the pseudo cross-spectrumC κg measured from the data: to the true power spectrum: However we cannot directly invert eq. 17 to get the true power spectrum, because for surveys covering only a small fraction of the sky the coupling matrix M becomes singular. To reduce the correlations of the C 's it is necessary to bin the power spectrum in . We used eight linearly spaced bins of width ∆ = 100 in the range 0 ≤ ≤ 800. Then, the estimator of the true bandpowersĈ κg L (hereafter C κg L denotes the binned power spectrum and L identifies the bin) is given bŷ Cross-power spectrum of simulated galaxy and lensing maps constructed with b = 3. The points connected by the solid blue line represent the binned input cross-spectrum while the average reconstructed spectrum from 500 simulations is shown by the orange points. Lower panel: Fractional difference between the input and extracted cross-spectra. Error bars obtained with the simulation covariance matrix (orange points) and with the analytical approximation (blue points) are shown for comparison. Middle. As in left plot, but for the galaxy auto-power spectrum. Right. As in left plot, but for the CMB convergence auto-power spectrum Here P L is the binning operator, Q L and B 2 are, respectively, the reciprocal of the binning operator and the pixel window function that takes into account the finite pixel size. Because of the small size of the sky area covered by the H-ATLAS survey the power spectrum for < 100 is very poorly estimated, we did not use it in our analysis. However, to avoid bias coming from the lowest order multipoles the first multipole bin is included in computation of the power spectrum i.e. the inversion of the binned coupling matrix K LL is performed including the first bin and pseudo-power spectrum for the first bin is used in the product given by eq. 18.
The main assumption in cross-correlation studies is that the noise levels related to the observables being analyzed are uncorrelated so that we do not need to debias the reconstructed cross-spectrum for any noise term. However, when dealing with auto-power spectra, such as C gg and C κκ , we have to correct the estimator given by eq. 18 in order to account for the noisê where Ñ gg MC and Ñ gg MC are the average noise pseudo-spectra estimated from the Monte Carlo (MC) simulations.

Covariance matrix
The errors on the cross power spectrum are described by the covariance matrix (Brown et al. 2005) where Cov κg is the pseudo-covariance matrix given by The corresponding covariance matrix of the galaxy autocorrelation is obtained by replacing in eq. (21) (23) The analytical expressions for the covariance matrices given above were used in the estimation of the galaxy bias and of the amplitude of the cross-correlation, presented in Section 6.
In order to validate the algorithms used for the computation of the estimators outlined in the previous section and to check that the cross-and auto-power spectra estimates are unbiased, we created 500 simulated maps of the CMB convergence field and of the galaxy overdensity field with statistical properties consistent with observations.
Using the theoretical spectra obtained with eqs. 7 and 8, we generated full sky signal maps injecting a known degree of correlation, so that the simulated CMB convergence and galaxy harmonic modes satisfy both the autoand the cross-correlations (Kamionkowski et al. 1997): For each value of and m > 0, ζ 1 and ζ 2 are two complex numbers drawn from a Gaussian distribution with unit variance, while for m = 0 they are real and normally distributed.
We also generated 500 noise realizations for both fields. To simulate Gaussian convergence noise maps we used the convergence noise power spectrum N κκ provided by the Planck team. Although this power spectrum is not sufficiently accurate for the estimation of the convergence power spectrum, as pointed out in the Planck Collaboration Products website, it should be sufficiently good for the cross-correlation analysis which is not biased by the noise term. For the same reason it is not crucial for our analysis to use the 100 simulations of the estimated lensing maps provided recently by the Planck team.
To take into account noise in the simulated galaxy maps, we proceeded in the following way. For each signal maps containing the galaxy overdensity, we generated a set of simulated galaxy number count maps, where the value in each pixel is drawn from a Poisson distribution -The CMB convergence -galaxy density cross-spectrum as measured from Planck and Herschel data. The data points are shown in blue, with error bars computed using the full covariance matrix obtained from Monte Carlo realizations of convergence maps. The theoretical spectra calculated with the bias values inferred from the likelihood analysis (as described in text) using the cross-correlation data only (solid red line) and the cross-correlation together with the galaxy auto-correlation data (dot-dashed green line) are also shown; we fix α = 3 in this analysis. The null (no correlation) hypothesis is rejected at the 20 σ level.
with mean λ(n) =n(1 + g(n)), wheren is the mean number of sources per pixel in a given H-ATLAS patch and g(n) is the corresponding simulated galaxy map containing only signal. The galaxy number counts map λ(n) was then converted into a galaxy overdensity map using eq. 15, substituting the real number of objects in a given pixel n(n) with the simulated one λ(n). Note that maps obtained in this way already include Poisson noise with variance N gg = 1/n. We applied the pipeline described above to our set of simulations in order to recover the input cross-and autopower spectra used to generate such simulations. The extractedĈ κg L ,Ĉ gg L , andĈ κκ L spectra averaged over 500 simulations are reported in Fig. 6. The mean band-power was computed as: where X, Y = {κ, g}, i refers to the i-th simulation, and N sim = 500 is the number of simulations. The errors were computed from the covariance matrix as and the covariance matrix Cov XY LL was evaluated from the simulations as (28) We also show, for comparison, the theoretical error bars obtained from eq. (10), modified to take into account the binning. They are in generally good agreement with the Monte Carlo error estimates, which, however, are slightly larger (by up to ∼ 25%).

CMB Convergence-Galaxy Cross-correlation
The recovered cross-spectrum is shown in Fig. 7. The error bars are estimated from cross-correlating 500 Monte Carlo realizations of simulated CMB convergence maps (containing both signal and noise) with the true H-ATLAS galaxy density map, as described in subsection 5.3. This method assumes that the two maps are uncorrelated; our error estimates are a good approximation, since both maps are very noisy and C κκ,tot C gg,tot (C κg ) 2 . We have also estimated the errors from crosscorrelations of 500 Monte Carlo realizations of simulated H-ATLAS galaxy density maps with the real Planck CMB convergence map. The former approach yields slightly smaller error bars, yet slightly larger than those estimated analytically (see Fig. 8).
We have exploited the simulations to build the covariance matrix, used to evaluate the probability that the measured signal is consistent with no correlation (our null hypothesis). As can be seen in Fig. 9, the covariance matrix is dominated by the diagonal components; however, off-diagonal components are non-negligible and have to be taken into account. The χ 2 was calculated as: For the analysis performed with the whole H-ATLAS sample we obtained χ 2 null = 83.31 for ν = 7 degrees of freedom, corresponding to a p-value p = 2.89 × 10 −15 . Since the χ 2 distribution has mean ν and variance 2ν, the null hypothesis is rejected with a significance of about 20 σ. This is the sum in quadrature of the significance of the correlation in each band-power, taking into account the correlations between different bins. The results of the χ 2 -Galaxy density auto-power spectrum for the whole sample of H-ATLAS galaxies. The data points are shown in blue, while the solid (red) line is the theoretical C gg evaluated for the best-fit value of the bias obtained using a likelihood analysis on the galaxy auto-spectrum data.
analysis for each patch are reported in Table 2.

Galaxy Auto-correlation
We also performed an analysis of the auto-correlation of Herschel galaxies on the different patches. The autopower spectrum measured for the complete H-ATLAS data set is shown in Figure 10. The error bars on the data points were evaluated from the diagonal part of the covariance matrix built from galaxy simulations with bias b = 3. The detected signal is highly significant (40 σ).

Null Tests
In order to verify our pipeline, and the reconstructed spectra, against residual systematic errors, we performed a series of null tests. The idea behind these null tests is to cross-correlate in turn the real map of one field with simulated maps of the other field. Since there is no common cosmological signal, the mean correlation must be zero.
To this end we cross-correlated our 500 simulated CMB lensing maps (containing both signal and noise) with the real H-ATLAS galaxy density contrast map and our 500 simulated galaxy maps constructed using b = 3 with the true Planck CMB convergence map. The error bars on the cross-power spectra were computed using the covariance matrices obtained from these simulations. As illustrated in Figure 11 in both cases no significant signal was detected. In the first test we obtained χ 2 = 7.2 with corresponding to a rejection probability of the null hypothesis (no correlation) p = 0.41, while in the second one we have χ 2 = 5.9 and p = 0.55.
A further test consisted in cross-correlating the galaxy distribution in one patch of the sky with the lensing map in another. We moved in turn the three H-ATLAS GAMA fields and the SGP field to the position of the NGP patch and shifted the NGP galaxies to the SGP area. Then we cross-correlated each shifted galaxy map with the convergence field in the same position. The errors on the cross-correlations were obtained as above. All of the cross-spectra are consistent with no signal.

CROSS-CORRELATION
We now discuss the origin of the cross-correlation signal, addressing its plausible cosmological origin. Following Planck Collaboration XVII (2013) we introduce an additional parameter, A, that scales the expected amplitude of the cross-power spectrum, C κg , of the Planck CMB lensing with the H-ATLAS galaxy overdensity map as A C κg L (b). Obviously, its expected value is 1. Because the theoretical cross-spectrum is also basically proportional to the galaxy bias, there is strong degeneracy between these two parameters. In order to break this degeneracy we use also the galaxy auto-power spectrum.
The best values of the amplitude and of the galaxy bias were obtained using the maximum likelihood approach. In the following we first describe the likelihood functions and present constraints on the redshift-independent galaxy bias and on the cross-correlation amplitude using galaxy auto-correlation data alone, cross-correlation data alone, and combining both datasets. In this analysis the cosmological parameters and the counts slope α are kept fixed to the fiducial values. In order to efficiently sample the parameter space, we use the Markov Chain Monte Carlo (MCMC) method assuming uninformative flat priors. For this purpose we employ EMCEE (Foreman-Mackey et al. 2013), a public implementation of the affine invariant MCMC ensemble sampler (Goodman & Weare 2010). In this paper, each quoted parameter estimate is the median of the appropriate posterior distribution after marginalizing over remaining parameters with uncertainties given by 16 th and 84 th percentiles (indicating the bounds of a 68% credible interval). For a Gaussian distributions, as is the case when combining both datasets, these percentiles correspond approximately to −1σ and +1σ values and median of the posterior is equal to mean and maximum likelihood value.
We assumed Gaussian likelihood functions for the cross-and auto-power spectra. For the galaxy autopower spectrum it takes the form b = 2.80 +0.12   Fig. 12.-Posterior distribution in the b − A plane with the 68% and 95% confidence contours (darker and lighter colour, respectively), together with the marginalized distributions of each parameter with 1σ errors shown by the dashed white lines, obtained combining the convergence-galaxy cross-correlation and the galaxy auto-correlation data for each patch. The solid red line represents the standard case in which A = 1 while α is set to 3 for the analysis.
where N L = 7 is the number of bins reconstructed and Cov gg LL is the covariance matrix given by eq.(23). Sampling this likelihood for the measured H-ATLAS galaxy power spectrumĈ gg L we obtained constraints on the galaxy bias. Estimated values of the bias for all as well as for separate patches are presented in Table 2. The results for the different patches are consistent within 2σ. The global value, b = 2.84±0.12, is consistent with earlier estimates (Xia et al. 2012).
We used the measured cross-spectra to constrain the b and A parameters in the same fashion. As noted above, the cross-spectra basically measure the product A × b. The likelihood function is given by: where Cov κg LL is the covariance matrix (eq. (21)). The results are shown in Table 2.
Finally, we studied the constraints on b and A by combining the cross-and galaxy auto-spectra data. For the joint analysis we used the Gaussian likelihood function that takes into account correlation between the cross and auto-power spectra in the covariance matrix.
The full 2-dimensional posterior distributions of the b and A parameters, as well as the marginalized ones obtained from this analysis, are shown in Fig. 12. Numerical values of the parameters are presented in Table  2. The constraint on the bias factor from the joint fit of the galaxy auto-correlation and of the cross-correlation power spectra, b = 2.80 +0.12 −0.11 , is consistent with earlier estimates (Xia et al. 2012). On the other hand, the crosscorrelation amplitude is A = 1.62±0.16 times larger than expected for the standard ΛCDM model for the evolution of large-scale structure. This is at odds with the results of the cross-correlation analyses presented in the Planck Collaboration XVII (2013) paper, which are consistent with A = 1 except, perhaps, in the case of the MaxBCG cluster catalog. Possible causes of the large value of A are discussed in the following section.

TESTING POTENTIAL SOURCES OF SYSTEMATIC ERRORS
The correlation between the CMB lensing potential and the distribution of high-z, sub-mm selected galaxies was found to be stronger than expected for the standard cosmological model. We consider now other possible (i.e. astrophysical) sources of the observed correlation. The first aspect we consider is the effect of the redshift dependence of the galaxy bias. Our choice of a constant b over the redshift range spanned by the H-ATLAS catalogue is obviously an approximation and the effective values of b may be different for the cross-and the galaxy auto-power spectra. To check the effect of this approximation on the estimates of C κg and C gg we have computed the effective values of the bias for the two cases using the bias evolution model b(z) from Sheth & Tormen (1999) for halo masses in the range 10 12 -10 13 M . We find that b κg eff is only slightly larger (by 6%) than b gg eff . Hence, considering a redshift-dependent bias factor would only marginally affect the expected crossspectrum.
Weak lensing by foreground structures modifies the observed density of background sources compared to the real one (magnification bias; Ho et al. 2008;Xia et al. 2009) and is especially important for high redshift objects. The effect on the galaxy over-density kernel is described by the second term on the right-side of eq. (6). The effect of the magnification bias on both C κg and C gg is illustrated in Fig. 13 where we show the expected power spectra for A = 1, b = 3 and three values of α: 1 (no magnification bias), 3, and 5. The impact of the magnification bias is clearly stronger for C κg .
Fitting the joint data for α = 1 we find b = 2.95 +0.12 −0.11 and A = 1.93 +0.18 −0.19 while for α = 5 b = 2.55 +0.13 −0.12 and A = 1.46 ± 0.14. The contours plots in the A − b plane are shown in Fig. 14. Higher values of α imply lower values of A, but even for α = 5 the data require A > 1.
Clusters of galaxies, which trace the large scale potential responsible for the CMB lensing, are visible at   mm and sub-mm wavelengths via the scattering of CMB photons by hot electrons (Sunyaev-Zel'dovich effect) and might therefore contaminate the cross-correlation signal to some extent. However, the redshift range populated by galaxy clusters only marginally overlaps with the redshift distribution of our sources, so that this contamination is negligible.
TABLE 2 H-ATLAS galaxy linear bias and cross-correlation amplitude as determined using both separately and jointly the reconstructed galaxy auto-and cross-spectra in the different patches. The best-fit values and 1σ errors are evaluated respectively as 50 th , 16 th , and 84 th percentiles of the posterior distributions. Chi-squares are evaluated as where b bf and A bf are the best-fit values. Note that posterior distributions of b and A obtained using only cross-correlation data are far from being Gaussian. Another systematic effect that can bias our measurement of the CMB convergence-galaxy cross-correlation is the leakage of cosmic infrared background (CIB) emission into the lensing map through the temperature maps used for the lensing estimation, as it correlates strongly with the CMB lensing signal (Planck Collaboration XVIII 2013). The 853 GHz map used as a Galactic dust template to clean the 143 and 217 GHz maps shows a particularly strong correlation.
H-ATLAS galaxies are well below the Planck detection limits and therefore are part of the CIB observed by Planck. Hence they are potentially contaminating the cross-correlation resulting in an enhanced observed signal. However, we could not find a feasible way to reliably quantify this possible contamination. We expect that with the next release of the Planck data, CMB lensing maps at different frequencies will become available. This will allow us to investigate the CIB leakage issue in more detail.
To check whether the CMB lensing maps in the H-ATLAS areas are affected by anomalously large noise we reconstructed the CMB convergence auto-power spectrum for each of the H-ATLAS patches as well as for four other patches, not covered by the H-ATLAS survey. The results of the analysis performed combining the five H-ATLAS patches show an excess of power for > 200-300 (Fig. 16). Considering the patches separately, we find that the main features of the CMB lensing power spectrum are recovered in the two largest patches, while the power spectrum in the three GAMA fields seems to be dominated by noise. The four non-H-ATLAS fields show a behaviour very similar to the H-ATLAS ones, indicating that the latter do not have anomalous noise properties.
To test the effect on the auto-and cross-spectra of errors on photometric redshift estimates we have redone the full analysis using the initial redshift distribution, dN/dz, i.e. the one represented by the dashed red line in Fig. 4. We get a slightly higher value of the crossspectrum amplitude (A = 1.70 +0.16 −0.17 ) and a somewhat lower value of the galaxy bias (b = 2.59 +0.11 −0.11 ). The reason for that is easily understood. As shown by Fig. 4, the convolution of the initial dN/dz with the smoothing kernel (representative of the uncertainties on estimated red-  shifts) results in a broadening of the distribution. This translates in a decrease of the expected amplitude for both the cross-and the auto-power spectra. Hence, in order to fit the same data, we need a higher value of the galaxy bias and, consequently, a lower value of the cross-spectrum amplitude A. Since the derived value of b is quite sensitive to the adopted redshift distribution, the agreement with other, independent determinations imply that our dN/dz cannot be badly off. Therefore it looks unlikely that the higher than expected value of A can be ascribed to a wrong estimate of dN/dz.
The tension between the estimated and the expected value of the amplitude A might be overrated because of an underestimate of the errors on the cross-spectrum. Thus, it is important to understand which term dominates the errors. The statistical error budget of the crosspower spectrum is plotted in Fig. 17. The auto-spectra contain a signal and a noise term asĈ XX L = C XX L +N XX L , so that the errors on the cross-spectra can be written as The first term represents the cosmic variance, the second one the pure noise and the remaining are mixed  signal-noise terms. As can be seen from Fig. 17, the main contribution to the C κg L variance is given by the noise only term. Moreover, the relative amplitude of the mixed terms is telling us that most of the error comes from the lensing noise. In order to reduce the errors of the reconstructed cross-spectrum, it is important to reach high-sensitivity in reconstructing the CMB lensing potential. This, of course, does not include the possible systematic errors discussed above.

SUMMARY AND CONCLUSIONS
We have presented the first measurement of the correlation between the lensing potential derived from the Planck data and a high−z (z ≥ 1.5) galaxy catalogue from the Herschel -ATLAS survey, the highest redshift sample for which the correlation between Planck CMB lensing and tracers of large-scale structure has been investigated so far. We have shown that the expected signal is remarkably strong, in spite of the small area covered by the H-ATLAS survey (about 1.4% of the sky), suggesting that cross-correlation measurements between CMB lensing maps and galaxy surveys can provide powerful constraints on evolution of density fluctuations, on the nature of the dark energy, and on properties of tracers of the matter distribution, provided that a good control of systematic errors for both data sets can be achieved.
The cross-correlation signal was detected with a significance of about 20 σ, substantially higher than expected. The reliability of this result was confirmed by several null tests. A joint analysis of the cross-spectrum and of the auto-spectrum of the galaxy density contrast yielded a galaxy bias parameter b = 2.80 +0.12 −0.11 , consistent with earlier estimates for H-ATLAS galaxies at similar redshifts. On the other hand, the amplitude of the cross-correlation was found to be a factor 1.62±0.16 higher than expected from the standard model and found by cross-correlation analyses with other tracers of the large-scale structure.
We have investigated possible reasons for the excess amplitude. Some of them, such as the redshift dependence of the bias parameter or the contamination by the Sunyaev-Zeldovich effect, were found to be negligible. Others, such as the higher magnification bias due to weak gravitational lensing or errors in the photometrically estimated redshifts, decrease the tension between the estimated amplitude and the theoretical prediction, but cannot fully explain the discrepancy.
We have also investigated the properties of the convergence maps. We find that in fact the convergence power spectrum shows a larger than expected dispersion around the theoretical prediction. This result is confirmed by an analysis of the power spectrum of the convergence map on other areas of the sky with geometry similar to that of the H-ATLAS regions. This suggests that the noise fluctuations in the observed patches may be higher than the average noise level on the sky. On the other hand, a similar noise properties were found analyzing four other, randomly chosen sky patches.
One possible source of systematic errors may be some residual contamination of convergence maps by unresolved point sources (Osborne et al. 2014;van Engelen et al. 2014), resulting in a stronger than expected correlation between the lensing convergence and the H-ATLAS high-z sources, which are unresolved by Planck.
We thank Karim Benabed for useful discussions and comments, Duncan Hanson and Michal Michalowski for a careful reading of the paper. FB would like to thank Simone Aiola, Matteo Calabrese and Giulio Fabbian for stimulating discussions. We gratefully acknowledge support from INAF PRIN 2012/2013 "Looking into the dust-obscured phase of galaxy formation through cosmic zoom lenses in the Herschel Astrophysical Terahertz Large Area Survey", and from ASI/INAF agreement 2014-024-R.0. FB acknowledges partial support from the INFN-INDARK initiative. LD, RJI and SM acknowledge support from the European Research Council (ERC) in the form of Advanced Investigator Program, COSMICISM. JGN acknowledges financial support from the Spanish CSIC for a JAE-DOC fellowship, co-funded by the European Social Fund. The work has been supported in part by the Spanish Ministerio de Ciencia e Innovacion, AYA2012-39475-C02-01, and Consolider-Ingenio 2010, CSD2010-00064, projects. The authors acknowledge the use of CAMB and HEALPix packages and of the Planck Legacy Archive (PLA). A.L. thanks SISSA for warm hospitality.