Exploring the Non-Gaussianity of the Cosmic Infrared Background and Its Weak Gravitational Lensing

Gravitational lensing deflects the paths of photons, altering the statistics of cosmic backgrounds and distorting their information content. We take the Cosmic Infrared Background (CIB), which provides plentiful information about galaxy formation and evolution, as an example to probe the effect of lensing on non-Gaussian statistics. Using the Websky simulations, we first quantify the non-Gaussianity of the CIB, revealing additional detail on top of its well-measured power spectrum. To achieve this, we use needlet-like multipole-band-filters to calculate the variance and higher-point correlations. Using our simulations, we show the 2-point, 3-point and 4-point spectra, and compare our calculated power spectra and bispectra to Planck values. We then lens the CIB, shell-by-shell with corresponding convergence maps, to capture the broad redshift extent of both the CIB and its lensing convergence. The lensing of the CIB changes the 3-point and 4-point functions by a few tens of percent at large scales, unlike with the power spectrum, which changes by less than two percent. We expand our analyses to encompass the full intensity probability distribution functions (PDFs) involving all n-point correlations as a function of scale. In particular, we use the relative entropy between lensed and unlensed PDFs to create a spectrum of templates that can allow estimation of lensing. The underlying CIB model is missing the important role of star-bursting, which we test by adding a stochastic log-normal term to the intensity distributions. The novel aspects of our filtering and lensing pipeline should prove useful for any radiant background, including line intensity maps.


INTRODUCTION
Gravitational lensing is the deflection of distant photons by intervening structure.Thus far, most attention of gravitational lensing has focussed on optical galaxies (Brainerd et al. 1996) and the Cosmic Microwave Background (CMB) (Lewis & Challinor 2006).In both cases, the changes to the distributions of the photons on the sky have been characterized, and observations of these changes have been used to extract critical science content on the evolution of structure (Hoekstra & Jain 2008).However, all extragalactic sources are gravitationally lensed by intervening large scale structure between the source and us.Upcoming intensity mapping surveys are gaining in sensitivity and extensive multi-line intensity mapping experiments are projected to probe more than 80% of the volume of the observable universe (Kovetz et al. 2019).Lensing of intensity mapping has been shown to be a challenge to detect with current surveys, but next-generation surveys might be able to detect lensing signals.This seems especially true when cross-correlation with other low-redshift tracers is utilized, potentially providing tighter constraints on cos-★ E-mail: astjason@sas.upenn.edumological parameters (Zahn & Zaldarriaga 2006;Pourtsidou et al. 2016;Foreman et al. 2018).
One example of intensity mapping with a long history of theory and observations is the the Cosmic Infrared Background (CIB), the emission from dust radiating down-shifted UV/optical radiation in star-forming galaxies to lower frequencies ( < 1000 GHz).Although the emission is from galaxy-scale objects, the finite resolution of the instruments used to detect it blend much of the emission in what appears as a diffuse component, albeit with non-Gaussian CIB anisotropies.It is possible that a significant CIB can arise from high redshift if there are luminous sources and dust, but the bulk of the emission is expected over the prime galaxy forming range, taken here to be from  < 4.2.The CIB traces well the clustering of galaxies throughout the evolution of the late-time universe.Its statistics were first elucidated in Bond et al. (1986), were put in terms of a "halo" model, with a Poissonian shot noise term and a continuous clustering term in Bond et al. (1991).These works focussed on amplitudes and 2-point statistics, with a first CIB-source map shown in Bond (1990) and maps with the full statistics using the Peak Patch algorithm done in Bond & Myers (1993) and Bond (1993), and further developed in Bond (1996) in the 90s, and in Stein et al. (2020) in the 2010s.
See also Gispert et al. (2000).The model used for the dust emission in these early works is of the same form as that used in the current literature, in particular as utilized in Planck analyses, but are based on mass-ordered emission that does not fully take into account the local perturbations that can cause galaxies of low mass to burst up in star formation activity compared with more massive ones, which can sometimes be more quiescent.The CIB model we adopt in this work is akin to the Planck CIB model (Ade et al. 2014b), which was also utilized in the Stein et al. (2020) Websky approach to extragalactic background mapping.Webskys use response functions of halos to luminous emissions, the CIB being one, to light-cone halos found using the peak patch algorithm and second order Lagrangian perturbation theory.The Websky (Stein et al. 2018;Stein et al. 2020) CIB maps are the direct descendants of the Bond & Myers (1993) and Bond (1996) maps of the early 90s.Webskys also give fully correlated lensing maps to complement the CIB maps, and are the basic tools used in this paper.
As with the CMB anisotropies, the CIB 2-point correlation and related power spectrum has been studied extensively, providing a window to the modeling of galaxy clustering at a wide band of redshifts (Ade et al. 2014b).Unlike the primary CMB radiation, the underlying localized nature of the CIB sources means it is intrinsically non-Gaussian so its monopole (average) emission and power spectrum, although valuable, misses the non-Gaussian clumpiness of the CIB.It is for this reason that the CIB bispectrum, the harmonicspace equivalent of the real-space 3-point correlation function, has been measured by multiple teams, such as with Planck (Ade et al. 2014b), the South Pole Telescope (SPT; Crawford et al. 2014), and the Atacama Cosmology Telescope (ACT; Coulton et al. 2018).Also, the bispectrum and higher order polyspectra have been modeled analytically, and it has been shown that complementary information from these measures can constrain the halo model further than the power spectrum on its own (Lacasa et al. 2014;Pénin et al. 2014).
Here, we present a formalism of analysing higher order spectra of extragalactic backgrounds through (angular or multipole) band filtering, with an application to maps of the CIB.By passing a CIB map through a contiguous series of filters defined in harmonic space by ranges of ℓ's characterized by a band centre ℓ  and a band-width Δℓ then calculating higher order map-statistics such as the skewness and kurtosis in each ℓ-band, we can quantify polyspectra that encode (albeit reduced or projected) information on the full n-point spectra in a straightforward manner.In particular, we explicitly compute the (reduced) bispectra at Planck frequencies.This is similar in spirit to quantifying non-Gaussianity with statistics like the bispectrumrelated power spectrum and the skew-spectra presented in Munshi & Heavens (2010) and Munshi et al. (2013).This method of band filtering proves to be a simple path to the quantification of the CIB non-Gaussianity.
Our main target in this paper is to develop a method to lens any radiation background that is made up of localized sources distributed over a broad redshift range, using the lensing associated with all mass structures below the object's redshift.Using this, we can consider the impact of gravitational lensing on CIB statistics.Although analytic calculations suggest that the gravitational lensing of the CIB does not change the CIB power spectrum substantially (Schaan et al. 2018), we show using our Websky simulations that its non-Gaussianity can be affected considerably.Prior to this work, weak lensing of the CIB has been investigated by Schaan et al. (2018), Mishra &Schaan (2019), andFeng &Holder (2019).Schaan et al. (2018) adapted the quadratic estimator (Hu & Okamoto 2002) to reconstruct the lensing of the CIB, and probed how the non-Gaussian nature of the CIB as well as its broad redshift extent biases lensing reconstruction.Mishra & Schaan (2019) quantified how much the lensing of foregrounds such as the CIB biases CMB lensing estimators, finding a small but potentially non-negligible effect.Feng & Holder (2019) found by comparing cross correlations between CIB lensing reconstructions and tracers such as the CIB at several frequencies and the CMB lensing potential for simulations and Planck data, that Planck CIB measurements contain excess non-Gaussianity consistent with CIB lensing.In this paper, we explicitly quantify the change in the CIB 3 and 4-point functions due to lensing relative to the unlensed CIB statistics.
Similar to the works of Schaan et al. (2018) and Mishra & Schaan (2019), our analysis pipeline can be adapted for various other non-Gaussian radiation fields associated source emissions, such as 21 cm, Lyman-, CO, CII, and other mm-wave intensity fields, to uncover further details of the 3-D structural evolution they trace.
An expansion to the Websky suite of extragalactic simulations1 (Stein et al. 2018;Stein et al. 2020) is used here to first show that the simulations capture the nearly-equilateral bispectra sufficiently well compared to the measured Planck values.We then investigate the change in CIB non-Gaussianity due to lensing, accounting for the fact that the CIB sources at different redshifts are lensed by different lensing convergences.Our simulation results suggest that weak gravitational lensing can change the 3-point and 4-point functions of the CIB to measurable degrees.This is potentially significant for observations, since the observed CIB necessarily includes gravitational lensing, while lensing has been ignored in theory predictions of the CIB bispectrum thus far.
The content of this paper is as follows.In Section 2, we describe the simulations used in our study, as well as the method we developed to accurately lens the CIB.In Section 3, we present our formalism for quantifying the n-point statistics.We then compare our simulation power spectra and bispectra with experimental data in Section 4, and discuss the redshift contribution to the CIB statistics in Section 5. We present our results on the effect of lensing on the CIB statistics in Section 6.In Section 7, we provide an alternative way to probe the effect of CIB lensing that includes additional non-Gaussian information using relative entropy.Next, we investigate stochastic effects on CIB statistics in Section 8. We wrap up with a conclusion and discussion in Section 9.

