Large-scale Maps of the Cosmic Infrared Background from Planck

TL;DR

This paper presents large-scale cosmic infrared background maps derived from Planck data, using HI4PI neutral hydrogen data to remove Galactic dust contamination, enabling studies of large-scale structure and galaxy evolution.

Contribution

We developed a method to produce unbiased, large-scale CIB maps from Planck data by effectively removing Galactic dust contamination using HI4PI data, and made these maps publicly available.

Findings

CIB maps cover approximately 25% of the sky.

Power spectrum constrained down to angular scales of about 70 degrees.

Data products are publicly accessible for further research.

Abstract

The cosmic infrared background (CIB) is a powerful probe of large-scale structure across a very large redshift range, and consists of unresolved redshifted infrared emission from dusty galaxies. It can be used to study the astrophysics of galaxies, the star formation history of the universe, and the connection between dark and luminous matter. It can furthermore be used as a tracer of the large-scale structure and thus assist in de-lensing of the cosmic microwave background. The major difficulty in its use lies in obtaining accurate and unbiased large-scale CIB images that are cleaned of the contamination by Galactic dust. We used data on neutral atomic hydrogen from the recently-released HI4PI Survey to create template maps of Galactic dust, allowing us to remove this component from the Planck intensity maps from 353 to 857 GHz for approximately $25\%$ of the sky. This allows us to constrain the CIB power spectrum down to $\ell\gtrsim 70$. We present these CIB maps and the various processing and validation steps that we have performed to ensure their quality, as well as a comparison with previous studies. All our data products are made publicly available at https://doi.org/10.7910/DVN/8A1SR3, thereby enabling the community to investigate a wide range of questions related to the universe's large-scale structure.

Figures

Tables

Table 1. Table 1: CIB cross correlation coefficients, averaged over the multipole range 150 < ℓ < 1000 150 ℓ 1000 150<\ell<1000

	353 GHz	545 GHz	857 GHz
353 GHz	1	$0.98\pm 0.01$	$0.91\pm 0.01$
545 GHz		1	$0.96\pm 0.01$
857 GHz			1

Equations

$${\kappa }_{\ell m}=\frac{1}{2}\ell (\ell +1){\varphi }_{\ell m}.$$

$$I_{\nu }(\alpha ,\delta )=\sum _i{\alpha }_{\nu }^i\cdot N_{\mathrm{H}\text{mathsc}\mathrm{i}}^i(\alpha ,\delta )+{\beta }_{\nu }$$

$$I_{\nu }(\alpha ,\delta )=\sum _{\mathrm{ch}}{ϵ}_{\nu }^{\mathrm{ch}}\cdot T_B^{\mathrm{ch}}(\alpha ,\delta )+{\beta }_{\nu }.$$

$$J(\theta )$$

$$L(f(X|\theta ),Y)$$

$$C(\theta )$$

$${\rho }_{\ell }^{{\nu }_1{\nu }_2}=\frac{C_{\ell }^{{\nu }_1{\nu }_2}}{\sqrt{C_{\ell }^{{\nu }_1{\nu }_1}{C}_{\ell }^{{\nu }_2{\nu }_2}}}$$

$${\left(\Delta C_{\ell }^{T\kappa }\right)}^2=\frac{{\left(C_{\ell }^{T\kappa }\right)}^2+\left(C_{\ell }^{TT}+N_{\ell }^{TT}\right)\cdot \left(C_{\ell }^{\kappa \kappa }+N_{\ell }^{\kappa \kappa }\right)}{(2\ell +1)f_{\mathrm{sky}}\Delta \ell }.$$

Full text

cosmology: large-scale structure of universe; infrared: diffuse background, ISM

††journal: ApJ

1 Introduction

††©2019 California Institute of Technology.

U.S. Government sponsorship acknowledged

The cosmic infrared background (CIB) consists of the integrated emission from unresolved dusty star forming galaxies (Puget et al., 1996; Gispert et al., 2000; Lagache et al., 2005; Dole et al., 2006). This emission stems from dust grains bound to these galaxies and heated by the ultraviolet emission from young stars. As a result, most of the CIB emission originates from the peak of the star formation epoch at $z=1-2$ and originates in galaxies living in dark matter (DM) halos with masses of $10^{11}$ to $10^{13}\,M_{\odot}$ (Béthermin et al., 2012; Schmidt et al., 2015). The CIB is an excellent tool with which to study the cosmic star formation history and the connection between dark and luminous matter.

Historically, the monopole of the CIB was first detected through measurements with the FIRAS and DIRBE instruments aboard the COBE satellite (Puget et al., 1996; Fixsen et al., 1998). It took another decade to also detect the anisotropies of the CIB (e.g. in IRAS data, Miville-Deschênes et al., 2002), which enabled the first measurement of the auto power spectrum and constraints on the bias of the CIB sources (Lagache et al., 2007). The latest milestone in understanding the CIB was a result of the Planck and Herschel missions (Pilbratt et al., 2010; Planck Collaboration, 2018 I). The combination of high resolution, high sensitivity, broad frequency coverage, and large area allowed for the extraction of unprecedented maps, and for constraints on a wide range of models (Planck Collaboration, 2011 XVIII, 2014 XXX; Serra et al., 2014; Planck Collaboration, 2016 XLVIII; Mak et al., 2017; Viero et al., 2019).

On top of being an excellent probe of star formation, the CIB is also a good tracer of the large-scale structures out to high redshift, and thus a good tracer of the lensing potential that affects the cosmic microwave background (CMB) (Planck Collaboration, 2014 XVIII). As a consequence, the CIB can be used to predict the lensing potential and to de-lens the CMB temperature map as demonstrated, for example, on Planck data (Larsen et al., 2016). Manzotti et al. (2017) used CIB data from Herschel to de-lens the CMB B-mode map from the South Pole Telescope (Carlstrom et al., 2011), thereby detecting for the first time a statistically significant de-lensing effect on the CMB B-mode. Forecasts conducted by Manzotti (2018) show that for next-generation CMB Stage-3 and Stage-4 experiments, the internal reconstruction of the lensing potential will be the dominant method to de-lens the CMB, but ancillary data will remain very helpful and serve as a very important and independent systematic check. Most notably, the CIB covers the lensing potential out to much higher redshift than what can be covered in current large-area galaxy surveys, albeit with much lower resolution. Lastly, Planck Collaboration (2018 VIII) combined the Planck lensing reconstruction at large scales, which was high signal-to-noise, with CIB information from the Planck generalized needlet internal linear combination (GNILC) maps (Planck Collaboration, 2016 XLVIII) at small scales to construct an optimal map of the lensing potential that combines the advantages of both tracers.

