Issue 
A&A
Volume 665, September 2022



Article Number  A52  
Number of page(s)  22  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202243710  
Published online  09 September 2022 
Cosmic star formation history with tomographic cosmic infrared backgroundgalaxy crosscorrelation
^{1}
Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum, Germany
email: yanza21@astro.rub.de
^{2}
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
^{3}
Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02668 Warsaw, Poland
^{4}
Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA
^{5}
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
Received:
4
April
2022
Accepted:
14
June
2022
In this work we present a new method for probing the star formation history of the Universe, namely tomographic crosscorrelation between the cosmic infrared background (CIB) and galaxy samples. The galaxy samples are from the KiloDegree Survey (KiDS), while the CIB maps are made from Planck sky maps at 353, 545, and 857 GHz. We measure the crosscorrelation in harmonic space within 100 < ℓ < 2000 with a significance of 43σ. We model the crosscorrelation with a halo model, which links CIB anisotropies to star formation rates (SFRs) and galaxy abundance. We assume that the SFR has a lognormal dependence on halo mass and that the galaxy abundance follows the halo occupation distribution (HOD) model. The crosscorrelations give a bestfit maximum star formation efficiency of η_{max} = 0.41_{−0.14}^{+0.09} at a halo mass log_{10}(M_{peak}/M_{⊙}) = 12.14 ± 0.36. The derived star formation rate density (SFRD) is well constrained up to z ∼ 1.5. The constraining power at high redshift is mainly limited by the KiDS survey depth. We also show that the constraint is robust to uncertainties in the estimated redshift distributions of the galaxy sample. A combination with external SFRD measurements from previous studies gives log_{10}(M_{peak}/M_{⊙}) = 12.42_{−0.19}^{+0.35}. This tightens the SFRD constraint up to z = 4, yielding a peak SFRD of 0.09_{−0.004}^{+0.003} M_{⊙} yr^{−1} Mpc^{−3} at z = 1.74_{−0.02}^{+0.06}, corresponding to a lookback time of 10.05_{−0.03}^{+0.12} Gyr. Both constraints are consistent, and the derived SFRD agrees with previous studies and simulations. This validates the use of CIB tomography as an independent probe of the star formation history of the Universe. Additionally, we estimate the galaxy bias, b, of KiDS galaxies from the constrained HOD parameters and obtain an increasing bias from b = 1.1_{−0.31}^{+0.17} at z = 0 to b = 1.96_{−0.64}^{+0.18} at z = 1.5, which highlights the potential of this method as a probe of galaxy abundance. Finally, we provide a forecast for future galaxy surveys and conclude that, due to their considerable depth, future surveys will yield a much tighter constraint on the evolution of the SFRD.
Key words: cosmology: observations / diffuse radiation / largescale structure of Universe / galaxies: star formation
© Z. Yan et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Understanding the star formation activity in galaxies is central to our understanding of the evolution of galaxies in the Universe (Tinsley 1980). Moreover, the observed relationship between star formation and other physical processes implies that there exists complex interactions within galaxies between gases, stars, and central black holes (for example, through feedback from supernovae or supermassive black holes). Star formation activity can be described by the star formation rate density (SFRD), defined as the stellar mass generated per year per unit volume. By studying the SFRD of galaxies at different redshifts, we can understand the cosmic star formation history. In the local Universe, the star formation rate (SFR) can be explored via imaging the molecular gas in nearby galaxies (Padoan et al. 2014; Querejeta et al. 2021). For distant galaxies, the SFRD is typically studied via multiwavelength observations (Gruppioni et al. 2013; Magnelli et al. 2013; Davies et al. 2016; Marchetti et al. 2016), which focus on the emission properties of galaxy populations over a wide range of wavelengths. By assuming the luminosity functions and the luminositystellar mass relations at different wavelengths, from ultraviolet (UV) to farinfrared, these works derive the SFR of galaxy populations from multiwavelength luminosities. Previous multiwavelength studies have reached a consensus (see Madau & Dickinson 2014 for a review). From these studies, we are confident that star formation in the Universe started between 6 ≲ z ≲ 20 and reached its peak at z ∼ 2 (corresponding to a lookback time of ∼10.3 Gyr), at a rate approximately ten times higher than what is observed today. Due to a subsequent deficiency of available gas fuel, star formation activity has been steadily decreasing since z ∼ 2.
Dust in stellar nurseries and in the interstellar medium of galaxies absorbs UV radiation produced by shortlived massive stars and reradiates this energy as infrared (IR) emission (Kennicutt 1998; Smith et al. 2007; Rieke et al. 2009). Therefore, extragalactic IR radiation can be used to study star formation throughout cosmic history. Moreover, one can also use the spectral energy distribution (SED) of IR galaxies to study the thermodynamics of interstellar dust in these galaxies (Béthermin et al. 2015, 2017). However, dusty starforming galaxies beyond the local Universe are typically highly blended, given the sensitivity and angular resolution of modern IR observatories, because they are both faint and numerous. This makes it more difficult to study them individually in terms of their SFRD at higher redshift (see, for example, Nguyen et al. 2010)^{1}.
Conversely, the projected allsky IR emission, known as the cosmic infrared background (CIB; Partridge & Peebles 1967), encodes the cumulative emission of all dusty starforming galaxies below z ∼ 6. It is therefore a valuable tool that can be leveraged to investigate the SFRD. However, measurements of the CIB itself are complicated: imperfect removal of point sources and foreground Galactic emission can lead to bias in the measured CIB signals. Nevertheless, measurement of the CIB is more robust to the various selection effects and sample variance uncertainties that affect galaxybased studies, and which depend highly on instrumental setup and observation strategies. The CIB was first detected by the Cosmic Background Explorer satellite (Dwek et al. 1998) and was then accurately analysed by Spitzer (Dole et al. 2006) and Herschel (Berta et al. 2010), via large IR galaxy samples, and by Planck (Planck Collaboration XXX 2014), via mapping of the CIB from its sky maps. As with the cosmic microwave background (CMB), the key to mapping the CIB is accurately removing the foreground Galactic emission. However, unlike the CMB, the peak frequency of the CIB is around ∼3000 GHz, which is dominated by the thermal emission from Galactic dust. Moreover, unlike the CMB and the thermal SunyaevZel’dovich (tSZ) effect (Sunyaev & Zeldovich 1972), the CIB has no unique frequency dependence. Therefore, the CIB is more difficult to extract from raw sky maps than the CMB or the tSZ effect and is generally restricted to high Galactic latitudes. There are multiple published CIB maps that are generated by various methods (differing primarily in the removal of the foreground). Planck Collaboration XXX (2014) and Lenz et al. (2019) remove Galactic emission by introducing a Galactic template from HI measurements, while Planck Collaboration XLVIII (2016) disentangles the CIB from Galactic dust by the generalised needlet internal composition (GNILC; Remazeilles et al. 2011) method. The mean CIB emission measured from maps gives the mean energy emitted from star formation activities, while the anisotropies of the CIB trace the spatial distribution of starforming galaxies and the underlying distribution of their host dark matter halos (given some galaxy bias). Therefore, analysing the CIB anisotropies via angular power spectra has been proposed as a new method to probe cosmic star formation activity (Dole et al. 2006).
Angular power spectra have been widely used in cosmological studies. By crosscorrelating different tracers, one can study cosmological parameters (Kuntz 2015; Singh et al. 2017), the cosmic thermal history (Pandey et al. 2019; Koukoufilippas et al. 2020; Chiang et al. 2020; Yan et al. 2021), baryonic effects (Hojjati et al. 2017; Tröster et al. 2021b), the integrated SachsWolfe (ISW) effect (Maniyar et al. 2019; Hang et al. 2021), and more. Cosmic infrared background maps have been extensively used to study dusty starforming galaxies via autocorrelations (Shang et al. 2012; Planck Collaboration XXX 2014; Maniyar et al. 2018) and crosscorrelations with other largescale structure tracers (Cao et al. 2020; Maniyar et al. 2021). Clusteringbased CIB crosscorrelation has been used to study star formation in different types of galaxies; for example, Serra et al. (2014) analyse luminous red galaxies (LRGs), Wang et al. (2015) analyse quasars, and Chen et al. (2016) analyse submillimetre galaxies (SMGs). The tracers used in these studies are either projected sky maps or galaxy samples with wide redshift ranges, leading to model parameters that describe the redshift dependence being highly degenerate. Schmidt et al. (2014) and Hall et al. (2018) crosscorrelate the CIB with quasars at different redshifts, yielding an extensive measurement of the evolution of the CIB signal in quasars. However, these studies are restricted to active galaxies and therefore may miss contributions from the wider population of galaxies.
This paper proposes a new clusteringbased measurement that allows us to study the cosmic star formation history with the CIB: tomographic crosscorrelation between the CIB and galaxy number density fluctuations. That is, crosscorrelating the CIB with galaxy samples in different redshift ranges (socalled tomographic bins) to measure the evolution of the CIB over cosmic time. Compared with other largescale structure tracers, galaxy number density fluctuations can be measured more directly. Firstly, galaxy redshifts can be determined directly via spectroscopy, although this process is expensive and must be restricted to particular samples of galaxies and/or small onsky areas. Alternatively, widearea photometric surveys provide galaxy samples that are larger and deeper than what can be observed with spectroscopy, and whose population redshift distribution can be calibrated to high accuracy with various algorithms (see Salvato et al. 2019 for a review). Successful models have been proposed to describe galaxy number density fluctuations. On large scales, the galaxy distribution is proportional to the underlying mass fluctuation; on small scales, its nonlinear behaviour can be modelled by a halo occupation distribution (HOD; Zheng et al. 2005) model. With all these practical and theoretical feasibilities, galaxy density fluctuations have long been used to study various topics in largescale structure, including reionisation (Lidz et al. 2008), cosmological parameters (Kuntz 2015), and the ISW effect (Hang et al. 2021). In the near future, the CanadaFrance Imaging Survey (CFIS; Ibata et al. 2017), the Rubin Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration 2009), and the Euclid (Laureijs et al. 2010) mission will reach unprecedented sky coverage and depth, making galaxy number density fluctuation a ‘treasure chest’ from which we will learn a lot about our Universe.
The CIB is generated from galaxies and so should correlate with galaxy distribution. Limited by the depth of current galaxy samples, CIBgalaxy crosscorrelations are only sensitive to the CIB at low redshift, but this will improve with future galaxy surveys. In this study, we crosscorrelate the galaxy catalogues provided by the Kilodegree Survey (KiDS; de Jong et al. 2013) with CIB maps constructed at 353, 545, and 857 GHz to study the SFRD. The galaxy samples are divided into five tomographic bins extending to z ∼ 1.5. Although the measurements are straightforward, modelling the CIB is more challenging than many other tracers. Firstly, SFRs and dust properties are different from galaxy to galaxy, and we do not have a clear understanding of both in a unified way. Previous studies take different models for the CIB: Planck Collaboration XXX (2014) and Shang et al. (2012) use a halo model by assuming a lognormal luminositytohalo mass (L − M) relation for the IR and a greybody spectrum for extragalactic dust; Maniyar et al. (2018) and Cao et al. (2020) use the linear perturbation model with empirical radial kernel for the CIB; and Maniyar et al. (2021) propose an HOD halo model for the CIB. In this work we use the Maniyar et al. (2021, M21 hereafter) model since it explicitly links the redshift dependence of the CIB with the SFR.
This paper is structured as follows. In Sect. 2 we describe the theoretical model we use for the crosscorrelations. Section 3 introduces the dataset that we are using. Section 4 presents the method for measuring crosscorrelations, as well as our estimation of the covariance matrix, likelihood, and systematics. Section 5 presents the results. Section 6 discusses the results and summarises our conclusions. Throughout this study, we assume a flat Λ cold dark matter cosmology with the fixed cosmological parameters from Planck Collaboration VI (2020) as our background cosmology: (h, Ω_{c}h^{2}, Ω_{b}h^{2}, σ_{8}, n_{s})=(0.676, 0.119, 0.022, 0.81, 0.967).
2. Models
2.1. Angular crosscorrelation
Both the galaxy and the CIB angular distributions are formalised as the lineofsight projection of their 3D distributions. This subsection introduces the general theoretical framework of the angular crosscorrelation between two projected cosmological fields. For an arbitrary cosmological field u, the projection of its 3D fluctuations (i.e. anisotropies) is written as
where is the 2D projection in the angular direction , and is the fluctuation of u in 3D space at the coordinate , where χ is the comoving distance. The kernel W^{u}(χ) describes the lineofsight distribution of the field^{2}.
We measure the angular crosscorrelation in harmonic space. In general, the angular crosscorrelation between two projected fields, u and v, at the scales ℓ ≳ 10 are well estimated by the Limber approximation (Limber 1953; Kaiser 1992):
where P_{uv}(k, z) is the 3D crosspower spectrum of associated 3D fluctuating fields u and v:
Generally, we can model a largescale cosmological field as a biasedtracer of the underlying mass, mainly in the form of dark matter halos (Cooray & Sheth 2002; Seljak 2000). In such a halo model, P_{uv}(k) is divided into the twohalo term, which accounts for the correlation between different halos, and the onehalo term, which accounts for correlations within the same halo. Smith et al. (2003) points out that simply adding the one and twohalo up yields a total power spectrum that is lower than that from cosmological simulations in the transition regime (k ∼ 0.5 h Mpc^{−1}). Mead et al. (2021) estimates that this difference can be up to a level of 40%, so one needs to introduce a smoothing factor α to take this into account. The total power spectrum is then given by
The redshift dependence of α is given by the fitting formula in Mead et al. (2021)^{3}. In Fig. 1, we plot the one and twohalo terms (dashdotted purple and red solid lines, respectively) of the CIBgalaxy crosscorrelation power spectrum (to be introduced below), as well as their sum (the dashed black line) and the smoothed power spectrum (the solid black line). It is clear that the smoothing changes the power spectrum in the transition regime.
Fig. 1. Halo model of the power spectrum of CIBgalaxy crosscorrelation at z = 0. The power spectrum in this plot only shows the spatial dependence of the correlation between the CIB and galaxy distribution, with all the irrelevant terms (redshift and frequency dependence) factored out, so the unit is arbitrary. The dashdotted purple line and the solid red line are one and twohalo terms, respectively; the dashed black line is the summation of one and twohalo terms, and the solid black line is the smoothed power spectrum defined in Eq. (4). 
Both one and twohalo terms are related to the profiles of u and v in Fourier space:
where the angled brackets ⟨ ⋅ ⟩ describe the ensemble average of the quantity inside. At a given redshift, P^{lin}(k) is the linear power spectrum, dn/dM is the halo mass function (number density of dark matter halos in each mass bin), b_{h} is the halo bias, and p_{u}(kM) is the profile of the tracer u with mass M in Fourier space:
where p_{u}(rM) is the radial profile of u in real space. In this work, we employ the halo mass function and halo bias given by Tinker et al. (2008, 2010), respectively, in accordance with M21.
2.2. Galaxy number density fluctuations
The 2D projected galaxy number density fluctuation is measured as
where is the surface number density of galaxies in the direction on sky, and is the average surface number density. Given the redshift distribution of a galaxy sample Φ_{g}(z) (determined by the true lineofsight galaxy distribution and any survey selection functions)^{4}, the projected galaxy density fluctuation is given by
where is the 3D galaxy density fluctuation. The radial kernel for galaxy number density fluctuation is then
The galaxy density fluctuation in a halo with mass M can be described by its number density profile p_{g}(rM), as
where δ^{3D} is the 3D Dirac delta function, N_{c}(M) and N_{s}(M) are the number of central galaxy and satellite galaxies as a function of the halo mass (M), respectively, and p_{s}(rM) is the number density profile of the satellite galaxies. Its Fourier transform will be given in Eq. (14). is the mean galaxy number density at redshift z, which is given by
Though we cannot say anything about galaxy counts for individual halos, their ensemble averages can be estimated via the HOD model (Zheng et al. 2005; Peacock & Smith 2000):
where M_{min} is the mass scale at which half of all halos host a galaxy, σ_{lnM} denotes the transition smoothing scale, M_{1} is a typical halo mass that consists of one satellite galaxy, M_{0} is the threshold halo mass required to form satellite galaxies, and α_{s} is the power law slope of the satellite galaxy occupation distribution. Θ is the Heaviside step function.
In Fourier space, the galaxy number profile is given by
where the dimensionless profile of satellite galaxies p_{s}(kM) is generally taken as the NavarroFrenkWhite (NFW) profile (van den Bosch et al. 2013; Navarro et al. 1996):
where q ≡ kr_{200}(M)/c(M), c is the concentration factor, and the functions {Ci, Si} are the standard cosine and sine integrals, respectively^{5}. The concentrationmass relation in this work is given by Duffy et al. (2008). Here r_{200} is the radius that encloses a region where the average density exceeds 200 times the critical density of the Universe. We take the total mass within r_{200} as the proxy of halo mass because in general r_{200} is close to the virial radius of a halo (Opher et al. 1998)^{6}.
The HOD parameters in Eq. (12) depend on redshift (Coupon et al. 2012). In this work, we fix σ_{lnM} = 0.4 and α_{s} = 1, consistent with simulations (Zheng et al. 2005) and previous observational constraints (Coupon et al. 2012; Ishikawa et al. 2020), and adopt a simple relation for {M_{0}, M_{1}, M_{min}} with respect to redshift. For example, we model M_{0} as in Behroozi et al. (2013):
where a is the scale factor, log_{10}M_{0, 0} is the value at z = 0, while log_{10}M_{0, p} gives the ‘rate’ of evolution^{7}. Therefore, in total we constrain six HOD parameters: {M_{0, 0}, M_{1, 0}, M_{min, 0}, M_{0, p}, M_{1, p}, M_{min, p}}. In practice, we find that the resolution of the CIB map is sufficiently low that this simple formalism fits the data well (Sect. 5).
2.3. Halo model for CIBgalaxy crosscorrelation
The intensity of the CIB (in Jy sr^{−1}) is the lineofsight integral of the comoving emissivity, j_{ν}:
Comparing with Eq. (1), one can define the radial kernel for the CIB to be
which is independent of frequency. Thus, the emissivity j_{ν} is the ‘δ_{u}’ for the CIB, which is related to the underlying galaxy population as
where L_{ν}(z) is the IR luminosity and dn/dL is the IR luminosity function. The second equation assumes that galaxy luminosity is also a function of the mass of the host dark matter halo. Furthermore, like galaxies, the model of the IR luminosity can also be divided into contributions from central and satellite galaxies (Shang et al. 2012; Planck Collaboration XXX 2014). We introduce the IR luminous intensity (i.e. the power emitted per steradian):
where the subscripts ‘c/s’ denote the central and satellite components, respectively. The profile of the CIB in Fourier space is formulated as
Comparing with Eq. (10), one recognises that the quantity f_{ν, c/s}(M) is directly analogous to N_{c/s}(M), and f_{ν}(kM) is the profile term p_{u}(kM) in Eq. (5) for CIB anisotropies. Following the standard practice of van den Bosch et al. (2013), we give the crosscorrelation between the Fourier profile of galaxies and the CIB that is needed for calculating the onehalo term:
We discuss how to model f_{ν, c/s} in Sect. 2.4.
2.4. CIB emissivity and star formation rate
Considering the origin of the CIB, j_{ν} should depend on the dust properties (temperature, absorption, etc.) of starforming galaxies, and on their SFR. In addition, the CIB also directly traces the spatial distribution of galaxies and their host halos. We take the halo model for the CIB from M21. The observed mean CIB emissivity at redshift z is given by
where ρ_{SFR}(z) is the SFRD, defined as the stellar mass formed per year per unit volume (in M_{⊙} yr^{−1} Mpc^{−3}), and S_{eff}(ν, z) is the effective SED of IR emission from galaxies at the given the restframe frequency ν and redshift z. The latter term is defined as the mean flux density received from a source with L_{IR} = 1 L_{⊙}, so it has a unit of Jy/L_{⊙}. We note that we change to the rest frame frequency by multiplying the observed frequency ν by (1 + z). The Kennicutt (1998) constant K is defined as the ratio between SFRD and IR luminosity density. In the wavelength range of 8 − 1000 μm, it has a value of K = 1 × 10^{−10} M_{⊙} yr^{−1} assuming a Chabrier initial mass function (IMF; Chabrier 2003). The derivation of this formula can be found in Appendix B of Planck Collaboration XXX (2014).
The SFRD is given by
where SFR_{tot}(M, z) denotes the total SFR of the galaxies in a halo with mass M at redshift z.
As is shown in Eq. (20), f_{ν}(M, z) can also be divided into components from the central galaxy and satellite galaxies living in subhalos. Following Shang et al. (2012), Maniyar et al. (2018, 2021), we assume that the central galaxy and satellite galaxies share the same effective SED. In the literature, the SED in Eq. (22) are given with different methods: Shang et al. (2012) parametrises the effective SED with a greybody spectrum, while Maniyar et al. (2018, 2021) use fixed effective SED from previous studies. In this work, we follow M21 and take the SED calculated with the method given by Béthermin et al. (2013); that is, we assume the mean SED of each type of galaxies (main sequence, starburst), and weigh their contribution to the whole population in construction of the effective SED. The SED templates and weights are given by Béthermin et al. (2015, 2017). Therefore, central and satellite components differ only in SFR, and the total SFR in Eq. (23) is given by SFR_{tot} = SFR_{c} + SFR_{s}. Combining Eqs. (18), (19), (22), and (23), we recognise that
The final piece of the puzzle for our model is in defining SFR_{c/s}. Following Béthermin et al. (2013), the SFR is given by the baryon accretion rate (BAR, measured in solar masses per year: M_{⊙} yr^{−1}) multiplied by the star formation efficiency η. That is,
For a given halo mass M at redshift z, the BAR is the mean mass growth rate (MGR; also measured in M_{⊙} yr^{−1}) of the halo multiplied by the baryonmatter ratio:
The MGR is given by the fits of (Fakhouri et al. 2010)
where Ω_{m}, Ω_{b}, and Ω_{Λ} are the density parameters of total mass, baryons, and dark energy of the Universe, respectively.
The star formation efficiency is parameterised as a lognormal function of the halo mass, M:
where M_{peak} represents the mass of a halo with the highest star formation efficiency η_{max}. An analysis of average SFRs and histories in galaxies from z = 0 to z = 8 shows that M_{peak} ought to be constant over cosmic time, at a value of M_{peak} ∼ 10^{12} M_{⊙}Behroozi et al. (2013). Therefore, in our model, we assume it to be a constant. And σ_{M}(z) is the variance of the lognormal, which represents the range of masses over which the star formation is efficient. Also, following M21, this parameter depends both on redshift and halo mass:
where z_{c} is the redshift above which the mass window for star formation starts to evolve, with a ‘rate’ described by a free parameter τ. In this work, we fix z_{c} = 1.5, as in M21, because our sample of KiDS galaxies is unable to probe beyond this redshift (see Sect. 3 and Fig. 2).
Fig. 2. Redshift distributions of the KiDS1000 gold galaxy sample (solid and dashed lines) and CIB emissions (dotted lines). Solid lines are the redshift distribution of the KiDS galaxies calibrated by the SOM without specifying lensing weight and are the redshift distributions we use in this work. For comparison, we show in dashed lines the redshift distributions from the SOM calibration with lensing weight, which are used in the standard KiDS1000 cosmological analyses. Coloured bands show the z_{B} ranges of corresponding tomographic bins. Dotted lines are dI_{ν}/dz at 353, 545, and 857 GHz calculated from Eq. (22) with the bestfit parameters in this work. The values follow the right y axis. 
For the central galaxy, the SFR is calculated with Eq. (25), where M describes the mass of the host halo, multiplied by the mean number of central galaxies ⟨N_{c}⟩ as given by Eq. (12):
For satellite galaxies, the SFR depend on the masses of subhalos in which they are located (Béthermin et al. 2013):
where m is the subhalo mass, and SFR is the general SFR defined by Eq. (25). The mean SFR for subhalos in a host halo with mass M is then
where dN_{sub}/dlnm is the subhalo mass function. We take the formulation given by Tinker & Wetzel (2010). Once we get the SFR for both the central and the subhalos, we can add them together and calculate the luminous intensity f_{ν} of a halo with Eq. (22), and then calculate the angular power spectra with the halo model and Limber approximation as discussed in Sects. 2.1 and 2.2.
There are a couple of simplifying assumptions in our model. First of all, we assume that the IR radiation from a galaxy is entirely the thermal radiation from dust, which is generated by star formation activity. However, part of the IR radiation may be generated from nondust radiation, including CO(3−2), freefree scattering, or synchrotron radiation (Galametz et al. 2014). We also assume that central and satellite galaxies have the same dust SED, which might be entirely accurate. In addition, we neglect the difference in quenching in central and satellite galaxies (Wang et al. 2015). However, the IR radiation is dominated by central galaxies, so the differences between central and satellite galaxies will not significantly affect our conclusion. In any case, these limitations need further investigation by future studies.
We note, though, that the measured power spectrum will also contain shotnoise resulting from the autocorrelated Poisson sampling noise. Therefore, the model for the total CIBgalaxy crosscorrelation is
where is the crosscorrelation predicted by the halo model, and S^{νg} is the scaleindependent shot noise. Shotnoise is generally not negligible in galaxy crosscorrelations, especially on small scales. There are analytical models to predict shot noise (Planck Collaboration XXX 2014; Wolz et al. 2018) but they depends on various assumptions, including the CIB flux cut, galaxy colours, galaxy physical evolution (Béthermin et al. 2012), and so on. Each of these assumptions carries with it additional uncertainties. Therefore, in this work, instead of modelling the shot noise for different pairs of crosscorrelations, we simply opt to keep their amplitudes as free parameters in our model^{8}. In practice, we set log_{10}S^{νg} to be free, where S^{νg} is in the unit of MJy sr^{−1}.
With the SFR and SED models introduced above, the redshift distribution of the CIB intensity can be calculated with Eq. (16). The redshift distributions of the CIB intensity at 353, 545, and 857 GHz are shown as dotted lines in Fig. 2. It is clear that CIB emission increases with frequency (in the frequency range we explore here). The peak redshift of CIB emission is z ∼ 1.5, which is close to the redshift limit of our galaxy sample.
In addition, once we have fixed the model parameters in SFR with our measurements, we can calculate ρ_{SFR} by adding up SFR of central and satellite galaxies and employing Eq. (23):
The primary goal of this work is to constrain this ρ_{SFR}(z) parameter with CIBgalaxy crosscorrelation.
3. Data
3.1. KiDS data
We used the lensing catalogue provided by the fourth data release (DR4) of KiDS (Kuijken et al. 2019) as our galaxy sample. KiDS is a widefield imaging survey that measures the positions and shapes of galaxies using the VLT Survey Telescope (VST) at the European Southern Observatory (ESO). Both the telescope and the survey were primarily designed for weaklensing applications. The footprint of the KiDS DR4 (and its corresponding galaxy sample, called the ‘KiDS1000’ sample) is divided into a northern and southern patch, with total coverage of 1006 deg^{2} of the sky (corresponding to a fractional area of the full sky of f_{sky} = 2.2%). The footprint is shown as the transparent patches in Fig. 3. Highquality optical images are produced with VSTOmegaCAM, and these data are then combined with imaging from the VISTA Kilodegree INfrared Galaxy (VIKING) survey (Edge et al. 2013), allowing all sources in KiDS1000 to be photometrically measured in nine optical and nearIR bands: ugriZYJHK_{s} (Wright et al. 2019). The KiDS1000 sample selects galaxies with photometric redshift estimates 0.1 < z_{B} ≤ 1.2. Although the sky coverage of KiDS is relatively small comparing to some galaxy surveys (such as the Dark Energy Survey; Abbott et al. 2016), galaxy photometric redshift estimation and redshift distribution calibration (especially at high redshift) is more reliable in KiDS thanks to the complimentary nearIR information from VIKING (which was codesigned with KiDS to reach complimentary depths in the nearIR bands). Each galaxy in the KiDS1000 sample has ellipticities measured with the Lensfit algorithm (Miller et al. 2013), which allows exceptional control for systematic effects such as stellar contamination (Giblin et al. 2021). The KiDS1000 sample is then further downselected during the production of highaccuracy calibrated redshift distributions (Wright et al. 2020; Hildebrandt et al. 2021) to produce the KiDS1000 ‘gold’ sample^{9}. We used the gold sample for this work as the redshift distributions are most accurately calibrated for these galaxies. We present the information of the galaxy sample that we use in Table 1. This resulting galaxy sample covers redshifts z ≲ 1.5, and is therefore a suitable dataset to trace the history of different components of the largescale structure into the intermediateredshift Universe.
Fig. 3. CIB map at 545 GHz for the Galactic north pole (left) and south pole (right). The transparent regions are the KiDS footprint. 
Information on the KiDS galaxy sample in each tomographic bin.
Although the KiDS survey provides highaccuracy shape measurements for galaxies, we do not use them in this analysis. As is argued in Yan et al. (2021), galaxy number density fluctuations are relatively easy to measure (compared to galaxy shapes) and are immune to the systematic effects inherent to the shape measurement process (including shape measurement accuracy, response to shear, shape calibration error, intrinsic alignment, etc). Moreover, the CIB is generated directly from galaxies, so we expect strong CIBgalaxy correlation signals, which reveal the star formation activity in these galaxies. Therefore, we focus on CIBgalaxy crosscorrelation in this work, allowing us to ignore shape information. However, we note that Tröster et al. (2021b) has made a significant detection of shearCIB crosscorrelation with the 545 GHz Planck CIB map, which can help us understand the connection between halos and IR emissions. We leave such an investigation for future work.
We perform a tomographic crosscorrelation measurement by dividing the galaxy catalogue into five bins, split according to the bestfit photometric redshift estimate z_{B} of each galaxy. These are the same tomographic bins used in the KiDS1000 cosmology analyses (Asgari et al. 2021; Heymans et al. 2021; Tröster et al. 2021a). The redshift distribution of each bin is calibrated using a selforganising map (SOM) as described in Wright et al. (2020), Hildebrandt et al. (2021). A SOM is a 2D representation of an ndimensional manifold, computed using unsupervised machine learning. For redshift calibration, the SOM classifies the distribution of galaxies in multidimensional colourmagnitude space into discrete cells. As galaxy colour correlates with redshift, cells on the SOM similarly correlate with redshift. Using tracer samples of galaxies with known spectroscopic redshift, the associations between galaxies and cells in a SOM can therefore be used to reconstruct redshift distributions for large samples of photometrically defined galaxies. We note that the SOMcalibrated redshift distributions in this study are not the same as Hildebrandt et al. (2021), in which the redshift distributions are calibrated with a galaxy sample weighted by the Lensfit weight. In this work the redshift distributions are calibrated with the raw, unweighted sample. The redshift distributions of the five tomographic bins are shown in Fig. 2. We also plot the SOMcalibrated Φ_{g}(z) with lensing weight as dashed lines. The absolute difference between the means of the two Φ_{g}(z) are at a level of ∼0.01, comparable to the mean Φ_{g}(z) bias given by Hildebrandt et al. (2021), and the difference is more evident in the higherredshift bins. We also show the mean CIB emissions (dotted lines) at 353, 545, and 857 GHz calculated from Eq. (22) with the bestfit parameters in this work.
We utilised maps of relative galaxy overdensity to encode the projected galaxy density fluctuations. These galaxy overdensity maps are produced for each tomographic bin in the HEALPIX (Gorski et al. 2005) format with NSIDE = 1024, corresponding to a pixel size of 3.4 arcmin. For the tth tomographic bin, the galaxy overdensity in the ith pixel is calculated as
where i denotes the pixel index, n_{t, i} is the surface number density of galaxies in the ith pixel and is the mean galaxy surface number density in the tth redshift bin of the KiDS footprint. We note that Eq. (35) is slightly different from Eq. (10) in that it is the mean galaxy overdensity in each pixel while Eq. (10) defines the galaxy overdensity on each point in the sky. In other words, Δ_{g, i} in Eq. (35) is the discretised with the window function corresponding to the HEALPIX pixel. The mask of the galaxy maps for the crosscorrelation measurement is the KiDS footprint, which is presented as the transparent regions in Fig. 3.
3.2. CIB data
In this work, we use the largescale CIB maps generated by Lenz et al. (2019) from three Planck High Frequency Instrument (HFI) sky maps at 353, 545, and 857 GHz (the L19 CIB map hereafter)^{10}. The IR signal generated from galactic dust emission is removed based on an HI column density template (Planck Collaboration XVIII 2011; Planck Collaboration XXX 2014). We use the CIB maps corresponding to an HI column density threshold of 2.0 × 10^{20} cm^{−2}. The CIB mask is a combination of the Planck 20% Galactic mask, the Planck extragalactic point source mask, the molecular gas mask, and the mask defined by the HI threshold. The CIB maps have an overall sky coverage fraction of 25%, and overlap with most of the KiDS south field and part of the north field. The overlapping field covers about 1% of the sky. The CIB signal in the maps is in the unit of MJy sr^{−1} with an angular resolution of 5 arcmin, as determined by the full width at half maximum (FWHM) of the beam. The original maps are in the HEALPIX format with NSIDE = 2048 and we degrade them into NSIDE = 1024 since we do not probe scales smaller than ℓ ∼ 1500.
The Planck Collaboration also makes allsky CIB maps (Planck Collaboration XLVIII 2016) in the three highest HFI frequencies. To make the allsky CIB map, the Planck Collaboration disentangle the CIB signal from the galactic dust emission with the GNILC method (Remazeilles et al. 2011). These maps have a large sky coverage (about 60%) and have been extensively used to constrain the CIB power spectra (Mak et al. 2017; Reischke et al. 2020) and to estimate systematics for other tracers (Yan et al. 2019; Chluba et al. 2017). However, Maniyar et al. (2018), Lenz et al. (2019) point out that when disentangling galactic dust from CIB, there is some leakage of the CIB signal into the galactic dust map, causing biases of up to ∼20% in the CIB map construction. Therefore, we opt to not use the Planck GNILC CIB map in this work at the expenses of sky coverage.
3.3. External SFRD data
In addition to the CIBgalaxy crosscorrelation power spectra we also introduce external SFRD measurements, estimated over a range of redshifts, as additional constraints to our model. The external SFRD measurements are obtained by converting the measured IR luminosity functions to ρ_{SFR} with proper assumptions of the IMF. We refer the interested reader to a review on converting light to stellar mass (Madau & Dickinson 2014). In this work, we use the ρ_{SFR} from (Gruppioni et al. 2013; Magnelli et al. 2013; Marchetti et al. 2016; Davies et al. 2016). We follow Maniyar et al. (2018) to account for different background cosmologies utilised by these studies: we first convert the external ρ_{SFR} values into cosmologyindependent observed intensity between 8 and 1000 μm per redshift bin, according to corresponding cosmologies, and then convert back to ρ_{SFR} with the cosmology assumed in this study.
4. Measurements and analysis
4.1. Crosscorrelation measurements
The crosscorrelation between two sky maps, which are smoothed with the beam window function b_{beam}(ℓ), is related to the real C_{ℓ} as
where denotes the smoothed C_{ℓ} between the sky maps u and v, and b_{pix}(ℓ) is the pixelisation window function provided by the HEALPIX package. In our analysis, we take the Gaussian beam window function that is given by
where . For the KiDS galaxy overdensity map, the angular resolution is much better than the pixel size, so we assume a FWHM = 0 and .
Both the galaxy and the CIB maps are partialsky. The sky masks mix up modes corresponding to different ℓ. The modemixing is characterised by the modemixing matrix, which depends only on the sky mask. We use the NAMASTER (Alonso et al. 2019)^{11} package to correct modemixing and estimate the angular crosspower spectra. NAMASTER first naively calculates the crosscorrelation between two masked maps. This estimation gives the biased power spectrum, which is called the ‘pseudoC_{ℓ}’. Then it calculates the modemixing matrices with the input masks and uses this to correct the pseudoC_{ℓ}. The beam effect is also corrected in this process. The measured angular power spectra are binned into ten logarithmic bins from ℓ = 100 to ℓ = 2000. The high limit of ℓ corresponds to the Planck beam, which has FWHM of 5 arcmin. The low limit is set considering the small skycoverage of KiDS.
The measurements give a data vector of crosscorrelation between five tomographic bins of KiDS and 3 CIB bands, resulting in 15 crosscorrelations bin1, bin2, bin3, bin4, bin5}; ν ∈ {353 GHz, 545 GHz, 857 GHz}. Given the covariance matrix to be introduced in Sect. 4.2, we calculate the squareroot of the χ^{2} values of the null line () and reject the null hypothesis at a significance of 43σ. With these measurements, we are trying to constrain the following model parameters with uniform priors: (i) SFR parameters: {η_{max}, log_{10}M_{peak}, σ_{M, 0}, τ}. (ii) HOD parameters: {log_{10}M_{0, 0}, log_{10}M_{0, p}, log_{10}M_{1, 0}, log_{10}M_{1, p}, log_{10}M_{min, 0}, log_{10}M_{min, p}}. (iii) Amplitudes of the shot noise power spectra: {log_{10}S^{νg}}. Here S^{νg} is in the unit MJy sr^{−1}. The prior boundaries are [ − 12, 8] for all the 15 shot noise amplitudes.
In total, there are 25 free parameters to constrain (see Table 2 for a summary). The number of data points is 3 (frequencies) × 5 (tomographic bins) × 10 (ℓ bins) = 150, so the degreeoffreedom is 125.
Summary of the free parameters, their prior ranges, and the equations that define them.
4.2. Covariance matrix
To estimate the uncertainty of the crosscorrelation measurement, we followed the general construction of the analytical covariance matrix in the literature. Compared with simulation or resamplingbased methods, an analytical method is free of sampling noise and allows us to separate different contributions.
Following (Tröster et al. 2021b), we decompose the crosscorrelation covariance matrix into three parts:
Here Cov is the abbreviation of . We note that both ℓ_{1} and ℓ_{2} are corresponding to ℓ bands rather than a specific ℓ mode. The first term Cov^{G} is the dominant ‘disconnected’ covariance matrix corresponding to Gaussian fields, including physical Gaussian fluctuations and Gaussian noise:
This is the covariance matrix for an allsky measurement. Sky masks introduce nonzero coupling between different ℓ as well as enlarge the variance. To account for this, we used the method given by Efstathiou (2004) and GarcíaGarcía et al. (2019) that is implemented in the NAMASTER package (Alonso et al. 2019). The angular power spectra in Eq. (39) are directly measured from maps so the contribution from noise is also included. We assume that the random noise in the map is Gaussian and independent of the signal.
The second term Cov^{T} is the connected term from the trispectrum, which is given by
where T_{uvwz}(k) is the trispectrum. Using the halo model, the trispectrum is decomposed into one to fourhalo terms. Schaan et al. (2018) shows that the onehalo term dominates the CIB trispectrum. As galaxies have a similar spatial distribution to the CIB, we only take the onehalo term into account for our CIBgalaxy crosscorrelation:
We will see that this term is negligible in the covariance matrix.
The third term Cov^{SSC} is called the super sample covariance (SSC; Takada & Hu 2013), which is the sample variance that arises from modes that are larger than the survey footprint. The SSC can dominate the covariance of power spectrum estimators for modes much smaller than the survey footprint, and includes contributions from halo sample variance, beat coupling, and their crosscorrelation. The SSC can also be calculated in the halo model framework (Lacasa et al. 2018; Osato & Takada 2021).
In this work, the nonGaussian covariance components Cov^{T} and Cov^{SSC} are calculated with the halo model formalism as implemented in the CCL package (Chisari et al. 2019)^{12}, and are then summed with Cov^{G} to get the full covariance matrix. Unlike Yan et al. (2021), who calculated covariance matrices independently for the different tomographic bins, the CIBgalaxy crosscorrelation in this work is highly correlated across galaxy tomographic bins and CIB frequency bands. Therefore, we calculate the whole covariance matrix of all the 15 crosscorrelations, giving a matrix with a sidelength of 5 × 3 × 10 = 150.
We note that, to calculate the analytical covariance matrix, we needed to use model parameters that we do not know a priori. So we utilised an iterative covariance estimate (Ando et al. 2017): we first took reasonable values for the model parameters, given by M21, to calculate these covariance matrices, and used them to constrain a set of bestfit model parameters. We then updated the covariance matrix with these bestfit parameters and fitted the model again. In practice, the constraint from the second step is always consistent with the first round, but we nonetheless took the constraints from the second step as our fiducial results.
Figure 4 shows the correlation coefficient matrix. The five diagonal golden blocks have very high offdiagonal terms, which means that the crosscorrelations between galaxies in the same tomographic bin with three CIB channels have a very high correlation (about 95%). This is because the CIB signals from different frequencies are essentially generated by the same galaxies. The correlation of offdiagonal golden blocks is weaker but still nonnegligible: this correlation is from the overlap of galaxy redshift distributions in different tomographic bins, as shown in Fig. 2. We also note that the SSC term contributes up to 8% of the total standard deviation, while the trispectrum term is insignificant (contributing < 0.1% to all covariance terms). This is in contrast to Tröster et al. (2021b), who study tSZ crosscorrelations and find that the SSC term was a more significant contributor to their covariance (contributing ∼20% to their offdiagonal covariance terms). The reason for this difference is that, compared to the tSZ effect, the galaxy distribution is more concentrated. This causes the nonGaussian term to remain insignificant until considerably smaller scales than the tSZ effect: beyond the scales probed in this study (ℓ> 2000).
Fig. 4. Correlation coefficient matrix of our crosscorrelation measurements. The colour scale is from 0 (black) to 1 (white). Each block enclosed by a white grid is the covariance between each pair of crosscorrelations indicated with ticks (bin p × ν GHz), while that enclosed by a golden grid corresponds to the covariance between the CIB crossgalaxies from each pair of tomographic bins. The matrix has nonzero elements at all cells, but the offdiagonal elements in each crosscorrelation are vanishingly small. 
Finally, an alternative estimation of the covariance matrix is shown in Appendix .
4.3. Systematics
4.3.1. CIB colourcorrection and calibration error
The flux given in the Planck sky maps follows the photometric convention that νI_{ν}=constant (Planck Collaboration XXX 2014). The flux therefore has to be colourcorrected for sources with a different SED. Therefore, the CIBgalaxy crosscorrelation should also be corrected as
where cc_{ν} is the colourcorrection factor at frequency ν. In this work, we adopt the colourcorrection factors given by Béthermin et al. (2012) with values of 1.097, 1.068, and 0.995 for the 353, 545, and 857 GHz bands, respectively.
Additionally, in Maniyar et al. (2018, 2021) the authors introduce a scaling factor as an additional calibration tool when working with the L19 CIB maps. However, they constrain this factor to be very close to one (at a level of ∼ ± 1%). As such, in this work we neglect the additional calibration factor.
4.3.2. Cosmic magnification
The measured galaxy overdensity depends not only on the real galaxy distribution but also on lensing magnification induced by the lineofsight mass distribution (Schneider 1989; Narayan 1989). This socalled cosmic magnification has two effects on the measured galaxy overdensity: (i) overdensities along the lineofsight cause the local angular separation between source galaxies to increase, so the galaxy spatial distributions are diluted and the crosscorrelation is suppressed; and (ii) lenses along the lineofsight magnify the flux of source galaxies such that fainter galaxies enter the observed sample, so the overdensity typically increases. These effects potentially bias galaxyrelated crosscorrelations, especially for highredshift galaxies (Hui et al. 2007; Ziour & Hui 2008; Hilbert et al. 2009). Cosmic magnification depends on the magnitude limit of the galaxy survey in question, and on the slope of the number counts of the sample under consideration at the magnitude limit. We follow Yan et al. (2021) to correct this effect, and note that this correction has a negligible impact on our bestfit CIB parameters.
4.3.3. Redshift distribution uncertainty
The SOMcalibrated galaxy redshift distributions have an uncertainty on their mean at a level of ∼0.01 Hildebrandt et al. (2021), which could affect galaxy crosscorrelations. To test the robustness of our results to this uncertainty, we run a further measurement including additional free parameters that allow for shifts in the mean of the redshift distributions. With these additional free parameters, the shifted galaxy redshift distributions are given by
where is the shifted galaxy redshift distribution in the ith tomographic bin, and δ_{z, i} is the shift of the redshift distribution parameter in the ith bin. The priors on δ_{z, i} are assumed to be covariant Gaussians centred at (i.e. the mean δ_{z, i}). Hildebrandt et al. (2021) evaluated both and the covariance matrix from simulations, but did so using the redshift distributions calculated with lensing weights. As previously discussed, however, the Φ_{g}(z) used in this work is from the SOM calibration without lensing weights. Therefore, neither the estimates of nor their covariance matrix from Hildebrandt et al. (2021) are formally correct for the Φ_{g}(z) in this work. To correctly estimate and the associated covariance matrix for this work, one needs to analyse another simulation suite for the SOM calibration without lensing weight, which is not currently available. However, given the similarity between the lensingweighted and unweighted redshift distributions (Fig. 2), we can alternatively adopt a conservative prior on δ_{z, i} with the mean values and uncertainties that are three times than the fiducial values given by Hildebrandt et al. (2021). Our fiducial δ_{z, i} covariance matrix is therefore simply defined as nine times the nominal KiDS1000 δ_{z, i} covariance matrix. This yields an absolute uncertainty at a level of 0.04, about two times the difference between the nominal KiDS1000 lensingweighted Φ_{g}(z) and the unweighted Φ_{g}(z) that we use in this work.
4.4. Likelihood
We constrained the model parameters in two ways. The first is crosscorrelation only (‘CIB × KiDS fitting’ hereafter), and the second is combining crosscorrelation and external ρ_{SFR} measurements (‘CIB × KiDS + ρ_{SFR} fitting’ hereafter). For the crosscorrelation only fitting, since we are working with a wide ℓ range, there are many degreesoffreedom in each ℓ bin. According to the central limit theorem, the binaveraged C_{ℓ}’s obey a Gaussian distribution around their true values. Thus, we assume that the measured power spectra follow a Gaussian likelihood:
where q stands for our model parameters given in Sect. 4. For the additional chain to test the robustness to redshift uncertainties, q also includes δ_{z, i} introduced in the previous subsection. The data vector is a concatenation of our 15 measured CIBgalaxy crosscorrelations; C(q) is the crosscorrelation predicted by the model described in Sect. 2 with parameters q.
The external ρ_{SFR} measurements are assumed to be independent of our crosscorrelation and similarly independent at different redshifts. Therefore, including them introduces an additional term in the likelihood:
where is the measured SFRD at the redshift z_{i}, while ρ_{SFR, i}(z_{i}q) is that predicted by our model (see Eq. (23)), and σ_{ρ, i} is the standard error of these SFRD measurements. We note that, we are still constraining the same HOD measurement as the crosscorrelationonly constraint.
We also perform fitting with only external ρ_{SFR} data and errors. This test serves as a consistency check between the CIB × KiDS measurement and previous multiwavelength studies, as well as a validation of our CIB halo model. The free parameters are the SFR parameters and HOD parameters (ten parameters in total).
We adopt the Markovchain MonteCarlo method to constrain our model parameters with the python package EMCEE (ForemanMackey et al. 2013). Best fit parameters are determined from the resulting chains, as being the sample with the smallest χ^{2} goodnessoffit. Marginal constraints on parameters, when quoted, are marginal means and standard deviations.
5. Results
5.1. Constraints on star formation rate density
We show our CIBKiDS crosscorrelation measurement in Fig. 5. Each panel presents the crosscorrelation between galaxies from a tomographic bin and CIB anisotropies in one frequency band. The points are the mean C_{ℓ} in each of the ten logarithmic ℓ bins. The error bars are the standard deviation calculated as the square root of the diagonal terms of the covariance matrix. We show the crosscorrelations calculated from the model given by Sect. 2 from the constrained parameters with CIB × KiDS fitting (see Table 3) as well as those calculated from constrained parameters given by the CIB × KiDS + ρ_{SFR} fitting. The reduced χ^{2} (RCS) for both fits are 1.14 and 1.10, with degreesoffreedom 125 and 142, respectively. In order to evaluate the goodnessoffit, we calculate the corresponding probabilitytoexceed (PTE) given the degreeoffreedom: our fits give PTE values of 0.13 and 0.2 for the CIB × KiDS and CIB × KiDS + ρ_{SFR} fitting, respectively. Heymans et al. (2021) adopts the criterion PTE > 0.001 (corresponding to a ∼3σ deviation) to be acceptable. Therefore, we conclude that our model generally fits the crosscorrelations well. We also notice that the fitting in lowredshift bins is worse than highredshift bins (see the ‘0.1 < z_{B} ≤ 0.3, 353 GHz’ panel in Fig. 5, for example), although correlation in the points makes ‘chibyeye’ inadvisable here.
Fig. 5. CIBgalaxy crosscorrelations with the five KiDS tomographic bins (rows) and the three CIB maps (columns). The grey points are measured from data, with standard deviation error bars calculated using the square root of the diagonal terms of the covariance matrix. The solid cyan lines show the bestfit crosscorrelation signals calculated using the CIBgalaxy crosscorrelation measurements alone, while the dashed blue lines show the bestfit crosscorrelations when jointly fitting the CIBgalaxy crosscorrelation measurements and the external SFRD. 
Summary of the prior ranges, the marginalised mean values, and the 68% credible intervals of SFR parameters.
We estimate the posterior with the last 100 000 points of each of our chains: CIB × KiDS, CIB × KiDS + ρ_{SFR}, and ρ_{SFR}only. The posterior of the four SFR parameters are shown in the triangle plot in Fig. 6. The distributions are marginalised over HOD parameters and shot noise amplitudes. In particular, we note the good constraint that our CIB × KiDS only results have over SFR parameters (cyan contours). The crosscorrelation only results are consistent with the results when analysing external SFRD data only (i.e. the red contours). This validates our CIBgalaxy crosscorrelation as a consistent yet independent probe of the cosmic SFRD, and further demonstrates the validity of the halo model (used in both studies) when applied to vastly different observational data. It is also encouraging that the crosscorrelation constraints are generally tighter than those made with the external SFRD data alone, demonstrating that the crosscorrelation approach provides a valuable tool for studying the cosmic SFR history. Our joint constraints are tighter still, demonstrating that there is different information in the two probes that can be leveraged for additional constraining power in their combination (the CIB × KiDS + ρ_{SFR} constraints shown in dark blue). The marginal parameter values and uncertainties are shown in Table 3, and are calculated as the mean and 68% confidence levels of the Gaussian kernel density estimate (KDE) constructed from the marginal posterior samples. The Gaussian kernel is dynamically adapted to the distribution of the marginal samples using the Botev et al. (2010) diffusion bandwidth estimator.
Fig. 6. Posterior of the SFR parameters. Contours show the 2D posteriors marginalised over all the other 25−2 = 23 parameters. The cyan contours show the constraints from the CIBKiDS crosscorrelation only, the dark blue contours show the constraints from a combination of the crosscorrelation and the external SFRD data, and the red contours show the constraints from the external SFRD data only. The dark and light regions in each contour show the 68% and 95% credible regions, respectively. The dashed lines show the bestfit model from M21. 
To evaluate the constraining power, we adopt the method from Asgari et al. (2021): we calculate the values of the marginalised posterior at both extremes of the prior distribution, and compare them with 0.135, the ratio between the peak of a Gaussian distribution and the height of the 2σ confidence level. If the posterior at the extreme is higher than 0.135, then the parameter boundary is not well constrained. We find that except the lower bound of M_{peak} for CIB × KiDS and the higher bound of sigma for ρ_{SFR}only, the other parameters are all constrained.
One of the parameters that is of particular interest is M_{peak}: the halo mass that houses the most efficient star formation activities. In Fig. 7, we summarise our results and recent observational results from the literature that have constrained this parameter. The three points above the dotted line are the constraints from this work. The other points are colourcoded according to their methods: the green point shows the result of SMG autocorrelations (Chen et al. 2016), the magenta point shows the measurement using LRGCIB crosscorrelations (Serra et al. 2014), and black points show measurements utilising CIB power spectra (Shang et al. 2012; Viero et al. 2013; Planck Collaboration XXX 2014; Maniyar et al. 2018, 2021). The purple band shows the 68% credible interval of our CIB × KiDS + ρ_{SFR} constraint. Except M21, our constraints are in agreement with previous studies within the 2σ level, but prefer a slightly lower M_{peak}. This may be due to the different data used in these studies, which would suggest that estimates of M_{peak} depend on galaxy types. Our results are in a mild tension with M21, which we hypothesis may be due to an inaccurate uncertainty estimate for their model parameters, driven by their assumption of a purely diagonal covariance matrix.
Fig. 7. Comparison of our constraints on M_{peak} with a number of recent results from the literature. The three points above the dotted line are the results from this work. The other points are colourcoded according to their methods: the green point shows the result from SMG autocorrelations (Chen et al. 2016), the magenta shows the measurements using LRGCIB crosscorrelations (Serra et al. 2014), and the black points show measurements using CIB power spectra (Shang et al. 2012; Viero et al. 2013; Planck Collaboration XXX 2014; Maniyar et al. 2018, 2021). The dark blue band shows the 68% credible interval of our CIB × KiDS + ρ_{SFR} marginal posterior constraint. 
It is interesting to compare our result with M21 because they also constrain SFRD history with the same halo model, but from CIB power spectra. According to Appendix A in M21, without introducing external SFRD data, CIB power spectra yield biased SFRD measurements at low redshift. Without the redshift information from SFRD measurements, the constraining power of model parameters are limited by degeneracy between them, which is reasonable because the CIB is a lineofsight integrated signal. In this regime, all parameters that describe the redshift distribution of CIB emission (M_{peak}, σ_{M, 0}, τ, z_{c}) should be degenerate. Therefore, it is remarkable that our CIB × KiDS constraints are able to constrain both σ_{M, 0} and τ. We attribute this increased precision to the use of tomography in our CIB × KiDS measurements.
We note that the crosscorrelationonly measurement does not constrain log_{10}M_{peak} well. This is because M_{peak} depends primarily on the CIB signal at high redshifts. We verify this by calculating the redshift distribution of mean CIB intensity defined in Eq. (16), while varying M_{peak} and fixing all other parameters. The result at 545 GHz is presented in Fig. 8; results using the other two channels show similar behaviour. It is clear that varying M_{peak} affects the CIB signal at high redshift more dramatically than at low redshift. In the redshift range of our galaxy sample, the mean CIB emission does not change significantly with M_{peak} as much as z > 1.5, especially for the lowest tomographic bins. Therefore, external SFRD measured at high redshifts, where the CIB intensity is more sensitive to M_{peak}, provides additional constraints on M_{peak}.
Fig. 8. CIB intensity at 545 GHz as a function of redshift while varying M_{peak} and keeping all other parameters fixed (solid lines, right x axis). We also plot the redshift distributions of the galaxy sample Φ_{g}(z) (dashed lines, left y axis). The CIB emissions are calculated from Eq. (22). The black line shows dI_{545}/dz, which corresponds to the CIB × KiDS bestfit parameters. 
The constraints of SFRD is shown in Fig. 9 with respect to lookback time in panel a and redshift in panel b. The CIB × KiDS fit is shown in cyan line while the CIB × KiDS + ρ_{SFR} in dark blue. We estimate the credible interval of our fits by calculating SFRD using 10 000 realisations of our model parameters, drawn from the posterior distribution shown in Fig. 6. The 68% credible interval is shown as bands with corresponding colours, which are calculated by computing the 10 000 samples from the model and deriving the SFRD at a range of redshifts. Credible intervals are computed at each redshift using a Gaussian KDE, and these constraints are connected to form the filled regions. The magenta and green lines are the bestfit SFRD from M21 and Davies et al. (2016), and the points with error bars are SFRD estimates from previous studies (which are included in our CIB × KiDS + ρ_{SFR} analysis). The SFRD fitting given by M21 is from a combination of the CIB auto and crosscorrelation power spectra and external SFRD measurements, while Davies et al. (2016) is an empirical model estimated using galaxy surveys. These two figures demonstrate that our measurements agree well with previous studies using different analysis techniques. The CIB × KiDS + ρ_{SFR} fitting agrees well with external SFRD in all redshift ranges. Notably, we are also able to obtain fairly accurate and precise constrain on SFRD up to z ∼ 1.5 (corresponding to a lookback time of 10 Gyr) using our CIB × KiDS crosscorrelations alone. Beyond z ∼ 1.5, CIB × KiDS fitting yields large uncertainties because our KiDS sample has few galaxies beyond this point. The constraint of SFRD drops below 3σ level beyond z ∼ 1.8 (see the dotted lines in Fig. 9b) As a result, our sample is not deep enough to directly constrain M_{peak}. We conclude that the CIB × KiDS constraint yields a peak SFRD of yr^{−1} Mpc^{−3} at , corresponding to a lookback time of Gyr, while CIB × KiDS + ρ_{SFR} constraint gives a peak SFRD of yr^{−1} Mpc^{−3} at , corresponding to a lookback time of Gyr, consistent with previous observations of the cosmic star formation history.
Fig. 9. Evolution of the SFRD with respect to lookback time (panel a) and redshift (panel b). The SFRD calculated from this work is presented as cyan (crosscorrelation only) and dark blue (crosscorrelation plus external SFRD) lines and shaded regions. The shaded regions enclose the 1σ credible region of the fits and are calculated from 10 000 realisations of SFR parameters from the posterior distribution. The 3σ credible region of the crosscorrelationonly SFRD is also shown in panel b with dotted cyan lines. We note that the lower 3σ limit crosses zero at z ∼ 1.8. The magenta and green lines are the bestfit SFRD from two previous studies, and the points with error bars are the SFRD from previous studies (which are included in our CIB × KiDS + ρ_{SFR} analysis). 
In parallel with observations, simulations and semianalytical models (SAMs) also give estimations on SFRD. In order to check the consistency between observation, simulations and SAMs, we compare our constraints on the SFRD with results given by simulations and SAMs in Fig. 10. The results are from the GALFORM SAM (Guo et al. 2016, with the GonzalezPerez et al. 2014 version; purple line), EAGLE (Guo et al. 2016, with the Schaye et al. 2015 simulation; khaki line), and the LGALAXIES SAM (Henriques et al. 2015, red line). These models adopt different simplifications for active galactic nucleus (AGN) and supernova feedbacks, different star formation thresholds, and different gascooling models. We see that EAGLE predicts a slightly lower SFRD (at essentially all times) than is predicted from our results. GALFORM, on the otherhand, agrees well with our CIB × KiDS fits at intermediatetoearly times, but predicts a higher SFRD than our model in the 1−5 Gyr range. As is discussed in Driver et al. (2017), this might be due to the high dust abundance at low redshift in the GonzalezPerez et al. (2014) version of GALFORM. LGALAXIES incorporates different methodologies to model the environment, cooling processes, and star formation of galaxies. Combining these modifications ensures that massive galaxies are predominately quenched at earlier times while low mass galaxies have star formation histories that are more extended into the latetime Universe. As can be seen, this results in a low SFRD prediction at intermediate redshifts compared to our data, but a better fit at low redshift. However, given the large uncertainties on our CIB × KiDS fits, the CIBgalaxy crosscorrelation is currently not precise enough to invalidate any of these simulations. Comparing to our CIB × KiDS + ρ_{SFR} analysis, none of the three simulations are able to reproduce our fits at all redshifts. This highlights the complexity of the underlying physics that determines the form of the observed SFRD.
Fig. 10. Evolution of the SFRD with respect to lookback time from this work (see Fig. 9), compared to simulations and SAMs. 
We also present the posterior of SFR parameters with varying δ_{z} in Fig. 11. The blue contour is our fiducial posterior, while the red contour is the posterior of the SFR parameters when allowing free variation of δ_{z} described in Sect. 4. Varying δ_{z} only slightly loosens the constraints on η_{max}, while all other posteriors are largely unchanged. The posterior distributions of our δ_{z, i} parameters follow the prior, demonstrating that there is no redshiftdistribution selfcalibration occurring in this analysis. Nonetheless, given the conservative choices made with our δ_{z, i} priors, we conclude that our constraints are robust to redshift distribution uncertainties in the galaxy sample. This result is largely expected, though, as CIB emission varies only slightly within the level of the δ_{z} uncertainties (see Fig. 2).
Fig. 11. Posterior of SFR parameters from the CIB × KiDS fit with fixed (blue) and freely varying (red) priors on the means of Φ_{g}(z). The blue contour is our fiducial posterior. 
5.2. Constraint on galaxy bias
In this subsection we present the constraints on HOD parameters described in Sect. 2 and the derived galaxy bias. This is not the main scope of this paper, but it is nonetheless an interesting outcome of our study. Galaxy bias parameters are typically constrained using galaxy power spectra; however, this is challenging to perform with KiDS (and indeed with any fulldepth galaxy survey) due to the complicated (artificial) variable depth between survey fields introduced by variable observing conditions. Yan et al. (2021) constrained the linear galaxy bias for KiDS using galaxyCMB lensing crosscorrelations, assuming a linear model. In this work, we derive the linear galaxy bias from the constrained HOD parameters.
The scaledependent galaxy bias is defined as
where the galaxy density fluctuation δ_{g} is the Fourier transform of Eq. (10). The linear galaxy bias is given by
The constrained HOD parameters are presented in Table 4. Similar to the calculation of the bestfit SFRD and its uncertainty, we calculate the bestfit and 1σ credible interval of and present it in Fig. 12. The resulting linear galaxy bias increases from at z = 0 to at z = 1.5. We also overplot constraints from previous studies on galaxy bias of starforming galaxies. The magenta line shows the bestfit model from Maniyar et al. (2018); the green line shows the bestfit ‘starforming’ model from Durkalec et al. (2015); the red points are the derived galaxy bias of starforming galaxies from Cochrane et al. (2017). We find good agreement between our result and these studies. The evolutionary trend of galaxy bias measured in this work is also in agreement with Lagache et al. (2007). It should be noted that the constraint of galaxy bias worsens at high redshift, as our galaxy sample is limited to z < 1.5.
Fig. 12. Linear galaxy bias constrained from CIBgalaxy crosscorrelation. The solid blue line shows the bestfit , and the band with dotteddashed boundary shows the upper and lower 1σ errors. We also overplot results from previous studies. 
Summary of the prior range, the mean, and 68% confidence level of the HOD parameters.
The galaxy bias constrained from CIBgalaxy crosscorrelation is slightly higher than that constrained from galaxyCMB lensing given by Yan et al. (2021) (see the orange points). Galaxy bias from this work also shows a stronger evolution. These might be due to the fact that in Yan et al. (2021) all the galaxies in the KiDS gold sample contribute to the constraint, while the CIBgalaxy crosscorrelation in this work is mainly sensitive to KiDS galaxies that are active in star formation. The fact that CIB crosscorrelation gives higher galaxy bias means that these starforming galaxies are (on average) more clustered than the galaxies detected by optical survey, especially at high redshift, which calls for further study.
6. Conclusions
In this work we measure the tomographic crosscorrelation between the KiDS1000 gold galaxy sample and the CIB maps made from Planck data. The motivation of this work is to provide a novel measurement of cosmic star formation history. We summarise our main conclusions in this section.