LENSED CIB SIMULATIONS
Our methods for performing the lensing of simulated CIB maps differ from lensing of the CMB in several aspects.First, while the CMB comes from a single redshift, the CIB comes from a broad range of redshifts, and the structures hosting CIB galaxies at one redshift will lens galaxies at a higher redshift.Second, while the CMB is a diffuse, smooth field, the CIB consists of many individual sources, and we must take care to treat the flux density from each galaxy appropriately, particularly since we aim to quantify non-Gaussian properties of the CIB both before and after lensing is performed.In the first part of this section we briefly describe the Websky simulations (Stein et al. 2020) our study is based on; we refer the reader to that paper for an in depth description.In the second part we describe how we calculate the lensed CIB signal.
The Websky simulations are based on the second order Lagrangian Perturbation Theory (2LPT) approach to nonlinear evolution of largescale structure, and follow the mass-Peak Patch approach (Bond & Myers 1996;Stein et al. 2018) to find regions that collapse into halos.First, peaks of density smoothed on a range of scales are flagged as potential halo candidates.For each such candidate, the 2LPT displacement field is then evaluated within the framework of the ellipsoidal collapse model to find regions that will gravitationally collapse into haloes.Overlapping halo candidates are either merged, or both kept with adjusted masses.
Each halo is modelled as a central subhalo with a central galaxy, for massive halos, with a number of satellite subhaloes that each hosts a single satellite galaxy; the number of subhaloes is generated using a subhalo mass function from Jiang & van den Bosch (2014).In the Websky simulation, there are ∼ 9×10 8 halos hosting ∼ 5×10 9 CIB galaxies.The flux density of each such galaxy is then calculated according to a CIB halo model (Shang et al. 2012), with parameters from Viero et al. (2013).In this model, the flux density of the galaxy depends only on the mass of its subhalo and its redshift.The key formulations of the model are given below as in Stein et al. (2020): The rest-frame spectral energy distrbution (SED) of a source is given by: where  is the frequency of the observation,  is the (sub)halo mass, and  is the redshift.]. ( eff is where the specific IR emissivity peaks, and  / describes the range of halo masses that produce the luminosity.log( eff / ⊙ ) = 12.3 and 2 / = 0.3 are used from model 2 of Viero et al. (2013). 0 is an overall frequency-independent normalization factor used to scale all galaxy flux densities to match the Planck 545 GHz CIB power spectrum measurements (Ade et al. 2014b) at the angular scale ℓ = 500.
This model describes a minimal parameter set that nevertheless provides a reasonable fit to the Planck and Herschel power spectra.One limitation is that every galaxy at a given redshift is assumed to have the same SED, regardless of its detailed properties, including its age, merger history, and environment.In principle, it could be possible to add more parameters, but these new parameters would be largely unconstrained by current CIB data.
To calculate a map of the unlensed CIB signal, we project galaxies onto a NSIDE = 4096 Healpix grid and add their flux densities within each pixel.The first panel of Figure 1 shows the CIB intensity in each Δ = 0.2 shell for the Websky simulations.The CIB intensity peaks around  = 1.4 to 1.6, which agrees with the CIB models in literature like Schmidt et al. (2015); Béthermin et al. (2012); Pullen et al. (2018).There is some uncertainty in where the CIB intensity actually peaks, as some models predict the peak to be at a lower redshift (Schmidt et al. 2015), but we note that it is the star-formation rate density that we expect to peak at 2 <  < 3 (Ade et al. 2014b).We have run our analysis on a model where the CIB intensity peaks around  = 1.1 and found that our results are not radically altered.
Although the Websky CIB maps, which go up to  = 4.2 as seen in Figure 1, capture the essence of the CIB, we note that the CIB extends further than  > 4.2 in reality.
Investigating the lensing of the CIB can be a complicated matter, especially when compared to CMB lensing.While the statistics of the unlensed CMB are very well known as it is essentially a Gaussian random field, the unlensed CIB itself is non-Gaussian and its statistics are not as well understood.Not only has the CMB power spectra been measured quite precisely (most recently, Aghanim et al. 2020a), the gravitational lensing effect on the CMB has been detected, and the matter between the observer and the CMB has been reconstructed to high precision (Aghanim et al. 2020b), utilizing quadratic estimators like the one given in Okamoto & Hu (2003).As mentioned in Section 1, the CIB power spectra and bispectra have been measured by a number of experiments, but the bispectrum measurements tend to have large uncertainties and its lensing reconstruction is in relatively early stages of development.Schaan et al. (2018) discuss methods to mitigate biases when reconstructing the lensing mass fluctuations from maps of the CIB.In particular, challenges arise both because of the intrinsic non-Gaussianity of the CIB itself, as well as the fact that there is redshift overlap between the lensing mass and the emissive sources.Such complications provide a motivation for our study using simulations.
Unlike the CMB, which is sourced at a narrow range of redshifts around  ∼ 1100 at the surface of last scattering, CIB sources are spread across a wide range of redshifts between  = 0 and  = 4.5.Using the full 3-D information from the Websky simulation, we generate a lensing convergence map for each redshift slice of CIB galaxies.This is a more involved treatment than that of Schaan et al. (2018), who took all the lensing matter to be at effectively a single redshift.Specifically, we split the galaxies into 21 redshift shells of width Δ = 0.2, and each galaxy in the -th redshift shell is lensed by all the matter in the first  − 1 shells 2 .This way we correctly account for the time evolution of the lensing effects and the fact that depending on the situation an individual galaxy can count as source of either CIB signal or lensing.To further simplify the calculation, when a galaxy acts as a source of CIB signal, we assume it is located at the central redshift of its redshift shell.This allows us to pre-calculate the lensing convergence for each redshift shell.
For each redshift shell, the lensing calculation starts with obtaining a NSIDE = 4096 map of the lensing convergence , obtained by integrating an appropriately weighted density field along the line of sight of the central positions of the individual map pixels.This matches what was done for the CMB lensing field in the released version of the Websky simulation, which was similarly calculated but for a source at  src = 1100.Within each halo, the density field is modeled as a NFW (Navarro-Frenk-White;Navarro et al. 1997) profile (Zhao 1996) and outside the haloes it is obtained from the 2LPT calculations.Following Stein et al. (2020), we also include a "field" component, representing the lensing matter that is in halos too small to be resolved by the simulation.
For the CMB, given the lensing potential map  and map of the unlensed CMB, the lensed map can be obtained by evaluating the CMB at the deflected positions given by ∇, using a pixel-based interpolation scheme.This is the approach that is codified in the commonly used and publicly available code Lenspix (Lewis 2005); it is also the approach that was taken for the lensed CMB map released with  src  being the midpoint of each redshift shell in -space, typically peak around half of its extent, although they display some skew towards where the CIB intensity is highest, especially as we integrate over more redshifts.The RMS deflection steeply increases at first from 0.34 ′ for the first lensed shell, up to 1.7 ′ for the last shell.
by the Websky team (Stein et al. 2020), in that case using the methods in the public pixell library3 .These interpolation-based codes have been shown to work accurately for the lensing of a diffuse field like the primary CMB that does not have significant fluctuations on scales near the pixel size.The primary CMB has this property thanks to its very red power spectrum, further enhanced by the washing out of structures on scales of a few arcminutes by diffusion damping at the last-scattering surface.However, since CIB galaxies are unresolved on the arcminute scales of the maps we consider, they appear as point sources the size of a pixel and Lenspix cannot be accurately used for this application.The smoothing of the signal map on small scales, which is an inherent part of Lenspix, leads to several artifacts such as nonphysical negative intensities we must avoid.Therefore, we developed our own lensing pipeline.
Rather than attempt to remap the given unlensed map with the deflection vectors, we opted to use the source galaxy catalogs and create entirely new CIB maps that include the effects of lensing.Given the convergence map    for a given source redshift   , which includes all the matter in the shells between the galaxy and the observer, we determine the lensed CIB intensity and position of each individual galaxy as follows: (i) We deflect the galaxy's angular position by an amount given by ì  = ∇   , with the lensing potential    given in terms of    by In terms of the 3D gravitational potential Φ( , ) at radial comoving distance  and conformal lookback time , we have More specifically, using spherical trigonometry, we invert Equations A15 and A16 of Lewis (2005), cos A patch of unlensed CIB and its corresponding lensing convergence for each redshift shell Δ = 0.2.In our simulations, each unlensed CIB shell is lensed by a convergence shell to create lensed CIB shells, which are then summed up to produce the total lensed CIB map.This method mitigates the 'self-lensing' effect substantially.Note that the CIB intensity visibly thins out by  = 3, while the integrated lensing convergence becomes brighter at higher redshifts.
and sin Δ = sin  sin  sin  ′ , where the lensed position is (, ) and the unlensed ( ′ ,  + Δ) on a sphere, with the deflection vector given as d ≡ ∇   =   e  +   e  =  cos e  +  sin e  .Inverting the equations above and reversing the deflection vector allows us to compute the lensed position of the galaxy from the unlensed position and the  corresponding to the unlensed pixel location of the galaxy following the Born approximation.
(ii) Since the deflection we have described will not change the apparent brightness of the galaxy, which is assumed to be a point source, we then multiply the galaxy's flux density by the magnification factor ) +    ,12 (e.g., Bartelmann & Schneider 2001).
Once we have the deflected and magnified galaxy positions, we combine lensed flux densities of all galaxies in a redshift bin in each pixel to obtain the total lensed CIB intensity from one redshift bin.We then add contributions from all redshifts to obtain the full lensed CIB map.
To account for the finite pixel size, we smooth each such  map with a Gaussian beam of  = ( √ 3 SIDE ) −1 , which roughly corresponds to the effective radius of a pixel.The unsmoothed  maps contain a substantial number of pixels where  > 0.1 (some are even larger than 0.3), which no longer lies in the weak-lensing regime.By smoothing the  map, we are able to reduce areas with large  values considerably.Smoothing the  maps suppress the  power spectra by about 50% at ℓ = 6000.Below, we will address how this choice impacts our results.
Because there still are non-negligible numbers of pixels in the  maps where the pixel values are larger than 0.1 even after the smoothing, we found that we must use the exact magnification factor [(1 −  2 ) −  2 ] −1 , rather than the often-used approximation (1 + 2).We use the transformation between  ℓ 's and  ℓ 's given in Eq. ( 11) of Jeffrey et al. (2021) as well as healpy's alm2map_spin() function to generate  1 and  2 maps.The shear factor that enters in the magnification is given by We note that our methodology of lensing each CIB shell with its corresponding  shell does not account for the fact that the halos and "field" components used to generate the lensing convergences get increasingly lensed as we move out to farther redshifts, in an effect known as lens coupling.It has been shown that incorporating this effect, along with other post-Born effects, suppresses the squeezed bispectrum and enhances the equilateral bispectrum of , especially at  > 1, although the impact on the  power spectrum is minimal (Pratten & Lewis 2016;Fabbian et al. 2019).Using lensing convergences not including some of these higher-order terms may impact our analysis, but we expect that not entirely capturing the non-Gaussianity of the lensing convergences does not significantly impact our results since the change in the CIB non-Gaussianity due to lensing is primarily caused by large scales of the lensing convergence maps.One potential way to incorporate some of the post-Born effects would be to lens each of the halo + field shell with its corresponding  before projecting it to the  shells, but this is beyond the scope of this work.
In Figure 1, we show the total CIB intensity in each redshift shell, their corresponding lensing kernels, and the RMS deflection angle.It can be seen that the lensing kernels and CIB intensity distributions are overlapping and trace the same matter fluctuations.The RMS deflection (bottom panel of Figure 1) in each redshift shell is calculated to range from 0.34 arcminutes for the very first lensed shell (0.2 <  < 0.4), to 1.6 arcminutes for the very last shell (4.0 <  < 4.2).The increase in the deflection is initally steep, but becomes more steady starting around  = 1, with the RMS deflection being about 1.3 arcminutes at  = 2. Overall, the mean RMS deflection for the CIB is around half of the amount (2.7 arcminutes RMS) that the CMB is deflected (Lewis & Challinor 2006).In Figure 2, we show all of the CIB shells as well as the corresponding  shell that "lenses" each CIB shell with their redshifts.It can be seen that the  shells are highly correlated, due to the overlapping lensing kernels, while the CIB shells are independent.Figure 3 shows a 0.5 • × 0.5 • patch of the CIB at 545 GHz for the redshift shell centered on  = 1.1.We show the unlensed, deflected (no magnification) and fully lensed signal.One can clearly see the downwards shift due to the deflection and the change in brightness after the magnification.In Figure 4, we display 1.5 • × 1.5 • patches of the  2 and  2 maps.We note that the  2 shows the clumps of matter inside the halos, while  2 shows the regions that surround the halos to have relatively high  2 values.

N-POINT CORRELATION FUNCTIONS WITH HARMONIC BAND FILTERS
Our main tool to study the effect of lensing on statistics of the CIB, especially its non-Gaussianity, will be the evaluation of variance, skewness and kurtosis of maps band limited to a certain range of spherical harmonic coefficients ℓ.As we show here, the former two are directly related to the power spectrum and equilateral bispectrum of the CIB map.Given a CIB map   at frequency  and its expansion into spherical harmonic coefficients we can define a band-filtered CIB map by only considering the  ℓ 's within a particular range of ℓ, where ℓ  is the center of an ℓ-band, We use Δℓ = 640 except for Section 4, where we compare with Planck's experimental data and use their binning of Δℓ = 128.
Given filtered maps, we can calculate their variance, skewness and kurtosis as where the average ⟨.⟩ is an average over map pixels, and Δ ℓ  () has the mean of the map subtracted.The subtraction in the last line ensures that  ℓ  4 is zero for a Gaussian random field.Notice alternative definitions can be found in the literature (e.g., Ben-David et al. 2015).
We checked that the sharp edges of the top-hat filters do not cause significant anomalies by running our maps through top-hat filters apodized at the edges with a sine-function as given in Eq. ( 9), where the (+) sign is for the left edge of the filter while the (−) sign is for the right edge of the ℓ band,  edge is either edge of the top-hat filter, and Δ is the width of apodization.With bandwidth Δℓ = 128 (Planck bins) and Δ = 5 to 15, we confirmed that the results are not qualitatively changed.
It is straightforward to show that the variance of the filtered map is related to the map's power spectrum through This equation can be inverted to give a prescription for calculating an estimate for the power spectrum in the bands, As usual, we define the angular bispectrum its angle-average as and the reduced bispectrum  ℓℓ ′ ℓ ′′ using where the matrices are Wigner-3  symbols.With these definitions, we can relate the skewness of the filtered map to  as follows (Komatsu & Spergel 2001): where We can invert Eq. ( 15) to evaluate an estimate of the equilateral bispectrum from the skewness of a band-filtered map as bℓ This is analogous to the equilateral binned bispectrum estimator in Bucher et al. (2016).We note that when ℓ  ≫ Δℓ, Eq. ( 17) reduces to: bℓ where we used the approximation for the Wigner-3  symbol (Bhattacharya et al. 2012) that if  = ℓ + ℓ ′ + ℓ ′′ is even and zero for odd .
The methods described above are computationally inexpensive  2 ) for a (1.5 • × 1.5 • ) patch of sky spanning up to  = 1.0 (so these are from the  and  shells that lens the CIB shell shown in Figure 3).The clumps of matter are evident from the bright pixels in the  2 map while the  2 closely traces the  2 but the bright regions are slightly different.The cosmic web structure can be seen more clearly in the  2 map.when using full-sky maps at NSIDE = 4096, with each run of passing a map through a series of top-hat filters and calculating  ℓℓℓ from Eq. ( 17) taking about 8 minutes on one node of the Niagara cluster4 .
As values of the CIB trispectrum are yet to be published to our knowledge, we use the kurtosis divided by ℓ 2  (ℓ −2   ℓ  4 ), which is proportional to the equilateral trispectrum at large ℓ  , as a proxy for how strongly gravitational lensing affects CIB 4-point functions.