Large-area maps of the CIB are maps of a highly biased tracer of dark matter at $z=1-2$. It is thus a promising avenue to constrain the scale-dependent bias imprinted by primordial non-Gaussianity (e.g. Dalal et al., 2008; de Putter & Doré, 2017, and references therein). Tucci et al. (2016) use Fisher forecasts to demonstrate that, even in the presence of Galactic dust residuals, an uncertainty of $\sigma(f_{\rm NL})$ of approximately 3.5 can be obtained for sky fractions between 0.2 and 0.6. This result would be competitive with CMB bispectrum-based measurements (Planck Collaboration, 2016 XVII) and would rely on different physical scales. Current CIB-based measurements of this effect are limited by the available area of the CIB maps with low enough residual galactic dust. Our work aims at directly addressing this limitation.

In the present work, we extend the H i column density-based approach to be more general, less subjective to human decisions, and to cover a larger fraction of the sky. At the same time, we now also have access to newer data, both from Planck for the far-infrared (FIR) data, and from the HI4PI Survey (HI4PI Collaboration et al., 2016) for the H i data. We present the underlying data sets in Section 2, our methodology in Section 3, the estimate of the power spectrum in Section 4, the map-based results in Section 5, and the results for the power spectra in Section 6. We present the validation of the results in Section 7, and finally conclude in Section 8.

2 Data and Preprocessing

We briefly describe the individual data products used throughout this study and the preprocessing that was performed on the publicly available data.

2.1 H i data

The H i data are based on the recently published HI4PI Survey (HI4PI Collaboration et al., 2016). This survey merges data from the Effelsberg–Bonn H i Survey (EBHIS, Winkel et al., 2010; Kerp et al., 2011; Winkel et al., 2016a) and the Galactic All-Sky Survey (GASS, McClure-Griffiths et al., 2009; Kalberla et al., 2010; Kalberla & Haud, 2015) to create a full-sky database of Galactic atomic neutral hydrogen.

Compared to its predecessor, the Leiden/Argentine/Bonn Survey (LAB, Kalberla et al., 2005), it offers full spatial sampling (instead of beam-by-beam sampling), higher angular resolution ($16.1^{\prime}$ instead of $35^{\prime}$), and higher sensitivity.

Our preprocessing of the HI4PI data is described below. This data set is publicly available on CDS.222http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/594/A116

We merge the individual smaller cubes (each containing data for one HEALPix pixel at Nside 4) into one large HDF5 table. 2. 2.

The Magellanic System, especially the Magellanic Stream, is a major contaminant due to its low dust content despite high H i column densities. We follow the procedure detailed in A.2 of Planck Collaboration (2014 XXX) to mask this emission in the 3D cube. We use the Milky Way model described in Kalberla & Dedes (2008) and mask all emission with model brightness temperatures less than 60 mK and Magellanic coordinates $240^{\circ}<\lambda<30^{\circ}$ and $|\beta|<10^{\circ}$. 3. 3.

Instead of working with the full spectral resolution of the HI4PI Survey, we bin the data along the spectral axis. We choose a non-uniform binning to capture the highly complex emission close to the velocity of the local standard-of-rest $v_{\rm LSR}=0$ while using large bins at the higher velocities that contain little emission. This binning scheme is illustrated in Figure 1. We also exclude emission with $|v_{\rm LSR}|>90\,\rm km/s$ because of its very low dust-to-gas ratio (Planck Collaboration, 2011 XXIV; Lenz et al., 2016). This binning scheme is applied to the full sky, and does not differ from sightline to sightline.

2.2 Planck data

For the FIR data, we used the Public Release 3 (PR3) data from the Planck satellite. We work with the high-frequency data at 353, 545, and 857 GHz. For the highest frequencies, the dust emission from the Milky Way and from the CIB dominate over other processes such as CMB, synchrotron, or free-free emission. For the lower frequencies, the Planck bands have a strong CMB component. For consistency, we remove the CMB from all frequencies.

At 217 GHz, we found the strong CMB contribution to be a limiting factor in the analysis. The resulting CIB map at this frequency depends strongly on the exact details of the CMB correction, and is consequently not presented here.

Our full preprocessing pipeline of the Planck data is described below. All the data products were obtained from the Planck Legacy Archive333http://pla.esac.esa.int/pla/. For most of the analysis, we used the healpy444https://healpy.readthedocs.io/en/latest/ implementation (Zonca et al., 2019) of the HEALPix pixelization (Górski et al., 2005). We refer to the $353-857\;\rm GHz$ data as intensity maps. Some of the operations listed below require converting the maps back and forth from $K_{\rm CMB}$ to $\rm MJy/sr$. This is done by using the conversion factors given in Planck Collaboration (2016 X, Table 3).

We start with the Planck intensity map in temperature at frequency $\nu$, at full resolution, and HEALPix Nside of 2048 (e.g. HFI_SkyMap_545_2048_R3.00_full-oddring.fits, using the I_STOKES data and II_COV for the covariance). 2. 2.

We convert the map from nested ordering to ring ordering. 3. 3.

We subtract from the intensity maps the CIB monopole that was added by hand to the Planck maps at 353, 545, and 857 GHz, following the model of Béthermin et al. (2012). The values are given in Planck Collaboration (2018 III, Table 12). 4. 4.

Given the beam window function in the Planck Reduced Instrument Model (RIMO), we convolve the maps to the window function of the SMICA CMB map. This means converting the map $a_{\ell m}$, multiplying these by $B^{\rm kernel}_{\ell}=B^{\rm\;SMICA}_{\ell}/B^{I_{\nu}}_{\ell}$, and converting back from $a_{\ell m}$ to map. 5. 5.

With the intensity maps now at the same 5’ resolution as the CMB map, we subtract the SMICA CMB map from the intensity map. The maps need to be in units of $K_{\rm CMB}$ for this. 6. 6.

