Cosmic star formation history with tomographic cosmic infrared background-galaxy cross-correlation

In this work, we probe the star formation history of the Universe using tomographic cross-correlation between the cosmic infrared background (CIB) and galaxy samples. The galaxy samples are from the Kilo-Degree Survey (KiDS), while the CIB maps are made from \planck\, sky maps. We measure the cross-correlation in harmonic space with a significance of 43$\sigma$. We model the cross-correlation with a halo model, which links CIB anisotropies to star formation rates (SFR) and galaxy abundance. We assume that SFR has a lognormal dependence on halo mass, while galaxy abundance follows the halo occupation distribution (HOD) model. The cross-correlations give a best-fit maximum star formation efficiency of $\eta_{\mathrm{max}}= 0.41^{+0.09}_{-0.14}$ at a halo mass $\log_{10}(M_{\mathrm{peak}}/M_{\odot})= {12.14\pm 0.36}$. The derived star formation rate density (SFRD) is well constrained up to $z\sim 1.5$. The constraining power at high redshift is mainly limited by the KiDS survey depth. A combination with external SFRD measurements from previous studies gives $\log_{10}(M_{\mathrm{peak}}/M_{\odot})=12.42^{+0.35}_{-0.19}$. This tightens the SFRD constraint up to $z=4$, yielding a peak SFRD of $0.09_{-0.004}^{+0.003}\,M_{\odot} \mathrm { year }^{-1} \mathrm{Mpc}^{-3}$ at $z=1.74^{+0.06}_{-0.02}$, corresponding to a lookback time of $10.05^{+0.12}_{-0.03}$ Gyr. Both constraints are consistent, and the derived SFRD agrees with previous studies and simulations. Additionally, we estimate the galaxy bias $b$ of KiDS galaxies from the constrained HOD parameters and yield 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$. 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.


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 multi-wavelength observations (Gruppioni et al. 2013;Magnelli et al. 2013;Davies et al. 2016;2021), the integrated Sachs-Wolfe (ISW) effect (Maniyar et al. 2019;Hang et al. 2020), and more. Cosmic infrared background maps have been extensively used to study dusty star-forming galaxies via auto-correlations (Shang et al. 2012;Planck Collaboration 2014;Maniyar et al. 2018) and cross-correlations with other large-scale structure tracers (Cao et al. 2020;Maniyar et al. 2021). Clustering-based CIB cross-correlation 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 sub-millimetre 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 clustering-based measurement that allows us to study the cosmic star formation history with the CIB: tomographic cross-correlation between the CIB and galaxy number density fluctuations. That is, cross-correlating the CIB with galaxy samples in different redshift ranges (so-called tomographic bins) to measure the evolution of the CIB over cosmic time. Compared with other large-scale 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 on-sky areas. Alternatively, wide-area 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. 2018 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 non-linear 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 large-scale structure, including re-ionisation (Lidz et al. 2008), cosmological parameters Kuntz (2015), and the ISW effect (Hang et al. 2021). In the near future, the Canada-France 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, CIB-galaxy cross-correlations are only sensitive to the CIB at low redshift, but this will improve with future galaxy surveys. In this study, we cross-correlate the galaxy catalogues provided by the Kilo-degree 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 (2014) and Shang et al. (2012) use a halo model by assuming a lognormal luminosity-to-halo mass (L − M) relation for the IR and a grey-body 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. Sect. 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 (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).

Angular cross-correlation
Both the galaxy and the CIB angular distributions are formalised as the line-of-sight projection of their 3D distributions. This subsection introduces the general theoretical framework of the angular cross-correlation between two projected cosmological fields. For an arbitrary cosmological field u, the projection of its 3D fluctuations (i.e. anisotropies) is written as where ∆ u (θ) is the 2D projection in the angular direction θ, and δ u (χθ, χ) is the fluctuation of u in 3D space at the coordinate (χθ, χ), where χ is the comoving distance. The kernel W u (χ) describes the line-of-sight distribution of the field 2 We measure the angular cross-correlation in harmonic space. In general, the angular cross-correlation between two projected fields, u and v, at the scales 10 are well estimated by the Limber approximation (Limber 1953;Kaiser 1992): Fig. 1. Halo model of the power spectrum of CIB-galaxy 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 dash-dotted purple line and the solid red line are one-and twohalo terms, respectively; the dashed black line is the summation of one-and two-halo terms, and the solid black line is the smoothed power spectrum defined in Eq. (4).
where P uv (k, z) is the 3D cross-power spectrum of associated 3D fluctuating fields u and v: Generally, we can model a large-scale cosmological field as a biased-tracer 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 one-halo term, which accounts for correlations within the same halo. Smith et al. (2003) points out that simply adding the one-and two-halo up yields a total power spectrum that is lower than that from cosmological simulations in the transition regime (k ∼ 0.5 hMpc −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 two-halo terms (dash-dotted purple and red solid lines, respectively) of the CIB-galaxy cross-correlation 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.
Both one-and two-halo 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 (k|M) is the profile of the tracer u with mass M in Fourier space: where p u (r|M) 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) and , respectively, in accordance with M21.