COMPARISON WITH EXPERIMENTAL DATA
A number of surveys, including the Balloon-borne Large Aperture Submillimeter Telescope (BLAST; Viero et al. 2009), Herschel/SPIRE (Amblard et al. 2011), ACT (Dunkley et al. 2011), Planck (Ade et al. 2014b) and SPT (Hall et al. 2010;Reichardt et al. 2021;Crawford et al. 2014), have measured the CIB power spectra and bispectra.In this section, we compare the statistics of unlensed CIB maps obtained from Websky with some of these experimental results.
The Planck team measured the CIB power spectra up to ℓ = 2000 and bispectra up to ℓ = 800.We focus on the three frequencies where Planck has both power spectrum and bispectrum measurements ( 217Figure 5. Statistics of the unlensed CIB maps from the Websky simulation at the three Planck frequencies; top (blue) is 545 GHz, middle (green) is 353 GHz, and bottom (red) is 217 GHz.As we go to higher-order statistics, the Poisson regime becomes more evident as the spectra flatten out at ℓ > 1000.We note that the Websky bispectra are mostly within Planck error bars even though only the power spectra were fit to match those of Planck's.While we do not plot the error bars for Websky values as there is only one realization, one can estimate the level of uncertainty from the scatter especially for the bispectra and kurtosis.
GHz, 353 GHz and 545 GHz); the corresponding experimental data are shown in Figure 5.
Data analysis by Ade et al. (2014b) included masking bright sources over a brightness threshold (225 mJy for 217 GHz, 315 mJy for 353 GHz and 350 mJy for 545 GHz).For a fair comparison, we thus also mask bright sources.Because the full analysis including incomplete sky coverage is rather involved, we decided to instead replace the pixels brighter than the corresponding experimental cutoff with the mean of the CIB map.We believe this treatment is sufficient, due to the rareness of the very bright sources.When we instead replace the pixels above the experimental cutoff with double the map mean, our results are not significantly changed.
After treating the bright pixels in this manner and calculating the power spectra and bispectra of the CIB maps according to Eqs. ( 11) and ( 17), we get results shown in Figure 5.We also show  ℓ  4 results for completeness.
As in Stein et al. (2020), we see good agreement with the Planck CIB power spectra, though one has to keep in mind we scaled the amplitude of the 545 GHz power spectrum to agree with Planck at ℓ = 500.On the other hand, the comparison with Planck CIB bispectra is a nontrivial test of our model.While the Websky bispectra are systematically below the measured values, we generally agree within the error bars.We also see the transitions from the clustering regime to the Poisson noise dominated regime at high ℓ for all three polyspectra.Planck measurements primarily capture the clustering regime, especially for the bispectrum, while the Websky simulations also yield predictions for smaller scales.
To check the high-ℓ limit of our bispectrum calculations, we compare with the SPT results (Crawford et al. 2014).In the limit of high ℓ, the statistics are dominated by the 1-halo contribution and the bispectrum converges to a constant (see e.g. Figure 5).For the SPT flux density cut (flux cut from hereon) of 22 mJy, this constant is measured to be (Crawford et al. 2014) This is in good agreement with the Websky-based theoretical expectation, given by the weighted sum of the third power of the galaxy flux densities, where   are flux densities of the Websky galaxies dimmer than the SPT cutoff and  pix is the number of pixels in a NSIDE 4096 map (12 × 4096 2 ).However, we note that the Websky  ℓ value at ℓ = 2940 is about 40.7 Jy 2 /sr for 217 GHz while the SPT value is about (26.7 ± 0.7) Jy 2 /sr for 220 GHz.This mismatch for the  ℓ implies that the similar  Poisson values between Websky and SPT do not necessarily signify that the Websky CIB model explains SPT measurements well.A possible reason for this is briefly discussed in Stein et al. (2020); the Websky simulations use the Planck CIB model and the CIB contribution from halos smaller than ∼ 1.2 × 10 12  ⊙ is not included.