The crosscorrelation has a significance of 43σ, which is impressive given that the CIB signal is relatively low in the redshift range of our galaxy sample. Our CIB model yields an RCS value of 1.14 when fitting to the KiDSCIB crosscorrelation. Given the degreesoffreedom, this corresponds to a PTE value of 0.13, meaning that the model fits the data well.

The constraints on the SFRD parameters from crosscorrelation agree with those measured from external data, demonstrating that crosscorrelation provides a novel, independent, and consistent probe of the cosmic star formation history. Moreover, this indicates that the halo model proposed by M21 is valid for both CIB crosscorrelations and multiwavelength studies of galaxies.

With our crosscorrelation measurement, the maximum star formation efficiency is , in agreement with M21. Our CIB × KiDSonly measurement is unable to yield tight constraints on M_{peak}, the halo mass that hosts the highest star formation efficiency, due to the galaxy sample being limited to z < 1.5. A combination of crosscorrelation and external SFRD measurements (the CIB × KiDS + ρ_{SFR} constraints), however, tightens the constraint on log_{10}M_{peak} to , in agreement with previous studies within the 2σ level (albeit with a preference for a slightly lower M_{peak} at posterior maximum). This may be due to the different data used in this study, which would imply that measurements of M_{peak} are dependent on galaxy types. We leave an exploration of this for future work. Moreover, the bestfit M_{peak} from both CIB × KiDS and CIB × KiDS + ρ_{SFR} are in mild tension with M21, which calls for further investigation.