Galaxy number density fluctuations
The 2D projected galaxy number density fluctuation is measured as where n(θ) is the surface number density of galaxies in the directionθ on sky, andn is the average surface number density. Given the redshift distribution of a galaxy sample Φ g (z) (determined by the true line-of-sight galaxy distribution and any survey selection functions) 4 , the projected galaxy density fluctuation is given by where δ g (χ(z)θ, χ) 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 (r|M), 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 (r|M) is the number density profile of the satellite galaxies. Its Fourier transform will be given in Eq. (14).n g (z) is the mean galaxy number density at redshift z, which is given bȳ 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, σ ln M 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 (k|M) is generally taken as the Navarro-Frenk-White (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 concentration-mass 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 1998). 6 The HOD parameters in Eq. (12) depend on redshift (Coupon et al. 2012). In this work, we fix σ ln M = 0.4 5 The cosine and sine integrals are defined as follows: 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): log 10 M 0 (a) = log 10 M 0,0 + log 10 M 0,p (a − 1), 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: In practice, we find that the resolution of the CIB map is sufficiently low that this simple formalism fits the data well (Sect. 5).

Halo model for CIB-galaxy cross-correlation
The intensity of the CIB (in Jy/sr) is the line-of-sight 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 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 ν (k|M) is the profile term p u (k|M) in Eq. (5) for CIB anisotropies. Following the standard practice of van den Bosch et al. (2013), we give the cross-correlation 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.

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 rest-frame frequency ν and redshift z. The latter term is defined as the mean flux density received from a source with L IR = 1L , 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 Jr (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 L −1 assuming a Chabrier initial mass function (IMF; Chabrier 2003). The derivation of this formula can be found in Appendix B of Planck Collaboration (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 sub-halos. Following Shang et al. (2012); Maniyar et al. (2018Maniyar et al. ( , 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 grey-body spectrum, while Maniyar et al. (2018Maniyar et al. ( , 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); Bethermin et al. (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 baryon-matter 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). 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 sub-halos in which they are located (Béthermin et al. 2013): where m is the sub-halo mass, and SFR is the general SFR defined by Eq. (25). The mean SFR for sub-halos in a host halo with mass M is then where dN sub /d ln m is the sub-halo mass function. We take the formulation given by Tinker & Wetzel (2010). Once we get the SFR for both the central and the sub-halos, 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 non-dust radiation, including CO(3-2), free-free 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 shot-noise resulting from the autocorrelated Poisson sampling noise. Therefore, the model for the total CIB-galaxy cross-correlation C νg,tot is where C νg,hm is the cross-correlation predicted by the halo model, and S νg is the scale-independent shot noise. Shotnoise is generally not negligible in galaxy cross-correlations, especially on small scales. There are analytical models to predict shot noise (Planck Collaboration 2014; Wolz et al. 2018) but they depends on various assumptions, including the CIB flux cut, galaxy colours, galaxy physical evolution , 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 cross-correlations, 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.
With the SFR and SED models introduced above, the redshift distribution of the CIB intensity dI ν dz 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 CIB-galaxy cross-correlation.

KiDS data
We used the lensing catalogue provided by the fourth data release (DR4) of KiDS ) as our galaxy sample. KiDS is a wide-field 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 weak-lensing applications. The footprint of the KiDS DR4 (and its corresponding galaxy sample, called the 'KiDS-1000' 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. High-quality optical images are produced with VST-OmegaCAM, and these data are then combined with imaging from the VISTA Kilo-degree INfrared Galaxy (VIKING) survey (Edge et al. 2013), allowing all sources in KiDS-1000 to be photometrically measured in nine optical and near-IR bands: ugriZY JHK s (Wright et al. 2019).
The KiDS-1000 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 near-IR information from VIKING (which was co-designed with KiDS to reach complimentary depths in the near-IR bands). Each galaxy in the KiDS-1000 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 KiDS-1000 sample is then further down-selected during the production of high-accuracy calibrated redshift distributions Hildebrandt et al. 2021) to produce the KiDS- 1000 '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 large-scale structure into the intermediate-redshift Universe.
Although the KiDS survey provides high-accuracy 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 CIB-galaxy correlation signals, which reveal the star formation activity in these galaxies. Therefore, we focus on CIBgalaxy cross-correlation in this work, allowing us to ignore shape information. However, we note that  has made a significant detection of shear-CIB 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 cross-correlation measurement by dividing the galaxy catalogue into five bins, split according to the best-fit photometric redshift estimate z B of each galaxy. These are the same tomographic bins used in the KiDS-1000 cosmology analyses Heymans et al. 2021;Tröster, Tilman et al. 2021). The redshift distribution of each bin is calibrated using a selforganising map (SOM) as described in ; Hildebrandt et al. (2021). A SOM is a 2D representation of A&A proofs: manuscript no. main  an n-dimensional manifold, computed using unsupervised machine learning. For redshift calibration, the SOM classifies the distribution of galaxies in multi-dimensional 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 SOM-calibrated 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 SOM-calibrated Φ 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 higher-redshift bins. We also show the mean CIB emissions (dotted lines) at 353, 545, and 857 GHz calculated from Eq. (22) with the best-fit 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 i-th pixel is calculated as where i denotes the pixel index, n t,i is the surface number density of galaxies in the i-th pixel andn KiDS,t is the mean galaxy surface number density in the t-th 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 ∆ g (θ) with the window function corresponding to the healpix pixel. The mask of the galaxy maps for the cross-correlation measurement is the KiDS footprint, which is presented as the transparent regions in Fig. 3.

CIB data
In this work, we use the large-scale 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 2011. 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 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 all-sky CIB maps (Planck Collaboration 2016) in the three highest HFI frequencies. To make the all-sky 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.

External SFRD data
In addition to the CIB-galaxy cross-correlation 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 cosmology-independent observed intensity between 8-1000 µm per redshift bin, according to corresponding cosmologies, and then convert back to ρ SFR with the cosmology assumed in this study.

Cross-correlation measurements
The cross-correlation between two sky maps, which are smoothed with the beam window function b beam ( ), is related to the real C aŝ whereĈ uv 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 σ = FWHM/ √ 8 ln 2. For the KiDS galaxy overdensity map, the angular resolution is much better than the pixel size, so we assume a FWHM = 0 and b g beam ( ) = 1.
Both the galaxy and the CIB maps are partial-sky. The sky masks mix up modes corresponding to different . The mode-mixing is characterised by the mode-mixing matrix, which depends only on the sky mask. We use the NaMaster ) 11 package to correct mode-mixing and estimate the angular cross-power spectra. NaMaster first naively calculates the cross-correlation between two masked maps. This estimation gives the biased power spectrum, which is called the 'pseudo-C '. Then it calculates the mode-mixing matrices with the input masks and uses this to correct the pseudo-C . 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 sky-coverage of KiDS.
The measurements give a data vector of crosscorrelation between five tomographic bins of KiDS and 3 CIB bands, resulting in 15 cross-correlations C νg : g ∈ {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 square-root of the χ 2 values of the null line (C νg = 0) 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 , τ}. The prior boundaries are given in Table. 3. (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 }. The prior boundaries are given in Table. 4. (iii) Amplitudes of the shot noise power spectra: log 10 S νg . Here S νg is in the unit MJy/sr. 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 degree-of-freedom is 125. Table 2. Summary of the free parameters, their prior ranges, and the equations that define them. The three blocks of the table correspond to three types of parameters: SFR parameters, HOD parameters, and amplitudes of the shot noise power spectra. Note that the last block actually contains 15 shot noise amplitudes.

Covariance matrix
To estimate the uncertainty of the cross-correlation measurement, we followed the general construction of the Fig. 4. Correlation coefficient matrix of our cross-correlation 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 cross-correlations indicated with ticks (bin p × ν GHz), while that enclosed by a golden grid corresponds to the covariance between the CIB cross-galaxies from each pair of tomographic bins. The matrix has non-zero elements at all cells, but the off-diagonal elements in each cross-correlation are vanishingly small.
analytical covariance matrix in the literature. Compared with simulation or resampling-based methods, an analytical method is free of sampling noise and allows us to separate different contributions. Following , we decompose the crosscorrelation covariance matrix into three parts: Here Cov is the abbreviation of Cov uv,wz 1 2 ≡ Cov C uv 1 , C wz 2 . 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 all-sky measurement. Sky masks introduce non-zero coupling between different as well as enlarge the variance. To account for this, we used the method given by Efstathiou (2004) and García-García et al. (2019) that is implemented in the NaMaster package . 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 four-halo terms. Schaan et al. (2018) shows that the one-halo term dominates the CIB trispectrum. As galaxies have a similar spatial distribution to the CIB, we only take the one-halo term into account for our CIB-galaxy cross-correlation: 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 cross-correlation. The SSC can also be calculated in the halo model framework (Lacasa et al. 2018;Osato & Takada 2021).
In this work, the non-Gaussian 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 CIB-galaxy cross-correlation 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 cross-correlations, giving a matrix with a side-length 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 best-fit model parameters. We then updated the covariance matrix with these best-fit 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 off-diagonal terms, which means that the cross-correlations 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 off-diagonal golden blocks is weaker but still non-negligible: 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 , who study tSZ cross-correlations and find that the SSC term was a more significant contributor to their covariance (contributing ∼ 20% to their off-diagonal covariance terms). The reason for this difference is that, compared to the tSZ effect, the galaxy distribution is more concentrated. This causes the non-Gaussian term to remain insignificant until considerably smaller scales than the tSZ effect: beyond the scales probed in this study ( > 2000).
Finally, an alternative estimation of the covariance matrix is shown in Appendix A.

CIB colour-correction and calibration error
The flux given in the Planck sky maps follows the photometric convention that νI ν =constant (Planck Collaboration 2014). The flux therefore has to be colour-corrected for sources with a different SED. Therefore, the CIB-galaxy cross-correlation should also be corrected as where cc ν is the colour-correction factor at frequency ν. In this work, we adopt the colour-correction factors given by  with values of 1.097, 1.068, and 0.995 for the 353, 545, and 857 GHz bands, respectively. Additionally, in Maniyar et al. (2018Maniyar et al. ( , 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.

Cosmic magnification
The measured galaxy overdensity depends not only on the real galaxy distribution but also on lensing magnification induced by the line-of-sight mass distribution (Schneider 1989;Narayan 1989). This so-called cosmic magnification has two effects on the measured galaxy overdensity: i) overdensities along the line-of-sight cause the local angular separation between source galaxies to increase, so the galaxy spatial distributions are diluted and the cross-correlation is suppressed; and ii) lenses along the line-of-sight magnify the flux of source galaxies such that fainter galaxies enter the observed sample, so the overdensity typically increases. These effects potentially bias galaxy-related crosscorrelations, especially for high-redshift 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 best-fit CIB parameters.

Redshift distribution uncertainty
The SOM-calibrated galaxy redshift distributions have an uncertainty on their mean at a level of ∼ 0.01 Hildebrandt et al. (2021), which could affect galaxy cross-correlations. 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 bỹ whereΦ g,i (z) is the shifted galaxy redshift distribution in the i th tomographic bin, and δ z,i is the shift of the redshift distribution parameter in the i th bin. The priors on δ z,i are assumed to be covariant Gaussians centred atδ z,i (i.e. the mean δ z,i ). Hildebrandt et al. (2021) evaluated bothδ z,i 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δ z,i nor their covariance matrix from Hildebrandt et al. (2021) are formally correct for the Φ g (z) in this work. To correctly estimateδ z,i 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 lensing-weighted 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 KiDS-1000 δ z,i covariance matrix. This yields an absolute uncertainty at a level of 0.04, about two times the difference between the nominal KiDS-1000 lensing-weighted Φ g (z) and the unweighted Φ g (z) that we use in this work.

Likelihood
We constrained the model parameters in two ways. The first is cross-correlation only ('CIB×KiDS fitting' hereafter), and the second is combining cross-correlation and external ρ SFR measurements ('CIB × KiDS + ρ SFR fitting' hereafter). For the cross-correlation only fitting, since we are working with a wide range, there are many degrees-of-freedom in each bin. According to the central limit theorem, the bin-averaged 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 vectorC is a concatenation of our 15 measured CIB-galaxy cross-correlations; 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 cross-correlation and similarly independent at different redshifts. Therefore, including them introduces an additional term in the likelihood: whereρ SFR (z i ) 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 cross-correlation-only 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 multi-wavelength 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 Markov-Chain Monte-Carlo method to constrain our model parameters with the python package emcee (Foreman-Mackey et al. 2013). Best fit parameters are determined from the resulting chains, as being the sample with the smallest χ 2 goodness-of-fit. Marginal constraints on parameters, when quoted, are marginal means and standard deviations.

Constraints on star formation rate density
We show our CIB-KiDS cross-correlation measurement in Fig. 5. Each panel presents the cross-correlation 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 degrees-of-freedom 125 and 142, respectively. In order to evaluate the goodness-of-fit, we calculate the corresponding probability-to-exceed (PTE) given the degree-of-freedom: 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 cross-correlations well. We also notice that the fitting in low-redshift bins is worse than high-redshift bins (see the '0.1 < z B ≤ 0.3, 353 GHz' panel in Fig. 5, for example), although correlation in the points makes 'chi-by-eye' inadvisable here.
We estimate the posterior with the last 100000 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 cross-correlation only results are consistent with the results when analysing external SFRD data only (i.e. the red contours). This validates our CIB-galaxy 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 cross-correlation constraints are generally tighter than those made with the external SFRD data alone, demonstrating that the cross-correlation 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.
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 colour-coded according to their methods: the green point shows the result of SMG auto-correlations (Chen et al. 2016), the magenta point shows the measurement using LRG-CIB cross-correlations (Serra et al. 2014), and black points show measurements utilising CIB power spectra (Shang et al. 2012;Viero et al. 2013;Planck Collaboration 2014;Maniyar et al. 2018Maniyar et al. , 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. 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 line-of-sight 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 × Fig. 5. CIB-galaxy cross-correlations 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 best-fit cross-correlation signals calculated using the CIB-galaxy cross-correlation measurements alone, while the dashed blue lines show the best-fit cross-correlations when jointly fitting the CIB-galaxy crosscorrelation measurements and the external SFRD.
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 cross-correlation-only 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 . 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 10000 realisations of our model parameters, drawn from the poste- 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 CIB-KiDS cross-correlation only, the dark blue contours show the constraints from a combination of the cross-correlation 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 best-fit model from M21. Table 3. Summary of the prior ranges, the marginalised mean values, and the 68% credible intervals of SFR parameters. The values and errors are calculated from the posteriors marginalised over all the other parameters. The last column shows the constraint from M21 as a reference. The last two rows show the χ 2 values, degrees-of-freedom, and PTE for our constraints.  Fig 6. The 68% credible interval is shown as bands with corresponding colours, which are calculated by computing the 10000 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 best-fit 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 cross-correlation power spec- 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 colour-coded according to their methods: the green point shows the result from SMG auto-correlations (Chen et al. 2016), the magenta shows the measurements using LRG-CIB 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 2014;Maniyar et al. 2018Maniyar et al. , 2021. The dark blue band shows the 68% credible interval of our CIB × KiDS + ρ SFR marginal posterior constraint. tra 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 cross-correlations 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 0.08 +0.03 −0.04 M year −1 Mpc −3 at z = 1.94 +0.1 −0.51 , corresponding to a lookback time of 10.42 +0.16 −1.07 Gyr, while CIB × KiDS + ρ SFR constraint gives a peak SFRD of 0.09 +0.003 −0.004 M year −1 Mpc −3 at z = 1.74 +0.06 −0.02 , corresponding to a lookback time of 10.05 +0.12 −0.03 Gyr, consistent with previous observations of the cosmic star formation history.

Parameter Prior
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 Gonzalez-Perez et al. 2014purple line), EAGLE (Guo et al. 2016, with the Schaye et al. 2015khaki line), and the L-GALAXIES 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 gas-cooling models. We see that EAGLE predicts a slightly lower SFRD (at essentially all times) than is predicted from our results. galform, on the other-hand, agrees well with our CIB×KiDS fits at intermediate-to-early 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 Gonzalez-Perez et al. (2014) version of galform. L-GALAXIES 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 late-time 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 CIB-galaxy cross-correlation 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.
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 redshift-distribution self-calibration 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).
A&A proofs: manuscript no. main (a) (b) 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 (cross-correlation only) and dark blue (cross-correlation plus external SFRD) lines and shaded regions. The shaded regions enclose the 1σ credible region of the fits and are calculated from 10000 realisations of SFR parameters from the posterior distribution. The 3σ credible region of the cross-correlation-only 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 best-fit 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). Fig. 10. Evolution of the SFRD with respect to lookback time from this work (see Fig. 9), compared to simulations and SAMs.

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 full-depth 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 galaxy-CMB lensing cross-correlations, assuming a linear model. In this work, we derive the linear galaxy bias from the constrained HOD parameters. The scale-dependent 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 Fig. 12. Linear galaxy bias constrained from CIB-galaxy crosscorrelation. The solid blue line shows the best-fit b lin g (z) , and the band with dotted-dashed boundary shows the upper and lower 1σ errors. We also over-plot results from previous studies. The constrained HOD parameters are presented in Table 4. Similar to the calculation of the best-fit SFRD and its uncertainty, we calculate the best-fit and 1σ credible interval of b lin g (z) and present it in Fig 12. The resulting linear galaxy bias increases from 1.1 +0.17 −0.31 at z = 0 to 1.96 +0.18 −0.64 at z = 1.5. We also over-plot constraints from previous studies on galaxy bias of star-forming galaxies. The magenta line shows the best-fit model from Maniyar et al. (2018); the green line shows the best-fit 'star-forming' model from Durkalec et al. (2015); the red points are the derived galaxy bias of star-forming 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.
The galaxy bias constrained from CIB-galaxy crosscorrelation is slightly higher than that constrained from galaxy-CMB 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 CIB-galaxy crosscorrelation in this work is mainly sensitive to KiDS galaxies that are active in star formation. The fact that CIB cross-correlation gives higher galaxy bias means that these star-forming galaxies are (on average) more clustered than the galaxies detected by optical survey, especially at high redshift, which calls for further study.