REDSHIFT ORIGIN OF CIB 𝑁-POINT FUNCTIONS
The CIB power spectrum is dominated by the redshifts around the CIB intensity peaks and lower.For the CIB at 545 GHz, as compiled by Schaan et al. (2018) using Béthermin et al. (2012); Schmidt et al. (2015); Pullen et al. (2018), the power spectrum is dominant over 0 <  < 2 with the CIB intensity peaking around 1 <  < 2 for various models (e.g., Ade et al. 2014a;Béthermin et al. 2017;Maniyar et al. 2018;Maniyar et al. 2021), whereas in the model Pénin et al. (2014) adopts, the CIB power spectrum is dominated by 0 <  < 3 with the CIB intensity peaking in the 1 <  < 3 range.In Figure 6 and Figure 7, we show the redshift accumulation of the CIB statistics as well as its ratio to the total CIB statistics.We confirm that the majority of the power spectrum comes from 0 <  < 2 for all scales.For the power spectrum, there is only a slight scale dependence; at ℓ < 2000, the contribution from  < 1.5 decreases at lower ℓ.For the bispectrum and kurtosis, we see a much larger contribution from  < 0.5, especially at ℓ > 1000, and even moreso for the kurtosis.At ℓ > 6000, more than 35% of the CIB bispectrum comes from  < 0.5 while more than half of the kurtosis comes from  < 0.5 for ℓ > 4000.While this is surprising, we attribute this to the fact that there are very bright galaxies at the closest redshifts which are not masked by the Planck flux cut.This is evident in Figure 8, where the Planck flux cut of 350 mJy at 545 GHz corresponds to roughly 5.6 MJy/sr.Because the bispectrum and kurtosis are higher powers of the fluctuations by definition, they are much more sensitive to these close-by bright galaxies.We also note that Planck has measurements with error bars up to around ℓ = 700.Pénin et al. (2014) also found that the CIB bispectra at frequencies lower than 857 GHz are dominated by low-redshift galaxies and Schaan et al. (2018) found that the CIB trispectrum is dominated by the lowest redshifts.As expected, most of the CIB statistics come from  < 2. For the power spectrum, the contribution from each Δ = 0.5 shell remain more or less constant throughout all scales.At subsequently higher-order statistics, the contributions from  < 0.5 increase considerably for ℓ > 4000.
Figure 7. Cumulative-redshift ratio of the CIB statistics relative to the total CIB, showing the same information as in Figure 6.It is clear that  < 0.5 contributes to the CIB bispectrum and kurtosis significantly at ℓ > 4000.However, higher redshifts are more important at low ℓ, especially for the bispectrum in the Planck measurement regime.Some ratios for kurtosis being larger than 1 at high ℓ is the result of imposing the same (Planck) flux cut for each cumulative shell.