Lastly, we downgrade the resolution of the HEALPix grid from an Nside of 2048 to an Nside of 1024, which is sufficient to fully sample the data at this resolution. We only conduct the generation of the foreground dust model at Nside 1024. This model is subsequently upsampled to Nside 2048 for all following analyses. The final published CIB maps are also distributed at Nside 2048.

We note that the beams are slightly asymmetric, which will introduce a 1-2% error at $\ell\approx 1500$, which we consider as negligible here. This effect has been investigated in Section 4.3 of Planck Collaboration (2018 VI).

Furthermore, the subtraction of the SMICA CMB map could be biased due to the contamination of that CMB map with residual CIB and Galactic dust. An illustration of the order of magnitude of these effects is given in Planck Collaboration (2016 IX, Figure E.6), showing that the contamination is very minor, in particular for the higher frequencies.

To avoid noise bias when measuring angular power spectra, we are measuring cross spectra between ring maps. The so-called odd- and even-ring maps are generated using only the first or the second half of each pointing period. These maps were used by the Planck team to test internal consistency and to characterize the noise by taking the difference of these half-ring maps. The periods are typically 20 minutes long, and there are two half-ring maps at each frequency. Each of the two half-ring maps at each frequency is processed according to the steps outlined above. We also ran the entire analysis on half-mission maps rather than half-ring maps, and found no significant difference.

2.3 The Planck lensing map

In order to quantify the effectiveness of our Galactic dust cleaning procedure, we will use the publicly available Planck lensing convergence map (Planck Collaboration, 2018 VIII), which is known to be strongly correlated with the CIB (Planck Collaboration, 2014 XVIII). To cross-correlate this lensing map with our CIB map, we generate a lensing map at an Nside of 2048, based on the publicly available lensing $a_{\ell m}$ and lensing mask. We note that this map is already corrected for the response function, including the beam function and the pixel window function.

In addition, the lensing angular power spectra are often expressed in terms of the lensing potential $\phi$ and not the lensing convergence $\kappa$. The two are related through

[TABLE]

In this paper, we express all lensing-related quantities using the lensing convergence $\kappa$.

2.4 Masks

To compute an accurate and reliable angular power spectrum of CIB anisotropies, we need to mask different components that would otherwise bias the analysis. Figure 2 shows the final mask that is used for our analysis. There are several components that contribute to this effective mask.

First of all, we limit our analysis to the low H i column densities. Based on our experience in Planck Collaboration (2014 XXX) and Lenz et al. (2017), we choose $N_{{\mathrm{H\,\mathsc{i}}}}=2.5\times 10^{20}\,\rm cm^{-2}$ as the threshold for our baseline results shown here. In addition to that, we also use the Planck 20% Galactic plane mask, so the final mask is the intersection of the two. The complexity of the interstellar medium (ISM) close to the Galactic disk makes the CIB measurements in this region almost impossible. We further explore different cuts in H i column density in Section 7.1.2, and publish maps with different sky fractions and thus different levels of Galactic dust residuals.

Second, we employ the public Planck masks of extragalactic point sources, which are limited to a signal-to-noise ratio of 5. For further reference and for the specific flux limits, see (Planck Collaboration, 2016 XXVI). These masks are slightly different for each frequency, with the highest, 857 GHz channel containing the most point sources. Additionally, we intersect them with the effective mask of the SMICA CMB map.

Third, we find that a major source of contamination is molecular gas at high Galactic latitudes, often associated with molecular intermediate-velocity clouds (MIVCs, Magnani & Smith, 2010; Röhser et al., 2016). On top of that, the linear correlation between H i column density and FIR dust emission only holds in the absence of CO-dark molecular gas, which cannot be observed directly (Planck Collaboration, 2011 XIX). A census of these objects was performed in the work of Röhser et al. (2016), and we use a HEALPix mask of all these sources (T. Roehser, private communication). With the $N_{\mathrm{H\,\mathsc{i}}}$ threshold already in place, this mask only excludes an additional $21\;\rm deg^{2}$.

Lastly, we mask residual Galactic dust emission in our final CIB maps. In most cases, this residual emission results from dust that is associated not with H i, but instead with CO-bright or CO-dark molecular gas. For the column densities investigated here, the influence of optically thick H i is negligible (e.g. Lee et al., 2015). We provide further details on this procedure in Section 3.

2.4.1 Mask apodization

To compute the power spectra of our maps more reliably, we apodize these masks following the procedure for cosine apodization described in the NaMaster software (Alonso et al., 2019). We choose an apodization kernel with an FWHM of $15^{\prime}$ for this procedure. The final boolean mask and the apodized mask are presented in Figure 2.

3 Building new CIB maps

To disentangle the FIR emission from dust in the Milky Way and that from unresolved CIB galaxies, two different methods have been effectively used to date. We describe them in the following sections. For completeness, we note that the spectral information from the FIR signal alone cannot be utilized to separate Galactic and CIB dust, unlike what is done to separate the CMB component (e.g. Commander, Eriksen et al., 2008; Planck Collaboration, 2016 X). This is due to the fact that both dust components are very well described by a modified blackbody, and a component separation based on this information alone would thus be strongly degenerate.

3.1 Disentangling via the spatial structure

The different spatial structures of the CIB and Galactic dust can be used to disentangle the two. This has been demonstrated with an implementation of the GNILC method (Planck Collaboration, 2016 XLVIII) or by subtracting the foreground dust at the power spectrum level (Mak et al., 2017).

Using the differences in angular power spectra of these two components introduces a powerful source of information, but is also limited when the power spectrum of the fields that are reconstructed is not known to some degree. Mak et al. (2017) assume a parametric form for the CIB power spectrum, which has been previously derived through H i-based component separation (Planck Collaboration, 2014 XXX, hereafter P14). An approach based on the GNILC method (Planck Collaboration, 2016 XLVIII) has shown that Galactic dust and CIB emission can be disentangled for large parts of the sky ($\gtrsim 60\%$), albeit with the caveat of oversubtracting the CIB (see Maniyar et al., 2019, Appendix A).

3.2 Disentangling via the H i column density

A commonly used approach to remove the foreground FIR intensity is to model its value based on the H i column density (e.g. Planck Collaboration, 2011 XVIII, 2014 XXX). In this case we model the observed FIR intensity in the following way:

[TABLE]