Conclusions
In this work we measure the tomographic cross-correlation between the KiDS-1000 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 cross-correlation 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 KiDS-CIB cross-correlation. Given the degreesof-freedom, 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 cross-correlation 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 multi-wavelength studies of galaxies. -With our cross-correlation measurement, the maximum star formation efficiency is η max = 0.41 +0.09 −0.14 , in agreement with M21. Our CIB × KiDS -only 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 cross-correlation and external SFRD measurements (the CIB×KiDS+ρ SFR constraints), however, tightens the constraint on log 10 M peak to log 10 M peak = 12.42 +0.35 −0.19 , 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 best-fit 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), cross-correlations give an SFRD that is consistent with previous galaxy-based multi-wavelength 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 signal-to-noise, this may become important.

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 cross-correlations between CIB and large-scale structure surveys. We make several assumptions to simplify the CIB halo model. For example, we assume that the SFR-IR luminosity connection can be described by a single Kennicutt parameter in the IR bands (Kennicutt Jr 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 stellar-mass 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 mid-IR signal at z > 2 ( Bethermin 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 cross-correlation 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 year-3 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 LSST-CIB 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 cross-correlating 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 In addition, dust properties should also depend on these properties (Wolf et al. 2009;Noll & Pierini 2005). A CIB cross-correlation 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 cross-correlation 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 star-forming 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 b lin g (z) . 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 CIB-galaxy cross-correlation 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 CIB-galaxy crosscorrelation. Larger, deeper datasets, coupled with more complex sample subdivisions, will allow us to leverage CIB cross-correlations to better understand the growth and evolution of galaxies in our Universe. . 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 cross-correlation data and the MCMC chain of this work will be shared on reasonable request to the corresponding author. Fisher forecast for CIB-galaxy cross-correlation 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.