LENSING OF CIB 𝑁-POINT FUNCTIONS
In Figure 9, we show fractional changes in the CIB power spectra, equilateral bispectra and kurtosis due to gravitational lensing.The changes in the power spectrum are small and below 2% at all scales, in agreement with Schaan et al. (2018).This is because, as visible from Figure 5, the CIB power spectra are relatively smooth and featureless.Compare this with the CMB, where the significant peak structure leads to up to ∼ 5% changes in power spectra due to gravitational lensing.We do not see any distinct difference between the three frequencies, with the peak effect around ℓ ∼ 2000 possibly related to the transition from a 2-halo to 1-halo dominated clustering regime.
The effects of gravitational lensing on the CIB bispectrum are more pronounced.The most significant effect is an increase in the bispectrum by about 15% at the largest scales, with the lensing influence gradually dropping to small values at ℓ ∼ 6000.The effects seem to be larger at higher frequencies.Overall, the lensing effects put our values closer to the Planck measurements, evident in Figure 10.
The situation is similar for the 4-point function, with lensing increasing the kurtosis of the CIB map by about 30% at the largest scales and the lensing importance dropping as we go to smaller scales.This time, the map frequency seems to matter much more, with the 217 GHz kurtosis showing relatively small (< 10%) lensing effects above ℓ = 1500 and 545 GHz map still showing ∼ 15% effects at ℓ = 5000.This however, may be attributed to the fact that the very first CIB shell spanning 0 <  < 0.2, where there is a significant contribution for the kurtosis at high ℓ, is not lensed in our analysis, together with the relatively high Planck flux cuts values for lower frequencies.In Figure 11, we investigate the sensitivity of these results to the strongly lensed regions.In the top panel we show results when we approximate the factor by which the galaxy flux densities are magnified as [(1 − ) 2 −  2 ] −1 ≈ 1 + 2 and find that regions where this approximation breaks down play crucial roles in the large lensing effects we see above.We see that including the  2 term has a smaller effect, impacting mostly large scales for the bispectrum and 4-point function.In the bottom we show results with no smoothing of the  map (cyan), or an alternative treatment of the high  tails where the values of  are capped at a certain maximum value  max (red and green).Again, we see sensitivity to how exactly the regions with high magnification are treated.While we only show results for 545 GHz, other frequencies show qualitatively similar behaviors.
We also show how the choice of flux cuts can impact our analysis in Figure 12.For the lensed case, we impose the flux cut after the lensing, as opposed to lensing the flux cut-imposed unlensed CIB, as this is equivalent to masking bright point sources in any observations of the CIB (which are lensed).While the change in the power spectrum due to lensing is essentially the same even when the flux cut values are increased or decreased by a factor of 2, the change in the bispectrum and kurtosis are significantly altered.This is because a substantial fraction of the CIB bispectrum and kurtosis come from bright sources at  < 0.2 (as seen in Figure 6), and increasing or decreasing the flux cut alters the contribution of these bright sources to the higher-order statistics.Indeed, when the flux cut is lowered to half of Planck's, the change in the bispectrum and kurtosis increase noticeably, while the change in higher-order statistics decrease accordingly when the flux cut is increased to twice of Planck's.We additionally note that the effect is more pronounced for the kurtosis, a significant portion of which comes from  < 0.2 at ℓ > 4000.Because the inclusion or exclusion of the brightest galaxies in the very first shell 0 <  < 0.2 has a major impact on the change in non-Gaussian statistics due to lensing, we further show in Figure 13 the change in CIB statistics due to lensing with and without the first shell.As expected, we find that removing the unlensed first shell altogether significantly increases the change in the bispectrum and kurtosis, especially at high ℓ.Hence, a more accurate portrayal of the change in CIB non-Gaussianity induced by lensing could be to break down the very first shell into finer shells (e.g.0 <  < 0.05 and 0.05 <  < 0.2), and lens all but the closest galaxies, although this is beyond the scope of this paper.