Here, $I_{\nu}$ is the observed FIR intensity at frequency $\nu$, $\alpha^{i}_{\nu}$ is the dust emissivity per hydrogen nucleon, $N^{i}_{{\mathrm{H\,\mathsc{i}}}}$ is the H i column density, and $\beta_{\nu}$ is the zero-point. The subscript $\nu$ indicates that this equation is solved at each FIR frequency individually and the arguments $(\alpha,\delta)$ indicate the position on the sky. The index $i$ indicates that the spectroscopic H i data can be binned into several column density maps with different velocity ranges. This approach is motivated by the fact that the dust emissivity $\alpha^{i}_{\nu}$ is different for multiple components because of a different composition or illumination by the interstellar radiation field (ISRF). By accounting for this difference, the FIR intensity can be modeled more accurately. The offset $\beta_{\nu}$ is required to better capture the spatially varying zero-point of dust and gas. In particular, the warm ionized medium (e.g. Lagache et al., 1999) contributes to large-scale fluctuations of this zero-point.

Most commonly, we distinguish between low-velocity, intermediate-velocity, and high-velocity clouds (LVCs, IVCs, and HVCs; see the review by Putman et al., 2012, and references therein). Usually, the shape of the H i spectrum is used to define these velocity ranges. An illustration of this procedure can be found in Figure 1 of Planck Collaboration (2011 XXIV). In practice, the separation of the H i data cube into multiple column density maps is done manually and can be subjective. Due to these limits, it is difficult to apply this method over areas of thousands of square degrees.

The most recent applications of this approach, particularly in light of the Planck data, are demonstrated in Planck Collaboration (2011 XVIII) and Planck Collaboration (2014 XXX). Here, H i data from the GHIGLS survey (Martin et al., 2015), from GASS (Kalberla & Haud, 2015), and from EBHIS (Winkel et al., 2010) have been used to clean the FIR intensity maps of foreground dust.

Similar to many previous studies on the subject (e.g. Planck Collaboration, 2014 XXX, for the latest application of this approach on FIR data), we use the H i data as a template for Galactic dust. However, to extend the sky coverage we use a slightly more automated and general approach. In the following, we describe the details of our component separation technique.

We focus on two aspects. First, we describe the spatial selection of the regions where the Galactic dust and CIB emission will be separated. The challenge here is to account for the spatial variability of the dust-to-gas ratio, while at the same time preserving the large-scale CIB fluctuations. Second, we detail how we utilize the full three-dimensional H i data set, and how the information on the radial velocity of the gas helps to optimally construct a model of Galactic dust.

3.3 Selecting the spatial structures

The correlation between gas and dust varies across different physical environments and thus across the sky. Thus, we divide the sky into patches and perform the foreground modeling on one patch at a time, allowing us to turn this into an embarrassingly parallel problem. This procedure has already been proven to work very effectively in Planck Collaboration (2014 XXX) and in Planck Collaboration (2014 XVIII).

One caveat to this procedure is that a local determination of the dust-to-gas ratios, and through that the CIB intensity, removes large-scale fluctuations of the CIB that extend beyond the size of the patches. In the power spectrum domain, we observe this as a sharp drop at low $\ell$. This effect has also been observed in Planck Collaboration (2014 XVIII, Figure 7, right panel) and Schmidt et al. (2015). We show an example of this effect in Figure 4. There, we generate a simulated CIB map at 545 GHz and subtract the mean value in each patch. For the patches, we use the HEALPix pixelization and demonstrate the effect for different Nsides. In the following, we refer to these larger HEALPix patches as superpixels.

Due to this limitation, we set up our procedure as follows. The HEALPix pixels on which we locally perform the modeling are chosen to have an Nside of 16 (we further discuss this choice in Section 7.1.3), giving us 4096 full-resolution pixels in one patch. After applying the model that is described in more detail in the following section, we convolve the full-sky maps of model parameters with a Gaussian kernel with an FWHM of $3^{\circ}$, thereby creating continuous maps of dust-to-gas ratios for different radial velocities, and for the zero-point in the H i/FIR relation. A similar version of this procedure has already been applied by Planck Collaboration (2014 XXX) in their analysis of the large GASS field.

To generate the dust model, these continuous parameter maps are then combined with the H i data (Eq. 3), which are subject to this spatial filtering. We further discuss the implications of this in Section 7.1.3.

3.4 Using the full spectral H i information

A common practice is to manually set the velocity ranges of the different H i column density maps for a single field. There are three caveats to this procedure. First, the separation is manual and subject to human opinion and error. Second, the typical number of column density maps for one field ranges from one to three, thus offering no possibility to account for higher-resolution substructures. Third, the separation may not be well defined for a given field, allowing different interpretations of where the low/intermediate/high-velocity gas intersections are.

To overcome this manual step and to allow the model to be more universal, we implement a generalized linear model (GLM), similar to the analysis done in Lenz et al. (2016). This means that we assume different emissivities for different radial velocity bins by allowing each spectral channel to have an individual dust emissivity. Building up on Eq. (2), the GLM for the FIR intensity $I_{\nu}(\alpha,\delta)$ at sky coordinates $(\alpha,\delta)$ can be written as a weighted sum over all spectral channels:

[TABLE]

Here, $T_{B}^{\mathrm{ch}}(\alpha,\delta)$ denotes the H i brightness temperature in each spectral channel and $\epsilon^{\mathrm{ch}}_{\nu}$ is the emissivity for each individual channel. Similar to Eq. (2), $\beta_{\nu}$ is a constant offset.

The HI4PI data have a total of 933 spectral channels with a channel width of $1.3\,\rm km\,s^{-1}$ each, which gives us far more spectral resolution than required. By spectrally binning the H i as described in Section 2.1, we still have a sufficient number of spectral H i components while at the same time reducing the degeneracy of the model, reducing the computational cost, and increasing the sensitivity of the data.

Many of the H i spectral channels contain no significant emission and are masked when solving Eq. (3). We use the 3$\sigma_{\mathrm{RMS}}$ of the H i data as a clip level. This avoids that CIB fluctuations being falsely fit by using a large number of degrees of freedom, thereby creating an artificial, false image of the background component.

To further constrain our model and avoid chance correlation of H i and the CIB, we use lasso regularization for the individual GLM coefficients $\alpha^{\mathrm{ch}}_{\nu}$ (Tibshirani, 1996). Thereby, we break the degeneracy between neighboring, correlated channels and avoid overfitting. In practice, this means that we define our cost function $J(\theta)$ as the sum of a square loss function $L(f(X|\theta),Y)$ and the total sum of the coefficients $C(\theta)$:

[TABLE]

Here, $||\cdot||_{1}$ denotes the $L^{1}$-norm and $\gamma$ is the strength of the regularization. The second term ensures that the fit yields a sparse set of coefficients.

To scale the regularization strength $\gamma$, we use cross validation (e.g. Picard & Cook, 1984). This technique is commonly used to optimize hyperparameters such as $\gamma$, which are not directly evaluated on the data, but need to be determined before the actual fitting procedure. For this purpose, the image is split into $n$ parts (so-called folds). For the $n$-fold cross validation, the data are fitted on $n-1$ folds (training sample) and the quality of the fit is evaluated on the $n$th fold (test sample). Each of these folds is randomly selected, and they do not represent spatially coherent features. This procedure is repeated $n$ times, so each fold serves as test sample exactly once. Therefore, the GLM is not evaluated just once for an image, but instead hundreds of times to find an appropriate solution. For our purpose, we work with values of $n=3...5$ and find that the exact choice of $n$ does not affect the results significantly.

An illustration of the GLM approach is given in Figure 5 for a random superpixel. We show the normalized GLM coefficients, the mean spectrum for that superpixel, and the maximum spectrum. The latter helps to identify spatially small features in an otherwise low-signal region. We further demonstrate the result of this technique in Figure 6, where we show the FIR data, model, and derived CIB map.

To remove the residual emission resulting from CO-dark molecular gas, we repeat iteratively the procedure detailed above and evaluate the Gaussian distribution of the resulting CIB maps. Outliers of more than $3\sigma$ are masked and the modeling step is repeated until convergence is reached (see also Figure 12 of Planck Collaboration (2011 XXIV) for further reference).

Finally, we obtain the CIB maps by subtracting the dust models from the total intensity maps.

We note that the discrete nature of the FIR emissivities that is assumed in this approach is primarily a schematic description and does not necessarily yield an accurate value of the physical dust emissivity for each binned spectral channel of the H i data. This is acceptable, however, because our main goal here is to create a robust foreground model, free of residuals from Galactic dust.

4 $C_{\ell}$ estimation

To estimate the angular power spectra, based on the CIB and CMB lensing fields, we use PyMaster, the Python implementation of NaMaster555https://github.com/LSSTDESC/NaMaster (Alonso et al., 2019). When estimating the $C_{\ell}$ at large scales, incomplete sky coverage leads to mode–mode coupling, often approximately described in the matrix $M_{\ell\ell^{\prime}}$ (Hivon et al., 2002). This effect is particularly strong at larger scales, and we found significant differences between the $C_{\ell}$ estimate using healpy.anafast (no correction for the mask) and that using NaMaster (partially correct for the mask), as expected.

We constrain the power spectra in the range $\ell<2000$ for two reasons. First, the resolution of the CIB maps is 5’, and any information above $\ell=2000$ would be sensitive to inaccuracies in the beam correction. Second, the underlying H i data that are used to build the dust model have a resolution of $16.2^{\prime}$. Hence, dust residuals on very small scales cannot be removed by the low-resolution H i data and contaminate the power spectrum.

To avoid ringing at the smaller scales in the power spectrum, we apodize the mask as described in Section 2.4. To measure the CIB angular power spectra at a given frequency, we measure the cross power spectra between two CIB maps built on the two ring halves, eliminating the correlated noise term that would impact the smaller scales. The power spectra are furthermore binned with a bin factor of 64 $\ell$ modes per bin. Lastly, we use the covariance of the power spectra as a measure of the uncertainty in our power spectra analysis. We present a more detailed discussion of this in Section 7.1.1. We validated the stage of estimation of the power spectrum extensively using simulations.

5 Map Results

In this section, we present the results exemplary for the 545 GHz channel as an example. The corresponding results for the other frequencies can be found in Appendix A.

We present the input FIR emission and the modeled foreground dust in Figure 7, which nicely shows how the H i-based dust model captures the large-scale cirrus features, while not picking up the CIB fluctuations that can be seen in the intensity map. Moreover, no imprint of the underlying patches in which the dust model is computed is visible.

Figure 8 presents the large-scale, full-resolution 545 GHz map. This is based on the full mission data and has a resolution of $5^{\prime}$. This map covers a total of 18.3$\%$ of the sky.

A smaller sky region that shows the total FIR intensity and the CIB signal after foreground subtraction is presented in Figure 9.

As described in Section 3, we convolved the parameter maps with a $3^{\circ}$ Gaussian kernel to avoid any edge effects in the dust model. An example of such a smoothed parameter map is given in Figure 10, where we show the spatially varying offset $\beta_{\nu}$ (see Eq. (3)). We note that the spatially varying offset is a critical parameter, without which the foreground model fails to provide us with an accurate CIB map (see also Section 7.1.5).

6 Results for the Power Spectra

In the following, we describe the results we have obtained from the angular power spectra of the CIB $C_{\ell}^{\nu_{1}\nu_{2}}$, and from the CIB–CMB lensing cross correlation $C_{\ell}^{T\kappa}$. The formalism to obtain these power spectra has been described in Section 4.

6.1 CIB auto power spectra

The resulting CIB angular auto and cross power spectra $C_{\ell}^{\nu_{1}\nu_{2}}$ for the different frequencies are shown in Figure 11. We include a comparison to the data points and model of Planck Collaboration (2014 XXX), as well as the CIB extracted from the model presented in Mak et al. (2017). We did not correct the results that are based on previous Planck releases for the slightly different calibration, because these are not explicitly studied and are only of the order of $\sim 1-2\%$.

We find an excellent agreement with these previous studies in general, but several differences are worth discussing. First, we observe that Mak et al. (2017) find more power on large scales than both Planck Collaboration (2014 XXX) and the results presented here. This excess signal hints at residual contamination by Galactic dust, and could be the result of the power spectrum-based dust modeling in Mak et al. (2017), as opposed by the map-based dust removal here and in Planck Collaboration (2014 XXX).

Second, toward smaller scales and toward higher frequencies, we find a slightly higher CIB level than what is found in Planck Collaboration (2014 XXX). Despite the steep power law that describes the angular power spectrum of Galactic dust Miville-Deschênes et al. (2002, $\propto\ell^{-2.7}$), we show in Section 7.1.2 that this excess results from dust contamination.