We derived the SFRD from our posterior constraints over our various model parameters. The CIB × KiDS constraint on the SFRD history is poor at high redshift because of sky coverage and depth limitations. Nevertheless, we note that in the redshift range probed by our sample of KiDS galaxies (z < 1.5, corresponding to a lookback time of ∼10 Gyr), crosscorrelations give an SFRD that is consistent with previous galaxybased multiwavelength studies, at comparable constraining power. The CIB × KiDS + ρ_{SFR} results tighten the constraint on the SFRD at all redshifts and are consistent with the SFRD from CIB × KiDS and previous studies. We also compare the SFRD from this work with simulations and SAMs and find that our CIB × KiDS constraints have the same trend as all the simulated SFRDs; however, our results are not sufficiently precise to invalidate any one model. Moreover, none of the simulations agrees with our CIB × KiDS + ρ_{SFR} constraint at all the times, highlighting the complexity of the physical processes that underpin star formation activity in galaxies.

We also constrain the linear galaxy bias for KiDS galaxies that have significant IR emission. As for the SFRD, we can only constrain galaxy bias below z ∼ 1.5, and with about 25% precision. The constraint is limited by both the sky coverage and angular resolution of the CIB map. We also note that we model the redshift dependence of HOD parameters as a simple linear model with respect to the scale factor a, which could be improved for future studies. We derived the linear galaxy bias from the constrained HOD parameters, yielding an increasing galaxy bias greater than one. The evolution of galaxy bias constrained in this work agrees with Durkalec et al. (2015), Cochrane et al. (2017), and Maniyar et al. (2018).