RELATIVE ENTROPY OF ℓ-BAND PDFS
While the 3 and 4-point functions are useful descriptions of the unlensed and lensed CIB non-Gaussianity, they do not capture the full CIB non-Gaussianity.Higher order -point statistics can be determined using the same method, determining the order  connected components of the ℓ-band-PDF, the probability distribution function determined in each of our ℓ-filtered maps.Or we could work with the ℓ-band-PDF directly, as a function of the CIB intensity.In practice, this involves constructing ℓ-band-histograms.In Figure 14, we show the unlensed and lensed CIB histograms at 545 GHz as well as their difference for selected bands.Although we can barely see the difference between the unlensed and lensed CIB when the two are overplotted, the difference between the lensed and unlensed clearly shows an interesting oscillation structure, caused by the deflection and hence smearing of the CIB at low-intensity pixels and magnification at high-intensity pixels.
To quantify the effects of the lensing signature in the PDFs more (1− ) 2 − 2 for different treatments of the  maps; no treatment, setting a hard  max , and smoothing at a pixel level.We see that the magnification significantly increases all 3 statistics, while using the approximation underestimates the non-Gaussian statistics in particular, and that setting a  max or smoothing the  maps substantially lowers the statistics.Smoothing the  map reduces the change especially at high ℓ.Results shown in the rest of the paper correspond to the choices in the thick blue lines.Lowering the flux cut masks considerably more pixels in the very first redshift shell  < 0.2, which we do not lens in our analysis.This causes the contribution to CIB statistics from  > 0.2 to increase, which is why the total lensing effect increases, especially for the skewness and kurtosis.Similarly, raising the flux cut decreases the lensing effect as the contribution from  < 0.2 becomes more important.
explicitly, we use two relative entropy quantities, one intensive and one extensive: The integral  rel is the (negative) of the Kullback-Leibler "distance" or KL divergence between the unlensed and lensed CIB PDF,  UNLENSED and  LENSED .In the following we often refer to  rel as the unweighted relative entropy; it portrays well the far tail dif-ference of the lensed and unlensed flux distributions.The extensive  rel , which we refer to as the weighted relative entropy, damps the extreme tails by the action of the  LENSED at large intensity, hence is more focussed on the mean and its vicinity, including the variance, skewness and kurtosis that we have concentrated on so far.
With the weighted version, one has a choice of multiplying by the lensed or unlensed PDFs, destroying the anti-symmetry between lensed and unlensed.(The small difference in the two PDFs means that the asymmetry is actually also small.)The information in the tails is visually enhanced by the logarithm, allowing for effective exploration of intensity regions far from the mean.The PDF-complexity  can thus encode much more of the non-Gaussianity than the low order connected components.Further extension would consider the PDFs to vary over coarse-grained space, with PDF-PDF spatial correlation analyses.We do not develop that extension here.
If we are trying to use relative entropy in practice, then we would be taking the observed histogram as the fully-lensed one, and would rather want to compare it with a model distribution that is lensed, but with a smaller or larger preferred amplitude than the observed PDF gives.We use a differential amplitude  lens, CIB to characterize the variation of the observed from the theoretical lensed CIB: + ln  th, CIB ( lens, CIB ,   , ℓ  ) .
Generally  lens, CIB will have a non-zero mean and a variance about it.A perfect model given the data would have a mean of zero.
In the differential between fully-lensed and lensing with a  lens, CIB ≠ 0 amplitude, the observed PDF drops out.Similarly, the unlensed theoretical distribution drops out if we are interested in small  lens, CIB variations from fully lensed.If our target is the explicit non-Gaussian nature of the PDF, as is often the case in cosmological applications, the KL distance of the lensed CIB PDF and a Gaussianized PDF with the same mean and variance are compared, and will have  lens, CIB dependence in both PDFs, complicating the template when comparing with the observed, but also focussing attention on what is truly non-Gaussian.In all cases we define a differential CIB-lensing template A similar approach to lensing was used in the Arcminute Cosmology Bolometer Array Receiver (ACBAR) power spectrum analysis (Reichardt et al. 2009), with the template as the relative entropy of Gaussians with different power associated with different lensing amplitudes,  rel = 1 2 Trace ln  (0) −1 ( lens, CMB ), where  ( lens, CMB ) is the CMB intensity correlation matrix at lensing amplitude  lens, CMB , taken relative to the fully-lensed  lens, CMB = 0 entropy.The terminology of relative entropy was not used in those days, but it does apply to previous works such as ACBAR.As in ACBAR and Calabrese et al. (2008), and ubiquitous in all subsequent lensed CMB power spectrum work, a different multiplier was used as the lensing parameter:  L,CMB , a measure of lensing strength multiplying the projected 2D gravitational potential.
At the linear level the two amplitudes are proportional to each other.We could also adopt an  L,CIB parameterization here for the CIB, but the small differences are such that  lens, CIB is adequate.
Observationally we do not know the unlensed spectrum, although much effort goes into trying to delens to isolate the unlensed.Thus the traditional expansion about unity for either  lens or  L has the classic problem of the number one in this context not being well determined since it is far from zero.Hence our emphasis here is on  lens about zero.
The model CIB depends upon a number of parameters, which we denote by   , which include dust temperature, dust density, the slope of the emission, and all of their redshift dependances, and other parameters which may appear in future improved CIB models.The relative entropy can then be expanded in basis elements (linear templates)   for each   , as well as the lensing basis element: We can then compute a spectrum of  lens, CIB and   in ℓ-bands preferred by observed data.One can also combine all ℓ-bands together since when the model is good there should be no ℓ dependence of the model parameters.The templates   and  lens, CIB will change as the preferred parameter values change.Subsequent iteration results in best fit (final)  ,f , and the   dependence on  ,f effectively turn the iterated linear expansion into a nonlinear one,  rel ( ,  ,  lens, CIB,  ).
A measure of how well the converged  rel does relative to the observations is the integral KL divergence,  rel ( ,  ,  lens, CIB,  ).If the templates are similar in shape to each other, then the associated parameters have near-degeneracies.Experimental noise can also enhance degeneracies, an issue that has been addressed in Horlaville et al. (2023) for the case of [C II ] line-intensity mapping.Also, with current CIB data, the high-ℓ regime may not be used due to the lack of experimental measurements of non-Gaussian statistics at high-ℓ.Nonetheless relative entropy analysis of PDFs is a promising avenue.For example, it could assist future CIB lensing reconstruction studies.
In Figure 15, we show the unweighted and weighted relative entropies plotted against  new , which is defined as: where Δ  =   − Ī is the CIB intensity with the mean CIB intensity subtracted out and  UNLENSED is the standard deviation of the unlensed CIB map (which is very close to that of the lensed CIB map).
We  UNLENSED | as |Δ  | → ∞), whereas when |Δ  | goes to zero, the scale is linear in intensity fluctuation.We also show on top of the figures the number of standard deviations away from the map mean Ī .We overplot 4 of 5 logarithmic ℓ-bands from ℓ = 1 to ℓ = 8000, leaving out the first band as it spans only a few ℓ's.These curves can serve as the aforementioned templates for the corresponding ℓ-bands.
In both the top and bottom panels, we see that the curves for bands 2 (ℓ = 38 to 145), 3 (ℓ = 145 to 552), and 4 (ℓ = 552 to 2101) are quite similar, but not so much for band 5 (ℓ = 2101 to 8000).We can interpret this to mean that the lensing amplitude is consistent throughout bands 2, 3, and 4, and that the CIB intrinsic parameters assumed for the Websky simulations are plausible.The discrepancy between band 5 and the other three bands show that the smearing and magnification effects are more apparent at the smallest scales, as we expect.
The application of relative entropy of PDFs to the quantification of the amount of lensing present is done here with templates characterized by band-amplitudes.As elucidated earlier, a general and ambitious use of relative entropy of PDFs is to expand it in a set of templates, each with an amplitude allowing for deviation from the "standard" dust emission model parameters.In the CMB parlance, these could be called nuisance parameters that we wish to marginalize over to obtain the lensing amplitudes.However, these parameters are of great physical interest, and one could determine best parameters using the templates, correct the templates according to the new amplitudes obtained, iterating until the relative entropy approaches zero.Thus, not only would one show the CIB is lensed, which of course it is, but with a self-consistent spectrum one would also have a refined CIB emission model determined by the data.Though we have shown a plausible path to a full PDF analysis as a function of ℓ  scale, the task on real data will be daunting.