For all CIB cross power spectra, we also estimate the cross-correlation coefficient spectrum $\rho_{\ell}^{\nu_{1}\,\nu_{2}}$, which is given by

[TABLE]

The resulting plot is shown in Figure 12. We furthermore compute the correlation coefficients averaged for the multipole range $150<\ell<1000$ (Table 1). The values are in line with the ones presented in Planck Collaboration (2014 XXX, their Table 11) and with Mak et al. (2017, their Table 5), which is expected given the agreement of the individual power spectra.

6.2 CIB–CMB lensing cross power spectra

The cross correlation of the CIB with the CMB lensing signal is a powerful tool to test the accuracy of the CIB maps, because any missing CIB flux will directly result in a weaker cross correlation. Furthermore, these two quantities are derived completely independently and their comparison serves as a good systematic and quality check. We present the result of this analysis in Figure 13.

We find that our results on the CIB–CMB lensing cross correlation are in good agreement with the results obtained in Planck Collaboration (2014 XVIII), and that the error bars are comparable as well for the higher frequencies. Consequently, our results are also in good agreement with the model presented in Maniyar et al. (2018), which is fit to the data from Planck Collaboration (2014 XVIII). On top of that, this model uses the CIB auto power spectra and measurements of the CIB mean brightness from Planck Collaboration (2014 XXX). These data sets are then used to constrain the cosmic star formation history, the redshift-dependent bias, and the mass of galaxies that contribute most to the CIB emission. We also note that we omitted the two data points at the largest scales from Planck Collaboration (2014 XVIII) because they are only upper limits. Hence, the model of Maniyar et al. is also only mildly constrained on scales $\ell\lesssim 400$.

The dust cleaning allows us to access the large scales at which the cirrus contamination makes any insights difficult. Interestingly, we find that the peak of the power spectrum is at smaller scales than found previously (Planck Collaboration, 2014 XVIII; Maniyar et al., 2018), and that it falls off more quickly toward larger scales. This result is robust and does not depend on the value of the different hyper parameters that we set for the analysis, which we demonstrate in Section 7.

A more detailed analysis of the advantages of this dust removal is shown in Section 7.1.6. For the CIB results of Planck Collaboration (2016 XVII) based on the GNILC component separation, we find a weaker cross correlation with the CMB lensing signal, presumably due to an oversubtraction of the CIB. This is in agreement with the findings in Maniyar et al. (2019, Appendix A).

7 Validation

The key challenge in the separation of CIB emission and Galactic dust emission is the validation: How can we convince ourselves that we removed all the Galactic dust but none of the CIB signal? How do we quantify the dust residuals?

7.1 Internal validation

7.1.1 Estimating uncertainties for the power spectra

The computation of the error bars for the power spectra is not straightforward, as our comparison of different measures shows (Figure 14). Here, we compare different estimators for the uncertainty. First, the analytical approximate estimate of the error on the CIB–CMB lensing cross power spectrum assuming a Gaussian signal and neglecting mode coupling is given by

[TABLE]

The $C_{\ell}$ and the $N_{\ell}$ describe the theoretical signal and noise power spectrum of the CIB temperature $T$ and the lensing convergence $\kappa$. $f_{\rm sky}$ is the effective sky fraction of the maps, and $\Delta\ell$ is the number of $\ell$-modes per bin. Note that since we use the theoretical CIB and lensing signal, we do not capture the uncertainties coming from foreground residuals this way.

We also compute the error bars based on the diagonal of the NaMaster-based covariance matrix of the power spectrum estimate, which we demonstrate to be a reliable measure and hence use throughout this work. Both the analytical and the covariance estimates agree well with the results we obtain from running 100 Gaussian sky simulations. For these, we simulate the CIB and the lensing signal and noise fields, and then recover the power spectra via PyMaster in the presence of the mask. The input power spectra for the CIB, the CMB lensing, and their correlation are taken from Planck Collaboration (2014 XXX), Planck Collaboration (2014 XVIII), and Maniyar et al. (2018), respectively.

On the other hand, simply measuring the error bars by binning the results from an estimator that does not correct for the mode–mode coupling underestimates the error significantly.

7.1.2 $N_{{\mathrm{H\,\mathsc{i}}}}$ thresholds

We investigated how different Galactic plane masks and $N_{{\mathrm{H\,\mathsc{i}}}}$ thresholds affect our results. The goal is to find the optimal trade-off between a large sky fraction and a small contamination by CO-dark molecular gas and Galactic dust (cirrus).

To this end, we test different threshold values and inspect the CIB auto power spectrum (Figure 15). We find that convergence is only reached for very low H i column densities, which does greatly reduce the available sky fraction. For all results presented here, we decided to limit the analysis to regions with $N_{{\mathrm{H\,\mathsc{I}}}}<2.5\times 10^{20}\rm cm^{-2}$, which is a good compromise of residual dust at small scales and a large sky fraction of $18.3$%. This is further confirmed in Section 7.1.6, where we study the uncertainty of the CIB–CMB lensing cross correlation in more detail.

7.1.3 Size of the patches

We also investigated the different choices of the Nside of the superpixels, in which we locally perform the component separation. If these superpixels are chosen to be too large, then the complexity of the H i/FIR relation will not be properly captured. If they are too small, the number of independent sight lines will be too small, and the model will eventually remove CIB structures due to chance correlations with the H i (see also Section 3.1).

We hence repeat the CIB component separation for different choices of the superpixel sizes. Here, we use HEALPix Nsides ranging from 8 to 32. We evaluate the results for the CIB auto power spectrum and for the CIB–CMB lensing cross correlation (Figure 16). The former is difficult to interpret, because there are two effects at play. First, a smaller patch enables us to capture the Galactic dust better, reducing the contamination of the CIB at large scales. Second, the high-pass filtering of the CIB angular power spectrum due to working on individual patches could also potentially remove power on these scales (see Section 3.1).

To overcome this degeneracy, the CIB–CMB lensing cross power spectrum can be conveniently used because residual Galactic dust will have a very different effect. Instead of biasing the power spectrum, the dust will only lead to a higher noise, because it is uncorrelated with the CMB lensing signal.