For systematics, we took the colour correction of CIB flux and the effects of cosmic magnification into account. We also investigated the robustness of our results to the uncertainties in redshift distributions by allowing for shifts in the redshift distribution means. We find that this does not affect our constraints and conclude that our results are robust to the uncertainty in redshift distribution calibration. However, for future studies with higher signaltonoise, this may become important.
7. Discussions and future prospects
In this work we adopt the halo model for the CIB from M21, with a minor modification such that the HOD model is consistent with Zheng et al. (2005). This model includes the information on the dust SED, star formation history, and galaxy abundance. Compared with other CIB models such as Shang et al. (2012) and Cao et al. (2020), it clearly formulates the redshift and halo mass dependence of SFR, which allows us to constrain the cosmic star formation history from crosscorrelations between CIB and largescale structure surveys.
We make several assumptions to simplify the CIB halo model. For example, we assume that the SFRIR luminosity connection can be described by a single Kennicutt parameter in the IR bands (Kennicutt 1998), with the assumption of a Chabrier IMF (Chabrier 2003). The SFR is modelled as the star formation efficiency times the baryonic accretion rate, which can be alternatively modelled by treating IR and UV components separately and introducing the stellarmass function (Béthermin et al. 2013). In addition, we have assumed that the SFR has the same mass dependence for central and satellite galaxies. In this work, we take the SED model from Béthermin et al. (2015), which does not include the AGN contribution that could contaminate the midIR signal at z > 2 (Béthermin et al. 2017). This is beyond the scope of this paper, but future studies at these redshifts and frequencies should consider the influence of such assumptions. We do not discuss the thermodynamic properties of extragalactic dust, which are encoded in the SED. Shang et al. (2012) provide an empirical model of the SED but do not model the SFR explicitly. Redshiftdependent SED changes and the evolution of the SFR are degenerate in the CIB signal, which might be resolved by introducing an additional tracer, or introducing the CIB at more frequencies (e.g. from the Herschel and Spitzer telescopes). We also note that the fit is worse at low redshift, which may indicate a limitation of our simplified model, or indicate the inaccurate measurement of our SED at low redshift. Finally, we have fixed all the cosmological parameters in this study. More sophisticated modelling is left for future studies that will use larger datasets.
The KiDS galaxy sample in this study has the advantage of extensive survey depth. Although our sample is not deep enough to directly constrain M_{peak}, it yields a sensible measurement of the star formation history out to z = 1.5. From Fig. 2, we see that the redshift distribution of the KiDS galaxies is mostly in the rising regime of the CIB signal, which peaks at z ∼ 1.5. For future galaxy surveys that will be able to reach z ∼ 3, such as the Rubin Observatory LSST and Euclid, one would expect a more pronounced CIBgalaxy crosscorrelation signal. In this context, we perform a forecast on the constraining power of the ongoing CFIS and the future LSST survey in Appendix B. We find that the improvement of sky coverage makes CFIS yield a similar constraining power as our CIB × KiDS + ρ_{SFR} constraint, while the more extensive LSST depth makes it possible to tightly constrain all the SFR parameters. The Dark Energy Survey has published its year3 data, and CFIS will provide data with larger sky coverage in the near future, which will certainly help in clarifying the statistical significance of our CIB × KiDS result without adding the external SFRD measurements. Furthermore, the LSSTCIB crosscorrelation will be a promising tool, yielding enough constraining power to validate different SFR models and give us insight into the physics underpinning the cosmic star formation history.
Other prospective CIB studies include crosscorrelating the CIB with galaxy samples as a function of brightness, colour, stellar mass, or morphology, as it is well known that star formation depends heavily on these properties (as examples, Yajima et al. 2012, for brightness dependence; Mahajan & Raychaudhury 2009 for colour dependence; Kusakabe et al. 2018 for stellar mass dependence; and Tojeiro et al. 2013 for colourmorphology combination dependence). In addition, dust properties should also depend on these properties (Wolf et al. 2009; Noll & Pierini 2005). A CIB crosscorrelation with different types of galaxies could therefore serve as a new independent probe of these properties.
Central and satellite galaxies are located in different physical environments, resulting in different SFRs. Specifically, AGN feedback and quenching are two interesting processes that can effect a galaxy’s SFR, and they are found to be different in central and satellite galaxies (Wang et al. 2018). In this work we do not separately study the SFR of central and satellite galaxies because the total SFR is dominated by central galaxies with higher halo mass, and the SFR of satellite galaxies is not well constrained. This may be improved with future surveys. Once more significant crosscorrelation measurements are available, we will be able to study quenching and AGN activities in central and satellite galaxies by adding more parameters that describe these effects to the SFR model.
The constraints on HOD parameters suggest an increasing galaxy bias through redshift. At high redshift, the high galaxy bias might be due to the triggering of star formation by mergers in dense environments at z ∼ 1 (Wetzel et al. 2007). At low redshift, however, star formation is quenched through gas stripping in dense regions (Postman et al. 2005), leading to a lower overall bias. In conclusion, a comparison between the galaxy bias of normal galaxies and starforming galaxies indicates the evolution of the link between star formation and environments. It should be noted that the constraining power of HOD parameters in this study is weak because of the limitation of both sky coverage and depth. Moreover, the linear formalism of HOD parameters is a simplified empirical model. Future studies with improved sky coverage and depth should improve the constraints of .
In summary, the CIB is a gold mine that encodes information about the star formation history of the Universe, extragalactic dust properties, and galaxy abundance. This work validates the CIBgalaxy crosscorrelation method as a valuable tool for understanding the cosmic star formation history. The success of our measurement here, despite the limitations discussed above, provides an exceedingly positive outlook for future analyses of the CIBgalaxy crosscorrelation. Larger, deeper datasets, coupled with more complex sample subdivisions, will allow us to leverage CIB crosscorrelations to better understand the growth and evolution of galaxies in our Universe.
The definition of the radial kernel can be quite arbitrary since, in practice, it can absorb any factors in the 3D field term that only depend on redshift. In the following subsections, we define the galaxy and the CIB kernels to emphasise the physical meaning of the 3D fields studied in this work.
Acknowledgments
We would like to thank Yunhao Zhang, Dr. Gary Hinshaw, and Dr. Joachim HarnoisDéraps for useful discussions. We thank the Planck Collaboration for the allsky data available at https://www.cosmos.esa.int/web/planck/pla and Dr. Daniel Lenz for making the CIB maps. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A3016, 177.A3017, 177.A3018 and 179.A2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWOM grants; Target; the University of Padova, and the University Federico II (Naples). Z.Y. acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max PlanckHumboldt Research Award endowed by the Federal Ministry of Education and Research (Germany). L.V.W. and S.G. acknowledge acknowledge the support support by the University of British Columbia, Canada’s NSERC, and CIFAR. M.B. is supported by the Polish National Science Center through grants no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698, 2018/31/G/ST9/03388 and 2020/39/B/ST9/03494, and by the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12. H.H. is supported by a Heisenberg grant of the Deutsche Forschungsgemeinschaft (Hi 1495/51) and, along with A.H. Wright, by an ERC Consolidator Grant (No. 770935). T.T. acknowledges support from the Leverhulme Trust The data in this paper is analysed with opensource python packages NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), ASTROPY (Astropy Collaboration 2018), MATPLOTLIB (Hunter 2007), HEALPY (Zonca et al. 2019), NAMASTER (Alonso et al. 2019), CCL (Chisari et al. 2019), EMCEE (ForemanMackey et al. 2013), and GETDIST (Lewis 2019). We also use WEBPLOTDIGITIZER (Rohatgi 2021) to digitize some external data from plots in the literature. Author contributions. All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (ZY & LvW) followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products. The second group covers those who have either made a significant contribution to the data products, or to the scientific analysis. Data availability. The crosscorrelation data and the MCMC chain of this work will be shared on reasonable request to the corresponding author.
References
 Abbott, T., Abdalla, F. B., Aleksić, J., et al. 2016, MNRAS, 460, 1270 [Google Scholar]
 Alonso, D., Sanchez, J., & Slosar, A. 2019, MNRAS, 484, 4127 [Google Scholar]
 Ando, S., BenoitLévy, A., & Komatsu, E. 2017, MNRAS, 473, 4318 [Google Scholar]
 Asgari, M., Lin, C.A., Joachimi, B., et al. 2021, A&A, 645, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Berta, S., Magnelli, B., Lutz, D., et al. 2010, A&A, 518, L30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23 [Google Scholar]
 Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [Google Scholar]
 Béthermin, M., Wu, H.Y., Lagache, G., et al. 2017, A&A, 607, A89 [Google Scholar]
 Botev, Z. I., Grotowski, J. F., & Kroese, D. P. 2010, Ann. Stat., 38, 2916 [Google Scholar]
 Cao, Y., Gong, Y., Feng, C., et al. 2020, ApJ, 901, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
 Chen, C.C., Smail, I., Swinbank, A. M., et al. 2016, ApJ, 831, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Chiang, Y.K., Makiya, R., Ménard, B., & Komatsu, E. 2020, ApJ, 902, 56 [Google Scholar]
 Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2 [Google Scholar]
 Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195 [NASA ADS] [CrossRef] [Google Scholar]
 Cochrane, R. K., Best, P. N., Sobral, D., et al. 2017, MNRAS, 469, 2913 [NASA ADS] [CrossRef] [Google Scholar]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [Google Scholar]
 Coupon, J., Kilbinger, M., McCracken, H. J., et al. 2012, A&A, 542, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Davies, L. J. M., Driver, S. P., Robotham, A. S. G., et al. 2016, MNRAS, 461, 458 [NASA ADS] [CrossRef] [Google Scholar]
 de Jong, J. T., Kleijn, G. A. V., Kuijken, K. H., et al. 2013, Exp. Astron., 35, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Dole, H., Lagache, G., Puget, J.L., et al. 2006, A&A, 451, 417 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2017, MNRAS, 475, 2891 [Google Scholar]
 Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64 [Google Scholar]
 Durkalec, A., Le Fèvre, O., Pollo, A., et al. 2015, A&A, 583, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1998, ApJ, 508, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32 [NASA ADS] [Google Scholar]
 Efstathiou, G. 2004, MNRAS, 349, 603 [Google Scholar]
 Fakhouri, O., Ma, C.P., & BoylanKolchin, M. 2010, MNRAS, 406, 2267 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Galametz, M., Albrecht, M., Kennicutt, R., et al. 2014, MNRAS, 439, 2542 [NASA ADS] [CrossRef] [Google Scholar]
 GarcíaGarcía, C., Alonso, D., & Bellini, E. 2019, JCAP, 2019, 043 [Google Scholar]
 Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GonzalezPerez, V., Lacey, C. G., Baugh, C. M., et al. 2014, MNRAS, 439, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
 Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23 [Google Scholar]
 Guo, Q., GonzalezPerez, V., Guo, Q., et al. 2016, MNRAS, 461, 3457 [Google Scholar]
 Hall, K. R., Crichton, D., Marriage, T., Zakamska, N. L., & Mandelbaum, R. 2018, MNRAS, 480, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Hang, Q., Alam, S., Peacock, J. A., & Cai, Y.C. 2021, MNRAS, 501, 1481 [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663 [Google Scholar]
 Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hilbert, S., Hartlap, J., White, S., & Schneider, P. 2009, A&A, 499, 31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124 [EDP Sciences] [Google Scholar]
 Hojjati, A., Tröster, T., HarnoisDéraps, J., et al. 2017, MNRAS, 471, 1565 [Google Scholar]
 Hui, L., Gaztanaga, E., & LoVerde, M. 2007, Phys. Rev. D, 76, 103502 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Ibata, R. A., McConnachie, A., Cuillandre, J.C., et al. 2017, ApJ, 848, 128 [Google Scholar]
 Ishikawa, S., Kashikawa, N., Tanaka, M., et al. 2020, ApJ, 904, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N. 1992, ApJ, 388, 272 [Google Scholar]
 Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189 [Google Scholar]
 Koukoufilippas, N., Alonso, D., Bilicki, M., & Peacock, J. A. 2020, MNRAS, 491, 5464 [Google Scholar]
 Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuntz, A. 2015, A&A, 584, A53 [EDP Sciences] [Google Scholar]
 Kusakabe, H., Shimasaku, K., Ouchi, M., et al. 2018, PASJ, 70, 4 [Google Scholar]
 Lacasa, F., Lima, M., & Aguena, M. 2018, A&A, 611, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lagache, G., Bavouzet, N., FernandezConde, N., et al. 2007, ApJ, 665, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Laureijs, R. J., Duvet, L., Sanz, I. E., et al. 2010, Int. Soc. Opt. Photon., 7731, 77311H [Google Scholar]
 Lenz, D., Doré, O., & Lagache, G. 2019, ApJ, 883, 75 [Google Scholar]
 Lewis, A. 2019, ArXiv eprints [arXiv:1910.13970] [Google Scholar]
 Lidz, A., Zahn, O., Furlanetto, S. R., et al. 2008, ApJ, 690, 252 [Google Scholar]
 Limber, D. N. 1953, ApJ, 117, 134 [NASA ADS] [CrossRef] [Google Scholar]
 LSST Science Collaboration 2009, LSST Science Book, Version 2.0 [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
 Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mahajan, S., & Raychaudhury, S. 2009, MNRAS, 400, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Mak, D. S. Y., Challinor, A., Efstathiou, G., & Lagache, G. 2017, MNRAS, 466, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Maniyar, A. S., Béthermin, M., & Lagache, G. 2018, A&A, 614, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maniyar, A., Lagache, G., Béthermin, M., & Ilić, S. 2019, A&A, 621, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maniyar, A., Béthermin, M., & Lagache, G. 2021, A&A, 645, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchetti, L., Vaccari, M., Franceschini, A., et al. 2016, MNRAS, 456, 1999 [Google Scholar]
 Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401 [Google Scholar]
 Miller, L., Heymans, C., Kitching, T., et al. 2013, MNRAS, 429, 2858 [Google Scholar]
 Narayan, R. 1989, ApJ, 339, L53 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Noll, S., & Pierini, D. 2005, ArXiv eprints [arXiv:astroph/0508549] [Google Scholar]
 Opher, R. 1998, in 19th Texas Symposium on Relativistic Astrophysics and Cosmology, eds. J. Paul, T. Montmerle, & E. Aubourg, 533 [Google Scholar]
 Osato, K., & Takada, M. 2021, Phys. Rev. D, 103, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Federrath, C., Chabrier, G., et al. 2014, Protostars and Planets VI (Tucson: University of Arizona Press) [Google Scholar]
 Pandey, S., Baxter, E. J., Xu, Z., et al. 2019, Phys. Rev. D, 100, 063519 [Google Scholar]
 Partridge, R. B., & Peebles, P. J. E. 1967, ApJ, 147, 868 [Google Scholar]
 Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144 [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Postman, M., Franx, M., Cross, N. J. G., et al. 2005, ApJ, 623, 721 [Google Scholar]
 Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, A&A, 656, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reischke, R., Desjacques, V., & Zaroubi, S. 2020, MNRAS, 491, 1079 [NASA ADS] [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J.F. 2011, MNRAS, 418, 467 [Google Scholar]
 Rieke, G. H., AlonsoHerrero, A., Weiner, B. J., et al. 2009, ApJ, 692, 556 [NASA ADS] [CrossRef] [Google Scholar]
 Rohatgi, A. 2021, Webplotdigitizer: Version 4.5, https://automeris.io/WebPlotDigitizer [Google Scholar]
 Salvato, M., Ilbert, O., & Hoyle, B. 2019, Nat. Astron., 3, 212 [NASA ADS] [CrossRef] [Google Scholar]
 Schaan, E., Ferraro, S., & Spergel, D. N. 2018, Phys. Rev. D, 97, 123539 [NASA ADS] [CrossRef] [Google Scholar]
 Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
 Schmidt, S. J., Ménard, B., Scranton, R., et al. 2014, MNRAS, 446, 2696 [Google Scholar]
 Schneider, P. 1989, A&A, 221, 221 [NASA ADS] [Google Scholar]
 Seljak, U. 2000, MNRAS, 318, 203 [Google Scholar]
 Serra, P., Lagache, G., Doré, O., Pullen, A., & White, M. 2014, A&A, 570, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shang, C., Haiman, Z., Knox, L., & Oh, S. P. 2012, MNRAS, 421, 2832 [NASA ADS] [CrossRef] [Google Scholar]
 Singh, S., Mandelbaum, R., & Brownstein, J. R. 2017, MNRAS, 464, 2120 [Google Scholar]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [Google Scholar]
 Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007, ApJ, 656, 770 [Google Scholar]
 Spitzer, I. 2022, PhD Thesis, University of Waterloo, Canada [Google Scholar]
 Sunyaev, R., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173 [Google Scholar]
 Takada, M., & Hu, W. 2013, Phys. Rev. D, 87, 123504 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J. L., & Wetzel, A. R. 2010, ApJ, 719, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [Google Scholar]
 Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878 [NASA ADS] [CrossRef] [Google Scholar]
 Tinsley, B. M. 1980, Fund. Cosm. Phys., 5, 287 [Google Scholar]
 Tojeiro, R., Masters, K. L., Richards, J., et al. 2013, MNRAS, 432, 359 [Google Scholar]
 Tröster, T., Asgari, M., Blake, C., et al. 2021a, A&A, 649, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tröster, T., Mead, A. J., Heymans, C., et al. 2021b, A&A, 660, A27 [Google Scholar]
 van den Bosch, F. C., More, S., Cacciato, M., Mo, H., & Yang, X. 2013, MNRAS, 430, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Viero, M. P., Wang, L., Zemcov, M., et al. 2013, ApJ, 772, 77 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
 Wang, L., Viero, M., Ross, N. P., et al. 2015, MNRAS, 449, 4476 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, E., Wang, H., Mo, H., et al. 2018, ApJ, 860, 102 [Google Scholar]
 Wetzel, A. R., Cohn, J. D., White, M., Holz, D. E., & Warren, M. S. 2007, ApJ, 656, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Wolf, C., AragónSalamanca, A., Balogh, M., et al. 2009, MNRAS, 393, 1302 [Google Scholar]
 Wolz, L., Murray, S. G., Blake, C., & Wyithe, J. S. 2018, MNRAS, 484, 1007 [Google Scholar]
 Wright, A. H., Hildebrandt, H., Kuijken, K., et al. 2019, A&A, 632, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wright, A. H., Hildebrandt, H., van den Busch, J. L., & Heymans, C. 2020, A&A, 637, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yajima, H., Umemura, M., & Mori, M. 2012, MNRAS, 420, 3381 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, Z., Hojjati, A., Tröster, T., Hinshaw, G., & Waerbeke, L. V. 2019, ApJ, 884, 139 [Google Scholar]
 Yan, Z., van Waerbeke, L., Tröster, T., et al. 2021, A&A, 651, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Ziour, R., & Hui, L. 2008, Phys. Rev. D, 78, 123517 [Google Scholar]
 Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Sour. Softw., 4, 1298 [Google Scholar]
Appendix A: The jackknife covariance matrix
An alternative method to estimate the covariance matrix is to use jackknife resampling, as is used in Yan et al. (2021). The idea is to resample the signal by removing a small patch of the sky and calculate the crosscorrelation from the remaining sky. We remove one patch at a time with replacement to get several crosscorrelation signals. The covariance matrix is estimated by calculating the covariance matrix from these resampled crosscorrelations. Except for assuming that all the patches are independent, this method does not depend on any other prior knowledge of the signal and is in principal able to take all the contributions into the account. However, there are two drawbacks of this method. Firstly, the KiDS footprint is small, so we can only have a small size of the jackknife sample. If we increase the sample by using small jackknife patches, we will lose accuracy in covariance corresponding to largescale modes. This problem is more severe in this study since we only use nearly half of the KiDS field. Secondly, removing patches from the sky changes the shape of the footprint, which will change the coupling matrix of crosscorrelations, thus biasing the covariance matrix. Therefore, we only use jackknife covariance as a consistency check of our analytical covariance matrix in this work. We generate 220 jackknife samples by removing HEALPIX pixels with NSIDE = 32, which correspond to a size of 1.8 degree. A comparison between standard deviations calculated from the analytical and jackknife covariance matrices is shown in Fig. A.1 from which we see that both covariance matrices agree well, with small deviations in large scales. We note that this scale (ℓ ∼ 100) is close to the size of removed pixels, so this deviation might be due to the poor sampling of large scales in our jackknife scheme.
Fig. A.1. Comparison between the standard deviation of obtained from the diagonal terms of the analytical (the beige lines) and jackknife (the dashed dark blue lines) covariance matrices. Each panel corresponds to one crosscorrelation signal. The three columns represent three CIB channels in increasing order from left to right; the five rows represent five KiDS tomographic bins in increasing order from top to bottom. 
Appendix B: Forecast for future surveys
As has been discussed in the main text, the SFRD constraint from CIBgalaxy crosscorrelation in this study is mainly limited by sky coverage and survey depth. Fortunately, future surveys will have much wider coverage and higher redshift depth. In this section we forecast the constraining power on the SFRD parameters from CIBgalaxy crosscorrelations with galaxies from future surveys. We adopt Fisher forecast by constructing the Fisher matrix:
where C(q) is again the model crosscorrelation with parameters q. The covariance matrix is calculated the same way as discussed in Sect. 4.2. Several factors in the covariance matrix are surveyspecified, namely sky coverage fraction, shot noises, and galaxy redshift distributions. The sky coverage fraction and galaxy redshift distributions can generally be found in survey proposal papers (which is to be introduced below). Shot noise (needed in Gaussian covariance matrix) for galaxygalaxy power spectra is calculated via
where is the mean number of galaxies per steradian, which can also be found in the survey proposal, and we estimate the galaxy density in each tomographic bin by weighting the total galaxy density by redshift distributions of each bin. The CIBCIB shot noise is taken from the constraints given by M21. CIBgalaxy shot noise are estimated as the bestfit shot noise corresponding to galaxy sample in the similar redshift bin from this work. For CIBgalaxy shot noise, one would expect that it should also be positively correlated with . For future surveys, this should be lower. However, without precise knowledge of this dependence, we use the CIBgalaxy shot noise for the KiDS survey as a conservative estimation.
We consider two ongoing and future galaxy surveys. The first is the LSST based at Vera Rubin Observatory (LSST Science Collaboration 2009). It is one of the major fourthgeneration cosmological surveys that cover about 20 000 deg^{2}. The overlap of LSST and the L19 CIB map is about 8% of the sky. The galaxy density is estimated to be 55.5 arcmin^{−2}. We adopt the overall redshift distribution and tomographic binning model from LSST Science Collaboration (2009). We take ten tomographic bins from z = 0 to z = 3.
The second is the CFIS (Ibata et al. 2017), an ongoing cosmological survey based at the CanadaFranceHawaii Telescope. When completed, it will yield a sky coverage of 3500 deg^{2}, overlapping with the L19 CIB map about 4% of the sky. The average galaxy density is about 10 arcmin^{−2}. The total redshift distribution is given by Spitzer (2022). We adopt the same tomographic binning model as LSST with five tomographic bins from z = 0 to z = 1.5.
We show the Fisher forecast contours in Fig. B.1. The green and the dark yellow contours correspond to CIBgalaxy crosscorrelations with future CFIS and LSST galaxy catalogues, respectively. The contours are 1σ and 2σ levels of confidence. We also plot the ‘forecast’ for our CIB × KiDS + ρ_{SFR} measurement with solid purple contours and the constraints from this work with dashed purple contours. We see that LSST yields much better constraints of all the SFR parameters while CFIS yields comparable constraints as CIB × KiDS + ρ_{SFR} (therefore both better than CIB × KiDS constraint).
Fig. B.1. Fisher forecast for CIBgalaxy crosscorrelation with LSST (green) and CFIS (dark yellow), as well as the ‘forecast’ for the CIB × KiDS + ρ_{SFR} measurement in this work (solid purple contour). The CIB × KiDS + ρ_{SFR} constraint is also shown as a dashed purple contour. Contours are the 1 and 2σ levels. 
We note that the maximum star formation efficiency η_{max} acts as an overall normalisation parameter, so its constraint is mainly related to sky coverage. The rest three parameters control the redshift dependence of SFRD (and thus the CIB emission). Therefore, a deeper survey would give a better constrain on them, especially M_{peak} (see Fig. 8). The constraints on these three parameters are not significantly improved with CFIS because CFIS has a similar depth as KiDS. LSST goes much deeper, so it can provide much better constraints.
The noise levels and redshift distributions in this forecast study are oversimplified. Moreover, the difference between the solid and the dashed purple contours shows that the posteriors are not Gaussian; therefore, the Fisher forecast might be inaccurate to predict the posterior. However, the parameter errors from the solid purple contours do not differ much from the measured errors, so the forecasts are still useful. The main information that we learn from the forecast is that sky coverage and survey depth significantly affects the SFR constraints from CIB tomography. The ongoing CFIS survey is promising to yield similar constraining power as CIB × KiDS + ρ_{SFR}, and LSST can achieve much tighter constraints. In this forecast study, we only consider the L19 CIB map with a relatively low angular resolution. LSST overlaps with Herschel/SPIRE (Viero et al. 2013), which detects IR galaxies with a higher resolution at a level of ∼ 10 arcsec. This makes it possible to study galaxy clustering on smaller scales. Nextgeneration CMB surveys (CMBS4, for example) will also achieve higher angular resolution. In addition, Roman Space Telescope and James Webb Space Telescope will thoroughly explore the nearIR sky. Euclid will also probe IRselected galaxies at high redshifts, which should be useful for CIB studies. When combined with other probes, these new data will provide us with better insight into the cosmic star formation history and the complicated physics behind it.
All Tables
Summary of the free parameters, their prior ranges, and the equations that define them.
Summary of the prior ranges, the marginalised mean values, and the 68% credible intervals of SFR parameters.
Summary of the prior range, the mean, and 68% confidence level of the HOD parameters.
All Figures
Fig. 1. Halo model of the power spectrum of CIBgalaxy crosscorrelation at z = 0. The power spectrum in this plot only shows the spatial dependence of the correlation between the CIB and galaxy distribution, with all the irrelevant terms (redshift and frequency dependence) factored out, so the unit is arbitrary. The dashdotted purple line and the solid red line are one and twohalo terms, respectively; the dashed black line is the summation of one and twohalo terms, and the solid black line is the smoothed power spectrum defined in Eq. (4). 

In the text 
Fig. 2. Redshift distributions of the KiDS1000 gold galaxy sample (solid and dashed lines) and CIB emissions (dotted lines). Solid lines are the redshift distribution of the KiDS galaxies calibrated by the SOM without specifying lensing weight and are the redshift distributions we use in this work. For comparison, we show in dashed lines the redshift distributions from the SOM calibration with lensing weight, which are used in the standard KiDS1000 cosmological analyses. Coloured bands show the z_{B} ranges of corresponding tomographic bins. Dotted lines are dI_{ν}/dz at 353, 545, and 857 GHz calculated from Eq. (22) with the bestfit parameters in this work. The values follow the right y axis. 

In the text 
Fig. 3. CIB map at 545 GHz for the Galactic north pole (left) and south pole (right). The transparent regions are the KiDS footprint. 

In the text 
Fig. 4. Correlation coefficient matrix of our crosscorrelation measurements. The colour scale is from 0 (black) to 1 (white). Each block enclosed by a white grid is the covariance between each pair of crosscorrelations indicated with ticks (bin p × ν GHz), while that enclosed by a golden grid corresponds to the covariance between the CIB crossgalaxies from each pair of tomographic bins. The matrix has nonzero elements at all cells, but the offdiagonal elements in each crosscorrelation are vanishingly small. 

In the text 
Fig. 5. CIBgalaxy crosscorrelations with the five KiDS tomographic bins (rows) and the three CIB maps (columns). The grey points are measured from data, with standard deviation error bars calculated using the square root of the diagonal terms of the covariance matrix. The solid cyan lines show the bestfit crosscorrelation signals calculated using the CIBgalaxy crosscorrelation measurements alone, while the dashed blue lines show the bestfit crosscorrelations when jointly fitting the CIBgalaxy crosscorrelation measurements and the external SFRD. 

In the text 
Fig. 6. Posterior of the SFR parameters. Contours show the 2D posteriors marginalised over all the other 25−2 = 23 parameters. The cyan contours show the constraints from the CIBKiDS crosscorrelation only, the dark blue contours show the constraints from a combination of the crosscorrelation and the external SFRD data, and the red contours show the constraints from the external SFRD data only. The dark and light regions in each contour show the 68% and 95% credible regions, respectively. The dashed lines show the bestfit model from M21. 

In the text 
Fig. 7. Comparison of our constraints on M_{peak} with a number of recent results from the literature. The three points above the dotted line are the results from this work. The other points are colourcoded according to their methods: the green point shows the result from SMG autocorrelations (Chen et al. 2016), the magenta shows the measurements using LRGCIB crosscorrelations (Serra et al. 2014), and the black points show measurements using CIB power spectra (Shang et al. 2012; Viero et al. 2013; Planck Collaboration XXX 2014; Maniyar et al. 2018, 2021). The dark blue band shows the 68% credible interval of our CIB × KiDS + ρ_{SFR} marginal posterior constraint. 

In the text 
Fig. 8. CIB intensity at 545 GHz as a function of redshift while varying M_{peak} and keeping all other parameters fixed (solid lines, right x axis). We also plot the redshift distributions of the galaxy sample Φ_{g}(z) (dashed lines, left y axis). The CIB emissions are calculated from Eq. (22). The black line shows dI_{545}/dz, which corresponds to the CIB × KiDS bestfit parameters. 

In the text 
Fig. 9. Evolution of the SFRD with respect to lookback time (panel a) and redshift (panel b). The SFRD calculated from this work is presented as cyan (crosscorrelation only) and dark blue (crosscorrelation plus external SFRD) lines and shaded regions. The shaded regions enclose the 1σ credible region of the fits and are calculated from 10 000 realisations of SFR parameters from the posterior distribution. The 3σ credible region of the crosscorrelationonly SFRD is also shown in panel b with dotted cyan lines. We note that the lower 3σ limit crosses zero at z ∼ 1.8. The magenta and green lines are the bestfit SFRD from two previous studies, and the points with error bars are the SFRD from previous studies (which are included in our CIB × KiDS + ρ_{SFR} analysis). 

In the text 
Fig. 10. Evolution of the SFRD with respect to lookback time from this work (see Fig. 9), compared to simulations and SAMs. 

In the text 
Fig. 11. Posterior of SFR parameters from the CIB × KiDS fit with fixed (blue) and freely varying (red) priors on the means of Φ_{g}(z). The blue contour is our fiducial posterior. 

In the text 
Fig. 12. Linear galaxy bias constrained from CIBgalaxy crosscorrelation. The solid blue line shows the bestfit , and the band with dotteddashed boundary shows the upper and lower 1σ errors. We also overplot results from previous studies. 

In the text 
Fig. A.1. Comparison between the standard deviation of obtained from the diagonal terms of the analytical (the beige lines) and jackknife (the dashed dark blue lines) covariance matrices. Each panel corresponds to one crosscorrelation signal. The three columns represent three CIB channels in increasing order from left to right; the five rows represent five KiDS tomographic bins in increasing order from top to bottom. 

In the text 
Fig. B.1. Fisher forecast for CIBgalaxy crosscorrelation with LSST (green) and CFIS (dark yellow), as well as the ‘forecast’ for the CIB × KiDS + ρ_{SFR} measurement in this work (solid purple contour). The CIB × KiDS + ρ_{SFR} constraint is also shown as a dashed purple contour. Contours are the 1 and 2σ levels. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.