STOCHASTIC EFFECTS
In section Section 6, we saw that while the bispectra derived from the Websky simulation are mostly within the Planck experimental error bars, they are consistently somewhat below the measurements.In this section, we test whether the situation can be alleviated by adding stochastic changes to galaxy flux densities.
As mentioned in section 2, in the CIB halo model used by Websky, the galaxy luminosity depends only on the mass of the corresponding subhalo.In reality, galaxy formation is affected by various environmental effects (e.g., Hearin et al. 2016) and galaxy luminosities will be consequently altered.
As a simple model of these stochastic effects, we multiply the flux density of each galaxy by exp [N (0,  ln  )], where N (, ) is a random number drawn from a Gaussian distribution with mean  and variance  2 .The parameter  ln  allows us to control the importance of these effects, and we consider three values:  ln  = 0.3, 0.4, and 0.5.More involved models are possible,  ln  could be different for central and satellite galaxies for example, but we do not consider such generalizations here.
We show changes relative to the case without the stochastic effects in Figure 16 for data at 545 GHz.In the left panel, we see that adding these stochastic effects increases the small scale power spectrum of the CIB at fixed large scale power.At ℓ = 5000, we see about 3%/7%/12% increase with  ln  of 0.3/0.4/0.5.For the bispectrum we see similar patterns, with larger relative changes as we go to smaller and smaller scales.Effects are a bit larger than for the power spectrum, with about 7%/15%/25% increase at ℓ = 5000 with  ln  of 0.3/0.4/0.5.The kurtosis shows the most interesting pattern, with a decrease at large scales and an increase at small scales, with the transition happening around ℓ = 1500.The amplitude of the change seems to increase with  ln  and is generally limited to below 15%.
While inducing stochastic effects on the flux densities of CIB galaxies does not increase the bispectra at low enough ℓ to compensate for the difference between Websky and Planck bispectra values (which was measured up to ℓ ≈ 800), our analysis shows that the power spectrum and bispectrum increase substantially due to stochastic effects at high ℓ.Thus we expect that future experiments might be able to distinguish between different models of stochastic effects as they can measure the CIB power spectra and bispectra down to small scales at high precision.

CONCLUSION AND DISCUSSION
In this paper, we have provided a simple and convenient formalism for estimating the non-Gaussianity of sky maps such as the Cosmic Infrared Background.Our method involves filtering the maps to isolate particular ranges of angular scales, and then computing the moments of the pixel distribution of the resulting map.Using this formalism, we first showed that the Websky simulations capture the equilateral bispectra of the CIB well, mostly within Planck error bars (although slightly lower).This is remarkable because no higher-order moments, such as bispectra, were considered when constructing the Websky CIB model.We also computed the kurtosis spectrum.All three types of spectra show clustering signals on large scales and a Poisson, or shot-noise, form on small scales.
We then lensed the CIB using a deflection-then-magnification method, by splitting the CIB into redshift shells and lensing each shell with its corresponding lensing convergence.This method accounts for the fact that the CIB is broadly distributed over a range of redshifts.We found that gravitational lensing causes the CIB power spectra to increase by less than 2% throughout all scales, the bispectra by 10 to 20% at large scales, and the kurtosis by 25 to 40% at large scales.
We found that the change in non-Gaussian statistics are quite sensitive to several factors: • the treatment of the convergence maps: setting a maximum value or smoothing; • the chosen magnification factor: using the commonly-employed weak lensing approximation  ≈ (1 + 2), rather than the more accurate  ≈ 1 (1− ) 2 or  = 1 (1− ) 2 − 2 ; and • the treatment of bright, nearby sources: the chosen flux cut threshold and whether the closest redshift shell is included.
The first two points are relevant to the fact that our calculations rely on the weak-lensing approximation, in which the lensing deflections are assumed to be small.As our results show, regions where this approximation is not completely valid play a non-negligible role in generating CIB non-Gaussianities.However, a full ray-tracing study is beyond the scope of this initial investigation.
On top of the effect of lensing on higher-order statistics of the CIB, we also laid out a procedure for examining the full non-Gaussian information of CIB lensing using relative entropy.We provided an example of how relative entropy and its differential can be expanded using the derivatives of the differential relative entropy with respect to a parameter as basis templates.Iterating over different values of the parameter allows us to constrain the amplitude of lensing as well as the intrinsic CIB parameters.Furthermore, we explored how stochastic effects not included in the halo model could change CIB statistics, and found that the statistics increase at small scales.
There are several reasons why this study is important.First, our lensing pipeline can be adapted for any 3-D intensity fields, such as 21cm, Lyman-alpha, as well as other mm-wave intensity fields like CO.In addition, CIB non-Gaussianity complicates the detection of primordial non-Gaussianity (Hill 2018;Coulton et al. 2022), which would be a very strong sign for inflation, and measurements of CMB lensing (e.g., van Engelen et al. 2014;Osborne et al. 2014).To isolate these effects from intrinsic CIB non-Gaussianity, the fact that CIB lensing changes the CIB non-Gaussianity should be considered (e.g., Mishra & Schaan 2019).Lastly, CIB non-Gaussianity provides additional information on top of the power spectrum in probing galaxy formation and clustering.In order to understand the underlying physics of galaxy formation and clustering, one needs to consider the fact that any observed CIB non-Gaussianity is lensed.We believe CIB lensing should be considered for upcoming highprecision surveys like the Simons Observatory (Ade et al. 2019) and Cerro Chajnantor Atacama Telescope-prime (CCAT-prime; Stacey et al. 2018;Aravena et al. 2022).SO is expected to observe CIB galaxies up to  ∼ 4 at high resolution in the 280 GHz frequency band, and potentially even observe what are now galaxy clusters at earlier stages of its evolution.CCAT-prime will be able to probe higher frequencies than SO and resolve up to 40% of the CIB at 850 GHz and detect a fraction of galaxies (< 0.5%) at redshifts up to 6, providing more insight into star formation rates using individuallydetected galaxies at mid-to-high redshifts than ever before.While Planck would not have been able to distinguish the lensing from the bispectra of the CIB itself given its large error bars, not accounting for CIB lensing in such upcoming surveys could potentially cause a bias.Additionally, CIB lensing could play a role in constraining intrinsic CIB parameters through the use of lensing templates.