Based on the analysis of the cross power spectrum of CIB and CMB lensing, we find that the additional noise from residual Galactic dust is very small, even for very large patches (Nside 8) for which the CIB auto power spectrum is strongly contaminated. We furthermore find that the high-pass filtering of the power spectrum due to working on individual patches is a minor effect and only affects the results below $\ell\lesssim 70$, where it biases the result by approximately $2\sigma$. Based on the present analysis, we choose an Nside of 16 for the component separation.

7.1.4 Histograms and Gaussianity

We also inspected the histograms of the recovered CIB at different frequencies to see whether our reconstruction is Gaussian-distributed. In particular, we look for signs of residual CO-dark molecular gas, which would be seen as a heavy tail toward positive residual values.

We present the histogram in Figure 17. To quantify the distributions, we fit a Gaussian to the histograms, both for the full data set and only for the rising flank that is uncontaminated by the dark gas.

As expected, we find the imprint of the dark gas, visible as an excess toward positive residuals. This scales with frequency, where the higher frequencies are more prone to contain Galactic dust instead of CIB signal. Especially for the 857 GHz band, we observe that we are probing the emission close to the peak of the modified blackbody spectrum of Galactic dust.

These histograms are also the basis for the residual masking scheme, which iteratively removes the foreground, evaluates the histogram of the resulting CIB, and masks outliers that exceed a threshold of $3\sigma$.

7.1.5 Offset in the FIR/H i relation

Lastly, we also omitted the offsets in the FIR/H I relation, not allowing a spatially varying zero-point. We find that this leads to widespread failure of the component separation, which is visible both in the map and in the auto power spectrum. This is consistent with the spatially varying offset that we presented in Figure 10, and which shows a large-scale variation of this zero-point. One advantage of this zero-point is that it can capture dust emission unrelated to the H i, such as that from the warm ionized medium (WIM) (Lagache et al., 1999). With Galactic dust dominating the CIB signal on large scales, this additional zero-point is also well suited to address spatial variations in the FIR/H i relation without biasing the CIB signal.

7.1.6 Uncertainty of the CIB–CMB lensing cross correlation

An important question to ask is how much of an impact the selected sky fraction has on the precision of the power spectrum. For the CIB–CMB lensing cross correlation, two effects are at play: a larger sky fraction implies a larger sample and hence smaller errors, but this additional area is also subject to dust residuals, which increase the errors. We present a comparison of the error bars for the CIB–CMB lensing correlation $C_{\ell}^{T\kappa}$ in Figure 18. These error bars are computed by simply binning the full-resolution power spectra and then evaluating the standard deviation in each bin. While we do not use this for any of the final error bars in our analysis, it is useful here because it allows us to compare the relative level of the errors without having to assume the underlying theory power spectra, which is particularly difficult for the residual dust contamination.

We find that for a wide range of sky fractions ($10\%$ – $34\%$), the two effects described above cancel out almost perfectly, and the derived uncertainties are mostly independent of the sky fraction. While this holds for the CIB–CMB lensing cross correlation, the effect is very different for the CIB auto power spectra. Here, the residual dust acts as a bias and not as an uncorrelated quantity that increases the noise (see Figure 15).

To quantify the effect of the foreground dust removal on the CIB–CMB lensing cross correlation, we compare the error bars for this cross correlation for different CIB maps. First, we use the CIB map derived here. Second, we use the raw, uncorrected FIR intensity within the footprint presented here. Lastly, we use the FIR intensity for the entire sky area for which lensing data are available ($66.9\%$). The result (Figure 19) shows that for the large scales ($\ell\lesssim 800$), the foreground removal is crucial to optimize the accuracy of the cross correlation. Even at smaller scales, this cleaning allows our results to be equivalent to or better than the raw intensity for a sky fraction that is 3.7 times larger (18.3$\%$ vs. $66.9\%$).

7.2 External validation

7.2.1 CIB maps from Planck XXX

Aside from the comparison of the power spectra derived here with previous analyses, we also analyzed the differences at the map level. Here, we use the 545 GHz channel and compare our results to the those obtained in Planck Collaboration (2014 XXX, P14). We inspected the smaller individual fields which are based on data from the Green Bank telescope (GBT), as well as the larger so-called GASS field.

To perform this comparison, we convolved the higher-resolution data from P14 to the $5^{\prime}$ resolution of the CIB maps presented here.

Small fields

The results of this comparison are shown in Figure 20. We find that the two CIB estimates show great agreement at the map level, with the differences being dominated by large scales. Several factors contribute to these differences.

First, our work and that conducted in P14 differ fundamentally in the way the dust cleaning is done spatially. For the 2014 results, each field was analyzed separately, and the dust-to-gas ratios were determined for the entire field. Our work uses smaller patches, based on the HEALPix grid, and performs the dust cleaning in each of these separately. Nonetheless, we point out that no sign of this underlying structure can be seen in the CIB difference maps.

Second, the underlying H i data, as well as the selection of the H i data in velocity space, differ between the present work and that in P14. For the small fields that are the basis for the comparison shown here, P14 used data from the GBT that was later published as the GHIGLS data set (Martin et al., 2015). In their analysis, the H i data are individually binned into maps of H i column density that represent the physical phases of the low- and intermediate-velocity clouds (LVC/IVC) (Putman et al., 2012). In our analysis presented here, we use the slightly lower-resolution HI4PI data, and use the GLM to select features in velocity space, without relying on a manual selection (Section 3.4).

Lastly, additional higher-order effects include different calibrations for the Planck data (PR1 vs. PR3), and slightly different resolutions due to different reprojection strategies.

The GASS field

We conduct the same comparison for the larger GASS field from P14 (Figure 21). Here, we bring both CIB maps to the GASS resolution of $16.2^{\prime}$. While our analysis of this field does not differ from the one in the previous section, the P14 approach differs from that applied to the smaller fields.

For the smaller fields, it was assumed that the dust-to-gas ratio is constant within the field, and only varies for the HVC/IVC/LVC phases. Moreover, the resolution of the GBT H i data is $9^{\prime}$, which is close to the Planck resolution at these frequencies ($4^{\prime}$–$5^{\prime}$). The larger GASS field in contrast is too large to assume a constant dust-to-gas ratio, hence a map of dust emissivities and offsets is constructed for patches with a diameter of $15^{\circ}$, centered on the HEALPix pixels with Nside 32. These parameter maps are then smoothed, very similar to the analysis presented here. The main difference from the present study is that P14 only used the map of total local H i column density map for the GASS field, while we use the velocity-resolved data and apply the GLM. This additional degree of freedom now allows us to account for different emissivities along the line of sight, thereby reducing the residual dust contamination.