Figure 1 .
Figure 1.CIB intensity, the corresponding lensing potential kernels, as well as the RMS deflection for each redshift shell (Δ = 0.2) within our Websky CIB model.The CIB intensity increases from  = 0 to  = 1.4,peaks around  = 1.4 to 1.6, then decreases until it is almost non-existent by  = 4.The lensing kernels, defined as    src  = 3 2 Ω   2 0

Figure 3 .
Figure 3. HealPix maps of unlensed, deflected (no magnification), and lensed CIB for a small (0.5 • × 0.5 • ) patch of sky centered on  = 1.1.Here, a "lensed" galaxy has both been deflected and has had its flux density magnified appropriately.The arrows denote the direction and magnitude of deflection.The light circled patch (Unlensed and Deflected) and the dark circled patch (Deflected and Lensed) are the same small patch of sky emphasized.One can clearly see the deflection by comparing the Unlensed and Deflected, and the magnification effect by comparing the Deflected and Lensed.

Figure 4 .
Figure 4. HealPix maps of the convergence-squared ( 2 ) and shear-squared ( 2 =  2 1 +  22 ) for a (1.5 • × 1.5 • ) patch of sky spanning up to  = 1.0 (so these are from the  and  shells that lens the CIB shell shown in Figure3).The clumps of matter are evident from the bright pixels in the  2 map while the  2 closely traces the  2 but the bright regions are slightly different.The cosmic web structure can be seen more clearly in the  2 map.

Figure 6 .
Figure 6.Contribution to the CIB statistics by redshift (cumulative) at 545 GHz.The build-up for the power spectrum and bispectrum are clearly visible.As expected, most of the CIB statistics come from  < 2. For the power spectrum, the contribution from each Δ = 0.5 shell remain more or less constant throughout all scales.At subsequently higher-order statistics, the contributions from  < 0.5 increase considerably for ℓ > 4000.

Figure 8 .
Figure 8. Breakdown of the CIB intensity histogram by redshift (no flux cut imposed).The closest redshifts contain very bright pixels from nearby galaxies, while the CIB dims significantly at  > 3.

Figure 9 .
Figure9.The effect of gravitational lensing on CIB statistics using lensing convergence maps smoothed at the pixel level.While the power spectra are changed by less than 2% for all three frequencies, the bispectra and kurtosis change by 10 to 40% at low ℓ.The apparent discrepancy between frequencies for the bispectra and kurtosis arises due to the relatively high flux cut values for Planck at lower frequencies.

Figure 10 .
Figure 10.Comparison of unlensed and lensed Websky bispectra with Planck measurements.Points are computed at the same set of central ℓ values but are horizontally offset for clarity.The lensed values are slightly larger than unlensed values and hence closer to Planck values.

Figure 11 .
Figure 11.Effect of lensing with various methods on CIB statistics at 545 GHz, illustrating our choice of methodology.The top figures show the changes in the variance (left), skewness (middle), and kurtosis (right) due to lensing with pixel-level smoothed  maps for proper and approximate magnifications.The bottom figures show the changes with  = 1

Figure 12 .
Figure12.Effect of lensing on CIB statistics at 545 GHz with different flux cuts.Lowering the flux cut masks considerably more pixels in the very first redshift shell  < 0.2, which we do not lens in our analysis.This causes the contribution to CIB statistics from  > 0.2 to increase, which is why the total lensing effect increases, especially for the skewness and kurtosis.Similarly, raising the flux cut decreases the lensing effect as the contribution from  < 0.2 becomes more important.

Figure 13 .
Figure 13.Effect of lensing on CIB statistics at 545 GHz with and without the very first shell  < 0.2, which remains unlensed.The change due to lensing is significantly more pronounced when  < 0.2 is excluded, especially for the kurtosis.Because the skewness and kurtosis are dominated by  < 0.2 starting around ℓ > 2000, removing the first shell dramatically increases the lensing effect on those statistics.

Figure 14 .
Figure 14.Top: Histograms of the unlensed and lensed Websky CIB maps at 545 GHz Bottom: The difference between the lensed and unlensed histograms, note that the number of pixels are roughly 100 times smaller for the difference histogram.While the two distributions look almost identical in the top panel, the difference of the two histograms uncover an oscillation structure in the probability distribution.

Figure 15 .
Figure15.The full information content of CIB lensing shown using relative entropy (which are essentially log lensing templates) with smoothed ; the top panel captures the tails while the bottom panel captures the distributions near the peaks.We show the  new axis, as defined in Eq. (25), on the bottom and the number of standard deviations away from the mean CIB intensity on the top of each panel, and Δ  =   − Ī .

Figure 16 .
Figure 16.Effect of stochasticity on CIB statistics at 545 GHz.The maps were renormalized by multiplying  −  2 ln  /2 with  = [0.997,0.994, 0.991] for  ln  = 0.3, 0.4, 0.5 respectively such that the power spectra at ℓ = 500 match the power spectrum of the unlensed Websky CIB map at 545 GHz.Stochasticity increases n-point statistics by various amounts: the increase becomes larger at high ℓ for all three statistics, but most dramatically for the skewness.
Viero et al. (2013) distribution Θ[,   ] is a greybody at low frequencies we consider for the CIB in this paper, Θ(, ) ∝     ( d ()) where   is the Planck function and  = 1.6 depends on the physical nature of the dust.The effective dust temperature is given by  d ≡  0 (1 + )  with  0 = 20.7 and  = 0.2.Φ() = (1 + )  CIB gives the redshift dependent global normalization of the  - relation with  CIB = 2.4 and a plateau at  = 2 is imposed as inViero et al. (2013).
lens, CIB e lens, CIB + ∑︁      , e  (  , ℓ  ) ≡ [ rel /  ] at   =   −  ,current .(24) use this nonlinear remapping of flux density values for the x-axis since it allows us to visualize the relative entropy more effectively at both small and large |Δ  |.When |Δ  | ≫  UNLENSED , | new | approaches the logarithm of the number of standard deviations from the mean (| new | → | ln