7.2.2 CIB maps from GNILC

A different approach to extract CIB maps from Planck data has been presented in Planck Collaboration (2016 XLVIII). Instead of using H i-based dust templates, the authors implement the GNILC component separation technique. This is based on the fundamentally different angular power spectra of Galactic dust and the CIB. These maps cover over 60% of sky and are available for the 353, 545, and 857 GHz bands. It should be noted that the primary goal of the GNILC maps is to clean Galactic dust from CIB emission to better determine the properties of Galactic dust such as temperature $T$ and dust spectral index $\beta$.

The most straightforward comparison can be obtained through cross-correlating the CIB against the CMB lensing convergence, which we present in Figure 13. There is no straightforward way to bias this cross correlation toward a higher correlation, but it can easily be lowered by not capturing the full CIB emission. This is the case for the GNILC maps, which were designed to clean the Galactic dust from CIB contamination. The cross correlation with the CMB lensing is systematically lower than what is found in Planck Collaboration (2014 XXX), as already shown in Maniyar et al. (2019), and also in this work.

8 Conclusions

We have presented new large-scale $5^{\prime}$ CIB maps, based on the Planck data and an H i-based foreground removal strategy. These maps are made publicly available for the 353, 545, and 857 GHz data, and for a range of sky fractions.

The foreground removal strategy has been advanced to be more robust, automatic, and to include more information on the three-dimensional H i data. We have validated our results through comparisons with previous studies and find strong agreements both for the images and for the angular power spectra. Additional internal validation shows that the Galactic dust residuals are present down to the smallest scales, depending on the frequency and on the H i threshold that is applied to the maps.

Consequently, different maps should be used for different purposes: studies that are based on the CIB auto power rely on the smallest dust contamination possible, and should therefore use the more aggressive masks. More specifically, we recommend using a thresholds of $2.5\times 10^{20}\;\rm cm^{-2}$ (353 GHz), $2.0\times 10^{20}\;\rm cm^{-2}$ (545 GHz), and $1.8\times 10^{20}\;\rm cm^{-2}$ (857 GHz). For cross correlation analyses, these dust residuals are only a source of noise and do not introduce a bias, hence a larger sky fraction might be preferable to increase the sample size.

For the CIB–CMB lensing cross correlation, we find that the peak is at smaller scales than presented in previous studies, and that the power spectrum falls off slightly faster toward larger scales. This is independent of the various parameters such as sky fraction and patch size for the dust modeling.

Despite the excellent agreement with previous studies, we note that the calibration differences between the three Planck releases complicate the comparison, and a complete cross calibration analysis between the different releases would be required. Part of the discrepancy could come from different absolute brightnesses in the different releases (up to 3.4% at 857 GHz, Planck Collaboration, 2016 VIII).

We anticipate that these maps will be particularly useful for cross-correlation studies with other tracers of the large-scale structure, such as the CMB lensing signal and galaxy surveys.

We make all our CIB data products available online 666https://doi.org/10.7910/DVN/8A1SR3. This includes the following:

The CIB maps at 353, 545, and 857 GHz. We provide these for the Planck full mission data, as well as for the odd-/even-ring data splits. 2. 2.

The H i-based dust model for all frequencies. 3. 3.

Each map comes with a mask, where we proved both the boolean and the apodized masks. 4. 4.

The effective window functions and the FWHM of the maps, which are required to deconvolve the angular power spectra. 5. 5.

Binned power spectra for the different CIB cross correlation spectra $C_{\ell}^{\nu_{1}\nu_{2}}$, and for the CIB–CMB lensing cross correlation $C_{\ell}^{T\kappa}$ 6. 6.

Python notebooks and scripts that illustrate how to work with the data are available online777https://github.com/DanielLenz/PlanckCIB.

8.1 Outlook

A central difficulty in developing component separation for the Galactic dust and the CIB is the validation of the results. Simulations would be a very powerful tool to have in this context, but would require an in-depth understanding of the interplay of the gaseous and dusty ISM for the full sky and with resolutions down to $1^{\prime}$.

The next step in advancing this work is to include multi-scale and multi-frequency information. The former has already been demonstrated in Planck Collaboration (2016 XLVIII), and a combination with the work presented here would be very promising. While the frequency information would not contribute much to disentangle the two components, a combination of the H i-based dust model with existing algorithms for CMB component separation could be a major advancement.

Acknowledgements

Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We acknowledge financial support from “Programme National de Cosmologie et Galaxies” (PNCG) of CNRS/INSU, France. D.L. acknowledges hospitality from the Laboratoire d’Astrophysique de Marseille, where part of this work was completed. Part of this work has been carried out thanks to the support of the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French government program managed by the ANR. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 788212).

We are very grateful to Paolo Serra, Brandon Hensley, Abhishek Maniyar, Bryan Steinbach, Agnès Ferte, Tzu-Ching Chang, François Boulanger, Eric Huff, Emmanuel Schaan, Simone Ferraro, and Matthieu Béthermin for long and very fruitful discussions. We thank Peter Kalberla for providing the Milky Way model of Galactic rotation.

This research has made use of NASA’s Astrophysics Data System, matplotlib (Hunter, 2007), SciPy (Jones et al., 2001–), NumPy (Van Der Walt et al., 2011), scikit-learn (Pedregosa et al., 2011), as well as Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013, 2018). Some of the results in this paper have been derived using the HEALPix (Górski et al., 2005) package

Part of this work was carried out under the program Milky-Way-Gaia of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02.

Appendix A CIB maps

We present the CIB maps at 353 and 857 GHz (Figs. 22 and 23). The 545 GHz is shown in Figure 8. These images are generated for an $N_{{\mathrm{H\,\mathsc{i}}}}$ threshold of $2.5\times 10^{20}\;\rm cm^{-2}$; further maps can be obtained through the released data.

Appendix B Comparison with P14 results

We present additional map-based comparisons of the 545 GHz CIB maps derived here with the ones presented in Planck Collaboration (2014 XXX). Figure 24 shows the SP and the Bootes field, Figure 25 the EBHIS field.