# The star formation history in the last 10 billion years from CIB cross-
correlations

Baptiste Jego1,2, Jaime Ruiz-Zapatero2, Carlos García-García2, Nick
Koukoufilippas2, and David Alonso2
1ENS Paris-Saclay, Gif-sur-Yvette, France
2Department of Physics, University of Oxford, Denys Wilkinson Building, Keble
Road, Oxford OX1 3RH, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Cosmic Infrared Background (CIB) traces the emission of star-forming
galaxies throughout all cosmic epochs. Breaking down the contribution from
galaxies at different redshifts to the observed CIB maps would allow us to
probe the history of star formation. In this paper, we cross-correlate maps of
the CIB with galaxy samples covering the range $z\lesssim 2$ to measure the
bias-weighted star-formation rate (SFR) density $\langle b\rho_{\rm
SFR}\rangle$ as a function of time in a model independent way. This quantity
is complementary to direct measurements of the SFR density $\rho_{\rm SFR}$,
giving a higher weight to more massive haloes, and thus provides additional
information to constrain the physical properties of star formation. Using
cross-correlations of the CIB with galaxies from the DESI Legacy Survey and
the extended Baryon Oscillation Spectroscopic Survey, we obtain high signal-
to-noise ratio measurements of $\langle b\rho_{\rm SFR}\rangle$, which we then
use to place constraints on halo-based models of the star-formation history.
We fit halo-based SFR models to our data and compare the recovered $\rho_{\rm
SFR}$ with direct measurements of this quantity. We find a qualitatively good
agreement between both independent datasets, although the details depend on
the specific halo model assumed. This constitutes a useful robustness test for
the physical interpretation of the CIB, and reinforces the role of CIB maps as
valuable astrophysical probes of the large-scale structure. We report our
measurements of $\langle b\rho_{\rm SFR}\rangle$ as well as a thorough account
of their statistical uncertainties, which can be used to constrain star
formation models in combination with other data.

###### keywords:

cosmology: large-scale structure of the Universe – galaxies: star formation

††pubyear: 2022††pagerange: The star formation history in the last 10 billion
years from CIB cross-correlations–5

##  1 Introduction

Dust in star-forming galaxies absorbs the stellar ultraviolet radiation from
newly-formed, short-lived, massive stars, obscuring it, and re-emits it in the
infrared (IR). The combined diffuse IR emission from all such sources thus
constitutes a tracer of star formation (Partridge & Peebles, 1967; Knox et
al., 2001), and is known as the Cosmic Infrared Background (CIB).

Since its first detection by COBE (Dwek et al., 1998) (although see Puget et
al. (1996)), the CIB has been both a nuisance and a blessing for cosmology and
astrophysics. On the one hand, it is the most significant extra-galactic
foreground affecting observations of the Cosmic Microwave Background (CMB)
temperature fluctuations at high frequencies (Planck Collaboration et al.,
2011; Dunkley et al., 2013). Unlike other foregrounds, the fact that the IR
spectral energy distribution (SED) is not universal across different galaxies
and cosmic times, makes the CIB particularly difficult to separate in multi-
frequency CMB observations. This complicates not only the cosmological
analysis of small-scale CMB data, but also the recovery of unbiased maps of
the CIB itself in the presence of Galactic dust contamination (Planck
Collaboration et al., 2014b, 2016; Lenz et al., 2019). Furthermore, since the
CIB fluctuations trace the non-Gaussian large-scale structures that develop at
low redshifts, it is an important systematic for CMB lensing analyses, which
exploit the non-Gaussian modifications caused by the gravitational lensing of
the otherwise Gaussian CMB fluctuations (Osborne et al., 2014; van Engelen et
al., 2014; Sailer et al., 2021; Darwish et al., 2021a). Given the central role
that small-scale primary CMB data and CMB lensing will play in cosmology in
the next decade (Schmittfull & Seljak, 2018; Fang et al., 2022), improving our
understanding of the CIB, and developing better models that allow us to
quantify its impact on these observations, is of paramount importance.

On the other hand, the CIB is in itself a treasure trove for both astrophysics
and cosmology. First, it is a high-sensitivity tracer of the star-formation
rate (SFR) density (i.e. the total mass of stars created per unit time and
volume) (Dole et al., 2006). Understanding star formation across time and
galaxy type is a central piece in the puzzle of galaxy formation and evolution
(Tinsley, 1980). Through observations of the infrared luminosity function at
high redshifts (Gruppioni et al., 2013; Magnelli et al., 2013; Marchetti et
al., 2016; Davies et al., 2016a), a qualitative picture has emerged where star
formation grows from the epoch of reionisation, peaks at around $z\sim 2$, and
then decreases as the gas that fuels it runs out (Madau & Dickinson, 2014b).
However, collecting large samples of star-forming galaxies at high redshift is
observationally challenging. Thus, since the CIB combines the emission of all
IR sources since $z\sim 6$, it encodes information about the full history of
star formation, making it a useful complementary probe for these studies.

Unfortunately, studies of the star formation history from maps of the CIB
alone (e.g. Shang et al., 2012; Viero et al., 2013; Planck Collaboration et
al., 2014b; Maniyar et al., 2018; Maniyar et al., 2021) are hampered by its
projected nature. Without additional information it is not possible to
disentangle, in a model-independent manner, the contribution to the CIB from
sources at different redshifts. Although cross-correlations with CMB lensing
can be used to improve constraints on different CIB models (Planck
Collaboration et al., 2014a; Maniyar et al., 2018; Cao et al., 2020; Darwish
et al., 2021b; McCarthy & Madhavacheril, 2021), they suffer from the same
problem. This, however, can be solved through the use of galaxy clustering
tomography (Chiang et al., 2019). Since the sources that give rise to the CIB
trace the same large-scale structure as any other galaxy sample, the amplitude
of the angular correlation between the CIB and a set of galaxies at a
sufficiently narrow range of redshifts will be proportional to the
contribution to the CIB map from sources at those redshifts. This technique
has been widely utilised in various areas of cosmology and astrophysics to
reconstruct the evolution, or the redshift dependence, of different
quantities. These applications include recovering the mean gas pressure from
maps of the thermal Sunyaev-Zel’dovich effect (Pandey et al., 2019;
Koukoufilippas et al., 2020; Chiang et al., 2020; Yan et al., 2021), the
cosmic UV background (Chiang et al., 2019), and the growth of structure
(García-García et al., 2021; Krolewski et al., 2021), as well as the
calibration of redshift distributions for imaging galaxy samples (the so-
called “clustering redshifts” method Newman, 2008; Ménard et al., 2013;
Morrison et al., 2017). Recently Yan et al. (2022) used cross-correlations
between CIB maps and galaxies from the Kilo-Degree Survey (KiDS) to constrain
a physical model of the CIB (see also Serra et al. (2014); Wang et al. (2015);
Chen et al. (2016)).

The aim of this paper is to use cross-correlations of the CIB with samples of
galaxies at different redshifts to reconstruct the star formation history
since $z\sim 2$. Our approach will be somewhat different from that of Yan et
al. (2022). First, we will show that, from the amplitude of the galaxy auto-
correlation and the cross-correlation with the CIB on large scales, robust and
model-independent constraints can be placed on the evolution of the bias-
weighted star-formation rate density $\langle b\rho_{\rm SFR}\rangle$. In this
quantity, the contribution to the total SFR density from each source is
weighted by the linear clustering bias of the dark matter halo it resides in.
$\langle b\rho_{\rm SFR}\rangle$ is therefore more sensitive to the
contribution from galaxies in massive haloes than direct measurements of the
SFR density $\rho_{\rm SFR}$, and thus provides a complementary probe for the
physics of star formation. The model-independent measurements of $\langle
b\rho_{\rm SFR}\rangle$ obtained can then be used to constrain the parameters
of a halo-based SFR model. This will allow us to propagate our measurements
onto constraints on $\rho_{\rm SFR}$ that can be compared with direct
measurements of this quantity from studies of the infrared luminosity
function, thus providing a useful consistency test of the physical models used
to describe star formation.

This paper is structured as follows. Section 2 lays out the theoretical
framework, based on the halo model, that will allow us to measure $\langle
b\rho_{\rm SFR}\rangle$ as a function of redshift. The galaxy samples and CIB
maps used in this analysis are then presented in Section 3. Section 4 presents
our model-independent measurements of $\langle b\rho_{\rm SFR}\rangle$,
quantifies their robustness, and propagates them onto constraints on
$\rho_{\rm SFR}$ and on two halo-based models for star formation. We then
present our conclusions and outlook in Section 5.

##  2 Methods

###  2.1 Theory

####  2.1.1 Projected probes

Let $u(\hat{\bf n})$ be a cosmological field projected on the sphere, and
assume that it is the radial projection of a 3D field $U$, i.e.

$u(\hat{\bf n})=\int d\chi q_{u}(\chi)\,U(\chi\hat{\bf n},z).$ (1)

Here $\chi$ is the comoving distance, $q_{u}(\chi)$ is the radial kernel
associated with $u$, and $z$ is the redshift corresponding to the comoving
distance $\chi$. The angular power spectrum of two such quantities, $u$ and
$v$, $C_{\ell}^{uv}$, defined as the covariance between their harmonic-space
coefficients, is then related to the power spectrum of the corresponding 3D
quantities, $P_{UV}(k,z)$, via a similar radial projection:

$C_{\ell}^{uv}=\int\frac{d\chi}{\chi^{2}}q_{u}(\chi)q_{v}(\chi)P_{UV}\left(k=\frac{\ell+1/2}{\chi},z\right).$
(2)

Here we have used the Limber approximation (Limber, 1953), which is
appropriate when the radial kernels are much broader than the typical
correlation length of the fields (this is true for the quantities studied
here).

####  2.1.2 Halo model

The halo model (Seljak, 2000; Peacock & Smith, 2000; Cooray & Sheth, 2002) is
a versatile formalism to model the 3D correlations of cosmological quantities.
In this framework, the power spectrum $P_{UV}(k)$ can be written as the sum of
two contributions, the one-halo and two-halo terms, $P^{1h}_{UV}(k)$ and
$P^{2h}_{UV}(k)$, with

$\displaystyle P^{1h}_{UV}(k)\equiv\int dM\,n(M)\langle U(k,M)V(k,M)\rangle.$
(3) $\displaystyle P^{2h}_{UV}(k)\equiv\langle bU(k)\rangle\langle
bV(k)\rangle P_{\rm lin}(k).$ (4)

Here, $n(M)$ is the halo mass function (the comoving number density of haloes
per mass interval), $U(k,M)$ is the profile of quantity $U$ around a halo of
mass $M$ in Fourier space, $P_{\rm lin}(k)$ is the linear matter power
spectrum, and $\langle bU(k)\rangle$ is

$\langle bU(k)\rangle\equiv\int dM\,n(M)\,b_{h}(M)\,\langle U(k,M)\rangle.$
(5)

On large scales, both $P^{1h}_{UV}$ and $\langle bU\rangle$ tend to a
constant. This can be seen easily by writing $U(k,M)$ in the limit
$k\rightarrow 0$ (i.e. when $e^{-i{\bf k}\cdot{\bf x}}\rightarrow 1$):

$U(k\rightarrow 0,M)=\int d^{3}x\,U({\bf x},M).$ (6)

Thus, on sufficiently large scales, the power spectrum can be approximated as

$P_{UV}(k)=\langle bU\rangle\langle bV\rangle P_{\rm lin}(k)+N_{UV},$ (7)

where $\langle bU\rangle$, $\langle bV\rangle$ and $N_{UV}$ are now scale-
independent quantities. Crucially, $\langle bU\rangle$, has a simple physical
interpretation: it is the cosmic mean value of quantity $U$ weighted by the
halo bias.

Extracting information from scales smaller than those for which the
approximation above is appropriate requires tackling several modelling
challenges. First, an accurate model of the scale dependence of the halo
profiles (i.e. the distribution of $U$ and $V$ around haloes of different
masses) is necessary. Secondly, the 1-halo term depends on the statistical
correlation of the two quantities under study, which in principle requires a
model of not only the mean profiles, but also of their statistical
correlation. Although these are often implicitly assumed to be independent
(i.e. $\langle UV\rangle=\langle U\rangle\langle V\rangle$), deviations from
this can be highly relevant and lead to incorrect conclusions (particularly
given the high precision with which small-scale correlations can be measured).
In turn, the 2-halo term requires no knowledge of the joint distribution of
the quantities. Finally, the halo model is known to fail at the $\gtrsim 10\%$
level on scales straddling the transition between the 1-halo- and 2-halo-
dominated regimes. This is likely caused by the overly simplified assumption
of linear biasing used in vanilla versions of the halo model to describe the
clustering of dark matter haloes (Mead et al., 2021; Mead & Verde, 2021).

Tackling these challenges is a complex problem that can be addressed through
more sophisticated astrophysical modelling, often at the cost of introducing
new nuisance parameters that must be marginalised over. The obvious benefit of
this approach is that it enables the use of smaller scales, potentially
achieving tighter constraints on the model parameters. This is indeed a common
approach used in previous analyses of CIB auto- and cross-power spectra (e.g.:
Maniyar et al. (2018); Yan et al. (2022)).

Here, we choose to avoid this additional complexity, and opt for a robust and
model-independent measurement of the phenomenological parameters in Eq. 7,
restricting our analysis to large scales. Once the bias-weighted quantities
$\langle bU\rangle$ have been measured, they can be used to constrain the
underlying physical model describing them (in this case the star-formation
rate history). All that remains is designing this model for the two tracers
under consideration: the projected galaxy overdensity and the CIB emissivity.

####  2.1.3 CIB and star formation

The CIB intensity at a given observed frequency $\nu$ is related to the
infrared emissivity (energy emitted per unit frequency, time and comoving
volume in the rest frame of the emitter) $j_{\nu(1+z)}$ via:

$I_{\nu}(\hat{\bf n})=\int d\chi\frac{j_{\nu(1+z)}(\chi\hat{\bf
n},z)}{4\pi(1+z)}.$ (8)

Given a model for the specific infrared luminosity of haloes of mass $M$,
$L_{\nu_{e}}(M)$, the mean comoving emissivity becomes

$\langle j_{\nu(1+z)}\rangle=\int dMn(M)\,L_{\nu(1+z)}(M),$ (9)

where $n(M)$ is the halo mass function. We now proceed as in Maniyar et al.
(2018); Maniyar et al. (2021). The specific infrared luminosity can be written
in terms of the observed source flux $S_{\nu}$

$L_{\nu(1+z)}=4\pi\chi^{2}(1+z)S_{\nu}.$ (10)

At the same time, the total infrared luminosity $L_{\rm IR}$ (integrated over
frequencies) is tightly correlated with the star-formation rate of the source:

${\rm SFR}=K\,L_{\rm IR},$ (11)

where the $K=10^{-10}\,M_{\odot}\,{\rm yr}^{-1}L_{\odot}^{-1}$ is the
calibration constant between the far infrared luminosity and star formation
rate Kennicutt (1998); Kennicutt & Evans (2012) for a Chabrier initial mass
function (Chabrier, 2003), and $L_{\rm IR}$ is integrated over the range
8-1000$\mu$m. Combining Eqs. 9, 10, and 11, the CIB intensity can be written
as

$I_{\nu}(\hat{\bf n})=\int d\chi\,\frac{\chi^{2}S_{\nu}^{\rm
eff}(z)}{K}\,\rho_{\rm SFR}(\chi\hat{\bf n},z),$ (12)

where $\rho_{\rm SFR}$ is the star-formation rate density, given by the
contribution of all haloes of different masses. $S^{\rm eff}_{\nu}(z)\equiv
S_{\nu}(z)/L_{\rm IR}$ is the mean flux of sources at redshift $z$ normalised
to a total luminosity of $L_{\odot}$. We use the estimates of $S^{\rm
eff}_{\nu}(z)$ for the three Planck frequency channels derived by Béthermin et
al. (2013, 2015, 2017), and made available by Maniyar et al. (2021)111See
https://github.com/abhimaniyar/halomodel_cib_tsz_cibxtsz..

Comparing this with Eq. 1, the CIB intensity $I_{\nu}(\hat{\bf n})$ can
therefore be understood as a projected tracer of the star-formation rate
density $\rho_{\rm SFR}$, with a radial kernel given by

$q_{\nu}(\chi)=\frac{\chi^{2}S_{\nu}^{\rm eff}(z)}{K}.$ (13)

The associated $\langle bU\rangle$ is the mean star-formation rate density
weighted by halo bias:

$\langle b\rho_{\rm SFR}\rangle\equiv\int dM\,n(M)\,b_{h}(M)\,{\rm SFR}(M),$
(14)

where ${\rm SFR}(M)$ is the mean star-formation rate in haloes of mass $M$.
The mean profile, in turn, would provide a measurement of the mean star-
formation rate density

$\rho_{\rm SFR}=\int dM\,n(M)\,{\rm SFR}(M).$ (15)

When interpreting the measured values of $\langle b\rho_{\rm SFR}\rangle$, we
will make use of a halo-model parametrisation similar to that of Moster et al.
(2018) and Maniyar et al. (2021) to describe the SFR-mass relation. The model
splits the total SFR of a given halo into the contribution from central (${\rm
SFR}_{c}$) and satellite galaxies (${\rm SFR}_{s}$):

${\rm SFR}(M,z)={\rm SFR}_{c}(M,z)+{\rm SFR}_{s}(M,z).$ (16)

As haloes accrete gas, part of it is transformed into stars. The model
parametrises this process as:

${\rm SFR}_{c}(M,z)=\eta(M,z)\,{\rm BAR}(M,z),$ (17)

where ${\rm BAR}$ is the baryon accretion rate, and $\eta$ is the efficiency
with which those baryons are transformed into stars. The BAR is modelled after
Fakhouri et al. (2010) as

${\rm
BAR}(M,z)=\dot{M}_{0}\frac{\Omega_{b}}{\Omega_{M}}\,\left(\frac{M}{10^{12}M_{\odot}}\right)^{1.1}(1+1.11z)\frac{H(z)}{H_{0}},$
(18)

with $\dot{M}_{0}=46.1M_{\odot}\,{\rm yr}^{-1}$.

We will consider two parametrisations for the efficiency $\eta(M,z)$:

  1. 1.

Model A: the model of Maniyar et al. (2021). In this case

$\eta(M,z)=\eta_{\rm max}\exp\left[-\frac{(\log M_{\rm max}-\log
M)^{2}}{2\sigma_{M}^{2}(M,z)}\right],$ (19)

where $M_{\rm max}$ is the halo mass of maximum efficiency $\eta_{\rm max}$.
The logarithmic width $\sigma_{M}$, describing the range of masses with
efficient star formation, depends on both mass and redshift as

$\sigma_{M}(M,z)=\sigma_{M,0}-\tau\Theta(M-M_{\rm max})\,{\rm
max}(0,z_{c}-z),$ (20)

where $\Theta$ is the Heaviside step function.

Fixing $z_{c}=1.5$ as in Maniyar et al. (2021), the model is described by 4
free parameters: $\eta_{\rm max}$, $\log_{10}M_{\rm max}$, $\sigma_{M,0}$, and
$\tau$.

  2. 2.

Model B: the parametrisation of Moster et al. (2018). In this case:

$\eta(M,z)=\frac{2\eta_{*}}{(M_{1}/M)^{\beta}+(M/M_{1})^{\gamma}},$ (21)

where $\eta_{*}$, $m\equiv\log_{10}M_{1}/M_{\odot}$, $\beta$, and $\gamma$ are
functions of redshift of the form:

$x(z)=x_{0}+x_{z}\frac{z}{1+z}$ (22)

Fixing $(\beta_{0},\beta_{z},\gamma_{0},\gamma_{z})$ to the best-fit values
found by Moster et al. (2018), $(3.344,-2.079,0.966,0)$, the model is
described by 4 free parameters: $\eta_{0}$, $\eta_{z}$, $m_{0}$, and $m_{z}$.

In order to compare this model with Model A, note that the peak efficiency in
this case is $\eta_{\rm max}=\eta_{*}f_{\eta}$, and it is attained at a mass
$M_{\rm max}=M_{1}\,f_{M}$, where

$f_{M}=\left(\frac{\beta}{\gamma}\right)^{\frac{1}{\beta+\gamma}},\hskip
6.0ptf_{\eta}=\frac{2\gamma}{\beta+\gamma}f_{M}^{\beta}.$ (23)

For the values of $(\beta_{0},\beta_{z},\gamma_{0},\gamma_{z})$ used here,
both factors are roughly constant and close to 1 ($f_{M}$ varies between 1.33
and 1.23, and $f_{\eta}$ varies between 1.17 and 1.04, in the range $0<z<4$).
Thus, to a good approximation, $\eta_{*}$ and $M_{1}$ provide the peak
efficiency and its associated halo mass.

Figure 1: Mass-dependent star-formation efficiency at $z=0.5$ (dotted lines),
$z=1$ (dashed lines) and $z=2$ (solid lines). Results are shown for the two
halo-based models described in the text: Model A, from Maniyar et al. (2021),
in black, and Model B, from Moster et al. (2018), in red. The curves shown
correspond to the best-fit parameters found in the two papers
($\log_{10}M_{\rm max}=12.94$, $\eta_{\rm max}=0.42$, $\sigma_{M,0}=1.75$,
$\tau=1.17$ for Model A, $m_{0}=11.34$, $m_{z}=0.692$, $\eta_{0}=0.005$,
$\eta_{z}=0.689$ for Model B).

Both models are qualitatively similar in that the star formation efficiency
peaks at a particular halo mass, and decreases for lighter and heavier haloes.
This is physically well motivated, based on the effect of low gravitational
potentials and supernova feedback at low masses (Silk, 2003), and an increased
gas cooling time at high masses (Kereš et al., 2005). There are, however,
important differences between them. On the one hand, the peak efficiency in
Model A is constant as a function of time, and all time evolution is
incorporated in the extent of the high-mass tail of $\eta(M,z)$. On the other,
the peak efficiency is explicitly time-dependent in Model B, while the slope
of the high-mass tail is time-independent. Furthermore, while the efficiency
falls exponentially away from the peak mass in Model A, this drop is slower
(power-law) in Model B. These differences are illustrated in Fig. 1, and they
will allow us to explore the dependence of our final constraints on $\rho_{\rm
SFR}$ from our measurements of $\langle b\rho_{\rm SFR}\rangle$ on the
assumptions of the underlying model.

Finally, the contribution from satellite galaxies will be calculated by adding
the contribution of all subhaloes for each parent halo of mass $M$:

${\rm SFR}_{s}(M,z)=\int_{M_{\rm min}}^{M}dM_{\rm sub}\,\frac{dN}{dM_{\rm
sub}}\,{\rm SFR}_{\rm sat}(M_{\rm sub},M,z),$ (24)

where $dN/dM_{\rm sub}$ is the subhalo mass function, calculated as in Tinker
& Wetzel (2010), and ${\rm SFR}_{\rm sat}(M_{\rm sub},M,z)$ is the SFR in a
satellite galaxy with subhalo mass $M_{\rm sub}$ in a parent halo of mass $M$.
As in Maniyar et al. (2021), we model ${\rm SFR}_{\rm sat}$ as

${\rm SFR}_{\rm sat}(M_{\rm sub},M,z)={\rm Min}\left[{\rm SFR}_{c}(M_{\rm
sub},z),\frac{M_{\rm sub}}{M}{\rm SFR}_{c}(M,z)\right].$

####  2.1.4 Projected galaxy clustering

The projected galaxy overdensity is related to its three-dimensional
counterpart via:

$\delta_{g}(\hat{\bf n})=\int d\chi H(z)p(z)\Delta_{g}(\chi\hat{\bf n},z),$
(25)

where $p(z)$ is the galaxy redshift distribution for the galaxy sample under
consideration, and $\Delta_{g}$ is its three-dimensional overdensity.
Therefore, the kernel for this tracer is $q_{g}=H(z)\,p(z)$. The associated
$\langle bU\rangle$ is the galaxy bias $b_{g}$.

A popular halo model parametrisation for galaxy clustering is the halo
occupation distribution (HOD) (Berlind & Weinberg, 2002), based on a model for
the number of galaxies as a function of halo mass $\bar{N}_{g}(M)$. In this
framework, $b_{g}$ would simply be given by

$b_{g}=\frac{\int dM\,n(M)\,b_{h}(M)\,\bar{N}_{g}(M)}{\int
dM\,n(M)\,\bar{N}_{g}(M)}.$ (26)

We will fit directly for the value of $b_{g}$, and thus our treatment is
insensitive to the details of the $\bar{N}_{g}(M)$ relation.

The observed overdensity of galaxies is also affected by lensing
magnification. Gravitational lensing by foreground structures induces a small
perturbation in the observed galaxy positions and in their observed flux. This
causes an additive contribution to $\delta_{g}$ of the form

$\delta_{g}^{M}(\hat{\bf n})=-\int
d\chi\,(2-5s)\,q_{L}(\chi)\,\Delta_{M}(\chi\hat{\bf n},z),$ (27)

where $\Delta_{M}$ is the 3D matter overdensity, $s=d\log_{10}N/dm$ is the
logarithmic slope in the distribution of galaxy magnitudes at the limiting
magnitude, and $q_{L}$ is the lensing kernel

$q_{L}(\chi)\equiv\frac{3H_{0}^{2}\Omega_{m}}{2}(1+z)\chi\int_{z}^{\infty}\,dz^{\prime}\,p(z^{\prime})\,\frac{\chi(z^{\prime})-\chi}{\chi(z^{\prime})}.$
(28)

Due to the cumulative nature of the weak lensing kernel, if relevant, the
contribution from magnification to the galaxy-CIB cross-correlation would bias
the tomographic measurement of $\langle b\rho_{\rm SFR}\rangle$ described in
the next section. To verify that this contribution can indeed be ignored, we
estimate the theoretical prediction with and without magnification for the
cross-correlation between our highest eBOSS redshift bin (for which
magnification is most relevant), and the CIB. We use $s=0.2$ for quasars
(Scranton et al., 2005). We find that magnification changes the cross-
correlation only at the sub-percent level on the scales used here, and can
therefore be ignored. Note that the impact of magnification on the DELS bins
should be significantly smaller than for the quasar samples (see Hang et al.
(2021)).

####  2.1.5 Tomography

Figure 2: Redshift distributions of the $6$ redshift bins used in this work.
The red lines show the distribution of the 4 DELS samples, while the dashed
green lines show the two eBOSS-QSO redshift bins. The dotted blue line shows
the effective radial kernel for the CIB (Eq. 13) at 545 GHz. All curves have
been rescaled by an arbitrary factor for visualisation purposes.

A crucial aspect of projected galaxy clustering is its local nature: galaxies
at redshift $z$ trace structures at the same redshift. This makes it possible
to directly measure the redshift evolution of astrophysical quantities by
cross-correlating a given projected tracer (in our case CIB intensity maps)
with galaxy samples at different redshifts, a procedure commonly called
“tomography” (Chiang et al., 2019).

Consider the Limber equation (Eq. 2), together with the simplified power
spectrum model of Eq. 7, and the radial kernels for the CIB (Eq. 13) and for
galaxy clustering (Eq. 25). Assuming that the radial/redshift dependence of
$b_{g}$, $\langle b\rho_{\rm SFR}\rangle$, and $q_{\nu}(z)$ is slow compared
with the width of the galaxy kernel $q_{g}(z)$, the auto-correlation of
galaxies and their cross-correlation with the CIB can be approximated as

$\displaystyle C_{\ell}^{gg}\simeq b_{g}^{2}\,M_{\ell}^{gg}+n_{gg},\hskip
12.0ptC_{\ell}^{g\nu}\simeq b_{g}\langle b\rho_{\rm
SFR}\rangle\,M_{\ell}^{g\nu}+n_{g\nu},$ (29)

where $b_{g}$ and $\langle b\rho_{\rm SFR}\rangle$ are evaluated at the mean
redshift of the galaxy sample, and we have defined the power spectrum
templates

$\displaystyle
M_{\ell}^{uv}\equiv\int\frac{d\chi}{\chi^{2}}q_{u}(\chi)q_{v}(\chi)P_{M}\left(\frac{\ell+1/2}{\chi},z\right),$
(30) $\displaystyle
n_{uv}\equiv\int\frac{d\chi}{\chi^{2}}q_{u}(\chi)q_{v}(\chi)N_{UV}(z),$ (31)

where $P_{M}(k,z)$ is the matter power spectrum, and $N_{UV}$ was defined in
Eq. 7. Assuming a fixed cosmology, and an accurate knowledge of the radial
kernels, $M^{gg}_{\ell}$ and $M^{g\nu}_{\ell}$ can be precomputed and treated
as fixed templates. In that case, the galaxy auto-correlation effectively
constrains the value of $b_{g}$, which then can be used to constrain $\langle
b\rho_{\rm SFR}\rangle$ at the redshift of the sample through the cross-
correlation.

Figure 2 shows the radial kernels of the different tracers used in this
analysis. The CIB kernel for the $\nu=545\,{\rm GHz}$ Planck channel, shown as
a blue dotted line (with arbitrary normalisation), varies very slowly in
comparison with the redshift distributions of the six galaxy samples used,
justifying our use of the approximate phenomenological model of Eq. 29 to make
tomographic measurements of $\langle b\rho_{\rm SFR}\rangle$.

Fig. 3 shows the cross-power spectrum between the galaxy overdensity and
$\rho_{\rm SFR}$ at $z\sim 0.8$. The solid red line shows the prediction using
the halo model from Maniyar et al. (2021). This prediction can be recovered to
$2\%$ accuracy up to $k=1\,{\rm Mpc}^{-1}$ using the simplified model of Eq.
29, given by the black dash-dotted line. This model is constructed by adding
the linear matter power spectrum scaled by the product of the scale-
independent biases associated with both quantities (dashed blue line), and the
constant large-scale limit of the 1-halo term (black dotted line). However, in
order to further improve the model and alleviate any residual theoretical
uncertainties associated with the 1-2-halo transition regime, we use the non-
linear matter power spectrum, computed using the model of Takahashi et al.
(2012), scaled by the scale-independent biases (blue dashed line). This
effective model is able to reproduce the halo model prediction to the same
accuracy on scales $k<0.2\,{\rm Mpc}^{-1}$ without the additional shot noise
terms. Therefore, we will by default use the non-linear matter power spectrum
as $P_{M}$ when calculating the templates in Eq. 30. We verified that the
results shown in Section 4.2 change by less than $0.25\sigma$ if we use the
linear power spectrum instead. Thus, the conservative choice of small-scale
cut $k_{\rm max}=0.15\,{\rm Mpc}^{-1}$ used in our analysis (see Section 4),
marked as a vertical dotted line in the figure, should allow us to obtain
unbiased constraints on $b_{g}$ and $\langle b\rho_{\rm SFR}\rangle$ from the
data.

Figure 3: Halo model power spectrum at $z=0.8$ for the galaxy-$\rho_{\rm SFR}$
cross-correlation. The solid red line shows the halo model prediction, while
the black dash-dotted line shows the simplified phenomenological model of Eq.
7, using the $k\rightarrow 0$ limit for the 1-halo term (black dotted line),
and for the bias-weighted halo profiles in the 2-halo term (black dashed
line). This approximation is accurate up to 2% on scales $k\leq 1{\rm
Mpc}^{-1}$, which justifies scale cuts used in Section 4.1. The dashed blue
line shows the simplified two-halo term using the non-linear matter power
spectrum instead of the linear one.

###  2.2 Angular power spectra and covariances

We estimate auto- and cross-power spectra from all datasets in this study
using the so-called pseudo-$C_{\ell}$ or MASTER algorithm as implemented in
NaMaster (Alonso et al., 2019). The method is described briefly here, and
further details can be found in Hivon et al. (2002); Alonso et al. (2019).

In all cases studied here (galaxy overdensity and CIB intensity maps), the
observed field $\tilde{a}(\hat{\bf n})$ is a masked version of the true
underlying field $a(\hat{\bf n})$, i.e.

$\tilde{a}(\hat{\bf n})=w(\hat{\bf n})a(\hat{\bf n}),$ (32)

where $w$ is the mask. $w$ can simply be a binary map defining the footprint
over which $a$ has been observed, but in general it can be thought of as a
local weight that can be designed to maximise the signal-to-noise ratio of the
resulting power spectrum estimator. Consider two such fields, $\tilde{a}$ and
$\tilde{b}$, with masks $v$ and $w$ respectively. Their pseudo-$C_{\ell}$ is
defined as $\tilde{C}^{ab}_{\ell}\equiv\sum_{m=-\ell}^{\ell}\tilde{a}_{\ell
m}\tilde{b}^{*}_{\ell m}/(2\ell+1)$, where $\tilde{a}_{\ell m}$ are the
harmonic coefficients of $\tilde{a}$. Due to the convolution theorem, the
$\tilde{a}_{\ell m}$ are a convolution of the harmonic coefficients of the
true map and the mask. This propagates to the power spectra, so that the
expectation value of $\tilde{C}_{\ell}^{ab}$ is

$\langle\tilde{C}^{ab}_{\ell}\rangle=\sum_{\ell^{\prime}}M^{vw}_{\ell\ell^{\prime}}C^{ab}_{\ell^{\prime}}.$
(33)

Here, $M^{vw}_{\ell\ell^{\prime}}$ is the mode-coupling matrix (MCM), which
depends only on the masks $v$ and $w$, and can be calculated analytically with
relative ease. In order to estimate the true underlying power spectra from the
observations we then need to invert the MCM. Since the MCM is not generally
invertible, the pseudo-$C_{\ell}$ estimator often involves binning
$\tilde{C}^{ab}_{\ell}$ into bandpowers before inverting the binned MCM. Here,
we used the same bandpowers used in García-García et al. (2021). Linear $\ell$
bins with $\Delta\ell=30$ were used up to $\ell=240$, and logarithmic bins
were used at higher $\ell$ with $\Delta\log_{10}\ell=0.055$.

The covariance matrix of angular power spectra is sourced by three
contributions: a disconnected (or “Gaussian”) part, a super-sample covariance
(SSC) due to modes larger than the survey footprint, and a connected non-
Gaussian (cNG) part caused by the intrinsic non-Gaussianity of the tracers
(Takada & Hu, 2013). On the large scales used here the Gaussian component
dominates, and we neglect the SSC and cNG contributions. We use an approximate
analytical method to compute the Gaussian covariance matrix accounting for the
effects of survey geometry, the so-called “narrow-kernel approximation” (NKA)
described in García-García et al. (2019), and implemented in NaMaster. The NKA
reduces the scaling of the exact calculation from $O(\ell_{\rm max}^{6})$ to a
more tractable $O(\ell_{\rm max}^{3})$ by assuming that the MCM is narrowly
peaked at the diagonal compared with the slow scale dependence of the true
power spectrum.

Schematically, the covariance between two pseudo-$C_{\ell}$s can be
approximated as

${\rm Cov}(\tilde{C}^{ab}_{\ell},\tilde{C}^{fg}_{\ell^{\prime}})={\cal
C}^{ag}_{(\ell,\ell^{\prime})}{\cal
C}^{bf}_{(\ell,\ell^{\prime})}\Xi_{\ell\ell^{\prime}}(w^{a}w^{g},w^{b}w^{f})+(a\leftrightarrow
b)\,$

where $\Xi_{\ell\ell^{\prime}}(w^{a}w^{g},w^{b}w^{f})$ are coupling
coefficients similar to those involved in the calculation of
$M^{vw}_{\ell\ell^{\prime}}$, but now dependent on the products of pairs of
masks (labelled $w^{a/b/f/g}$ here). The modified spectra ${\cal
C}^{ab}_{\ell}$ are related to the true underlying power spectra via (Nicola
et al., 2021)

${\cal C}^{ab}_{\ell}\equiv\frac{\langle\tilde{C}^{ab}_{\ell}\rangle}{\langle
w^{a}\,w^{b}\rangle_{\rm pix}},$ (34)

where $\langle\cdots\rangle_{\rm pix}$ denotes averaging over all pixels in
the sky (as opposed to ensemble averaging). In general
$\langle\tilde{C}^{ab}_{\ell}\rangle$ can be computed from a theoretical model
and convolved with the MCM. However, since the power spectra estimated here
have a high signal-to-noise ratio, we simply use the estimated
pseudo-$C_{\ell}$s in the formula above (i.e. replace
$\langle\tilde{C}_{\ell}^{ab}\rangle\rightarrow\tilde{C}_{\ell}^{ab}$). In
addition to simplicity, this has the added benefit of automatically accounting
for the noise contribution to the power spectrum with no additional modelling.

This estimate of the covariance matrix was validated against jackknife
resampling. As described in Section 4.1, both estimates are in reasonable
agreement, although a small modification was applied to the final covariance
matrix elements involving two different CIB frequency maps to accurately
account for their map-level correlation.

###  2.3 Likelihood

We derive parameter constraints from the power spectrum data in two different
ways. First, as described in Section 4.2, we derive model-independent
tomographic measurements of $\langle b\rho_{\rm SFR}\rangle$ from the galaxy
auto-correlations and their cross-correlation with the CIB. We will then use
these measurements, as described in Section 4.3, to place constraints on the
free parameters of the halo models for the SFR described in Section 2.1.3. In
both cases, we will make use of Gaussian likelihoods of the form

$-2\log p({\bf d}|\Theta)=({\bf d}-{\bf m}(\Theta))^{T}{\sf C}_{\bf
d}^{-1}({\bf d}-{\bf m}(\Theta))+K,$ (35)

where ${\bf d}$ is the data vector (power spectra or measurements of $\langle
b\rho_{\rm SFR}\rangle$), ${\sf C}$ is the covariance of ${\bf d}$, ${\bf m}$
is the theoretical model for ${\bf d}$, and $\Theta$ denotes the set of
parameters used to define ${\bf m}$. The posterior parameter distribution is
then given by the product of this likelihood and the parameter priors, which
we discuss in Section 4.2. In all cases we will use the affine-invariant
Markov-Chain Monte-Carlo (MCMC) algorithm implemented in emcee (Foreman-Mackey
et al., 2013).

The assumption of a Gaussian likelihood for the power spectra stems from the
central limit theorem: for sufficiently high $\ell$, the estimate of
$C_{\ell}$ involves averaging over a large number of independent $m$ modes,
and thus the distribution of $C_{\ell}$ values tends to a Gaussian. This is
valid on the scales used in this analysis Hamimeche & Lewis (2008). The
validity of the Gaussian approximation for the likelihood of the measured
$\langle b\rho_{\rm SFR}\rangle$ values will be demonstrated in Section 4.2.

All cosmological theory predictions were computed using the Core Cosmology
Library (CCL, Chisari et al. (2019)222The source code can be found at
https://github.com/LSSTDESC/CCL). Cosmological parameters were fixed to the
Planck best-fit values
$(\Omega_{c},\Omega_{b},h,n_{s},\sigma_{8})=(0.261,0.049,0.677,0.9665,0.8102)$.
When using the halo model, we use the halo mass function parametrisation of
Tinker et al. (2008), and the halo bias of Tinker et al. (2010). Halo masses
are defined using a spherical overdensity $\Delta=200$ times the critical
density.

##  3 Data

Figure 4: Sky footprint of the CIB and galaxy surveys used in this analysis in
Equatorial coordinates.

###  3.1 Galaxy clustering

Bin | Redshift range | Mean redshift | Density (${\rm deg}^{-2}$)
---|---|---|---
DELS
1 | [$0.10,0.30$) | $0.21$ | 808
2 | [$0.30,0.15$) | $0.37$ | 651
3 | [$0.45,0.60$) | $0.50$ | 760
4 | [$0.60,0.80$] | $0.63$ | 409
eBOSS
5 | [$0.80,1.50$) | $1.12$ | 34
6 | [$1.50,2.20$] | $1.87$ | 35
Table 1: Mean redshift values and approximate redshift ranges for the $6$
redshift bins shown in Figure 2. For DELS, the bin edges are in photo-$z$
space, while the mean redshift corresponds to the mean of the true redshift
distribution. The 4th column shows the angular number density of each sample.

In order to recover the star formation history in the range $0\lesssim
z\lesssim 2$ we make use of photometric galaxy samples from the DESI Legacy
Imaging Survey (Dey et al., 2019, DELS), and quasar samples (QSOs) from the
extended Baryon Oscillation Spectroscopic Survey (eBOSS), derived the from
Sloan Digital Sky Survey Data Release 16 (DR16).

At low redshifts, we use the photometric sub-sample of DELS selected by Hang
et al. (2021)333The data can be found at https://gitlab.com/qianjunhang/desi-
legacy-survey-cross-correlations. We thank the authors for the remarkable care
with which these data were made publicly available.. DELS combines photometry
from three different telescopes: the DECam Legacy Survey (Flaugher et al.,
2015; Blum et al., 2016) at declinations ${\rm DEC}<33^{\circ}$, the Mayall
$z$-band Legacy Survey (Silva et al., 2016), and the Beijing-Arizona Sky
Survey (Zou et al., 2019). The final survey covers 17739 deg2. In Hang et al.
(2021), each galaxy was assigned a photometric redshift based on a multi-
dimensional matching in colour space with a set of spectroscopic samples. We
use these redshifts to separate the sample into the four tomographic bins used
in Hang et al. (2021), following the same procedure described there to model
the redshift distributions of the resulting samples. Moreover, we impose low-
declination cut of ${\rm DEC}>-36^{\circ}$, discarding a region exhibiting
significant variations in completeness.

At higher redshifts we use the homogeneous eBOSS quasar sample used for the
cosmological power spectrum analysis of Neveux et al. (2020); Hou et al.
(2021), and presented in Ross et al. (2020). The catalogue comprises $343708$
objects with measured redshifts in the range $0.8\leq z\leq 2.2$, covering
over $4800$ deg2. We combine the measurements from the North and South
Galactic Caps into one single catalogue that we split into two different bins
with redshifts above and below $z=1.5$ respectively. The redshift distribution
of each bin is estimated directly from the data as a histogram of the measured
spectroscopic redshifts.

We construct maps of the galaxy overdensity from both samples using the
methods described in García-García et al. (2021). These include the
construction of galaxy masks, the use of random catalogues for eBOSS, the
correction of large-scale systematic fluctuations for DELS, and the estimation
of the noise bias for both samples. Details can be found in Section 3 of
García-García et al. (2021). The impact of large-scale systematics in the
galaxy auto-correlation is minimal given the scale cuts used for this analysis
(see Section 4.2). Likewise, our analysis is not sensitive to inaccuracies in
the estimation of the shot-noise contribution, since we fully marginalise over
a noise-like white component ($n_{xy}$ in Eq. 29) for all power spectra
analysed. All maps are constructed and processed using HEALPix (Górski et al.,
2005), using a resolution parameter $N_{\rm side}=2048$, corresponding to a
pixel size of $\sim 1.7^{\prime}$. Our final galaxy sample thus comprises 6
different redshift bins, shown in Fig. 2, and described in Table 1.

###  3.2 CIB maps

We use the full mission CIB maps and associated sky masks constructed by Lenz
et al. (2019). The maps were created from the Planck 353, 545 and 857 GHz
intensity maps. Galactic neutral hydrogen data from the HI4PI Survey (Bekhti
et al., 2016) were used to create template maps of galactic dust. These
templates were then used to subtract potential galactic dust contamination in
the maps.

We apply the sky mask constructed by Lenz et al. (2019), based on the $20$%
Planck Galactic plane mask, removing small-scale Galactic filamentary
structure via HI column densities above $N_{\rm HI}>2.5\times 10^{20}\,{\rm
cm}^{-2}$. The combined mask was then apodised using a $15^{\prime}$ FWMH
kernel.

When estimating the galaxy-CIB cross-correlations, we de-convolve the
effective window function provided with the data products of Lenz et al.
(2019). This accounts for the Planck Reduced Instrument Model (RIMO) beam, the
window function of the SMICA CMB map, and the pixel window function for
$N_{\rm side}=2048$ (see section 2.2 of Lenz et al. (2019) for further
details).

As discussed in Lenz et al. (2019), the use of small sky patches to remove the
Galactic dust component from the Planck maps leads to an over-subtraction of
the underlying CIB fluctuations that results in a lack of power on large
scales. The effect is limited to scales $\ell\lesssim 70$, and therefore we
will omit all such scales from our analysis. We applied color corrections as
described in Planck Collaboration et al. (2014b).

The footprints of the three datasets used here are shown in Fig. 4. The CIB
maps overlap with the DELS and eBOSS footprints over $\sim 13\%$ and $\sim
6\%$ of the sky respectively.

##  4 Results

###  4.1 Power spectra and covariance

Figure 5: Galaxy auto-correlation spectra (left column) and galaxy-CIB cross-
correlation spectra (right column) for the six tomographic bins presented in
Fig. 2. Each panel shows the measurements and their uncertainties (black
points with error bars), and their best-fit predictions (red lines), which are
extrapolated outside of the range of scales used in our analysis. Our fiducial
scale cuts are marked by the light grey vertical bands, while the dark grey
bands mark the alternative, less conservative scale cuts corresponding to
$k_{\rm max}=0.2\,{\rm Mpc}^{-1}$. The bottom part of each panel shows the
difference between the data and the best-fist curve normalised by the
$1\sigma$ uncertainties. Figure 6: Correlation matrix $(r_{i,j})_{{}_{1\leq
i,j\leq 200}}$ of the total data vector. The six blocks along the diagonal
correspond to the different redshift bins used in the analysis. Within these
blocks, the first sub-block contains the covariance of the galaxy auto-
correlation, while the other blocks correspond to the different cross-
correlations with the CIB. The different CIB frequency channels are strongly
($\sim 95\%$) correlated. Figure 7: Top: variance of the galaxy-CIB cross-
correlation for the fourth DELS redshift bin with the three CIB frequencies
(solid, dashed, and dotted lines respectively). Results are shown in red for
the analytical covariance, and in black for the jackknife estimate. Bottom:
ratio between jackknife and analytical estimates for the correlation
coefficient of pairs of cross-power spectra involving different CIB
frequencies. Results are shown for the $(353-545)$ GHz, $(545-857)$ GHz and
$(353-857)$ GHz pairs in red, blue and yellow respectively. We find that, on
average, the jackknife estimate predicts a $\sim 2\%$ lower correlation
coefficient between frequencies. We correct our analytical estimate for this
factor (shown as a horizontal dashed line).

We measure the galaxy auto-correlation power spectrum for the six redshift
bins considered, as well as their cross-correlations with the three CIB
frequency maps, for a total of 24 auto- and cross-spectra. As an example, the
galaxy clustering auto-correlations and the cross-correlations with the $545$
GHz CIB map are shown in Fig. 5. The measured power spectra are shown in black
together with their statistical uncertainties. The light grey bands show the
fiducial scale cuts used in our analysis. We impose a default low-$\ell$ cut
of $\ell_{\rm min}=100$. This is mostly motivated by the loss of large-scale
power in the CIB maps caused by the removal of Galactic contamination on
different large-area patches, as reported by Lenz et al. (2019). Although this
should not affect the galaxy auto-correlation, observational systematics
leading to artificial contamination of the galaxy number counts across the sky
may lead to extra power on the largest scales. $\ell_{\rm min}=100$ is likely
a conservative choice in this case (Hang et al., 2021; García-García et al.,
2021), but we choose to keep the same cut for both power spectra for
simplicity. On small scales, we remove all multipoles above $\ell_{\rm
max}=k_{\rm max}\bar{\chi}-1/2$, where $\bar{\chi}$ is the radial comoving
distance to the mean redshift of each galaxy sample, and $k_{\rm
max}=0.15\,{\rm Mpc}^{-1}$. As we argued in Section 2.1.5, we expect the
phenomenological model of Eq. 29 to provide a good fit to the data on these
scales. To test the robustness of our results to this choice, we will also
present results for a less conservative cut of $k_{\rm max}=0.2\,{\rm
Mpc}^{-1}$ (shown by the dark grey bands in the Figure). After these scale
cuts, the complete data vector, comprising six galaxy auto-correlations and
6$\times$3 galaxy-CIB cross-correlations contains $200$ elements.

Figure 5 also shows, in red, the best-fit theoretical prediction following the
model of Eq. 29. The residuals between data and prediction normalised by the
diagonal errors are shown at the bottom of each panel. As discussed and
quantified in the next section, overall we find that the model is able to
provide a good fit to the data.

Figure 6 shows the correlation matrix, defined in terms of the covariance
$C_{ij}$ as $r_{ij}\equiv C_{ij}/\sqrt{C_{ii}C_{jj}}$. The blocks along the
diagonal correspond to the six different redshift bins, with the galaxy auto-
correlation followed by the 3 CIB cross-correlations. In each block, the
cross-covariance between $C_{\ell}^{gg}$ and the $C_{\ell}^{g\nu}$s is
significantly smaller than the cross-covariance between CIB channels. The
latter is at the level of $\sim 95\%$. This is because the three Planck maps
trace essentially the same underlying CIB fluctuations with little frequency
decorrelation. This makes an accurate calculation of the full covariance
crucial since this high level of correlation requires the model residuals
between different frequencies to follow the same trend almost exactly. The
overlaps between the DELS redshift distributions shown in Figure 2 also lead
to an appreciable cross-covariance between the power spectra in adjacent
redshift bins.

We validate the analytical estimate of the covariance matrix described in
Section 2.2 through jackknife resampling, dividing the footprint into 413
jackknife regions. In general, we find good agreement between both estimates
of the covariance matrix, and we are able to recover the diagonal errors
within $\sim 4-5\%$ (see top panel of Fig. 7). However, we find that the
jackknife estimate recovers cross-correlation coefficients between power
spectra involving different CIB frequencies that are consistently lower than
the analytical prediction by $\sim 2\%$ (see bottom panel of Fig. 7). Although
this is a small difference, the tight correlation between different CIB
frequencies makes it necessary to accurately characterise these off-diagonal
covariance elements. We find that, otherwise, the goodness of fit to the full
dataset combining all frequencies was unacceptable, even if the same model
parameters were able to provide a good fit to the three different frequencies
separately, which clearly points to an inconsistent covariance matrix. To
correct for this, we therefore multiply all off-diagonal covariance matrix
elements involving two different CIB frequencies by an overall factor of 0.98.
The difference is likely due to inhomogeneity in the map noise across the
footprint, which is not accounted for in the analytical covariance.

###  4.2 Model-independent measurements of $\langle b\rho_{\rm SFR}\rangle$

Figure 8: Inferred values of $\langle b\rho_{\rm SFR}\rangle$ as a function of redshift for different analysis choices, with their $1\sigma$ error bars. Our fiducial measurements are shown in black. The results assuming less restrictive scale cuts are shown in blue, and the measurements carried out in each redshift bin independently are shown in grey. The constraints found from different frequency maps independently are shown in orange, red, and burgundy (for 353, 545, and 857 GHz respectively). The yellow and green bands show the prediction and 1-$\sigma$ confidence interval for $\langle b\rho_{\rm SFR}\rangle$ obtained from the two halo models described in Section 2.1.3 (Model A and Model B respectively). Figure 9: Posterior probability distributions for our measurements of $\langle b\rho_{\rm SFR}\rangle$ in each redshift bin. The 1- and 2-$\sigma$ contours calculated from the MCMC chains are shown in red. The blue dash-dotted line shows the same contours for a multivariate normal distribution with the same mean and covariance as the chains. This justifies the use of a Gaussian likelihood when using these measurements to place constraints on halo-based models of the SFR in Section 4.3. |  | $k_{\rm max}=0.15\,{\rm Mpc}^{-1}$ |  |
---|---|---|---|---
Bin | $\bar{z}$ | $b_{g}$ | $\langle b\rho_{\rm SFR}\rangle$ |
($M_{\odot}{\rm yr}^{-1}{\rm Mpc}^{-3}$) |
DELS
1 | 0.21 | $1.068\pm 0.044$ | $-0.005\pm 0.008$ |
2 | 0.37 | $1.382\pm 0.026$ | $0.026\pm 0.006$ |
3 | 0.50 | $1.306\pm 0.017$ | $0.041\pm 0.005$ |
4 | 0.63 | $1.733\pm 0.019$ | $0.048\pm 0.005$ |
eBOSS
5 | 1.12 | $2.057\pm 0.149$ | $0.110\pm 0.014$ |
6 | 1.87 | $2.165\pm 0.191$ | $0.260\pm 0.032$ |
$\chi^{2}/N_{\rm d.o.f}=201.2/188=1.07,\hskip 12.0pt{\rm PTE}=0.245$ |
Table 2: Summary of the obtained values for the galaxy bias and the bias
weighted star-formation rate within each redshift bin as presented in Figure
2, when all parameters are sampled jointly. We display the mean values and the
$1\sigma$ uncertainties. The last row show the $\chi^{2}$ per degree of
freedom value and PTE for our constraints.

We obtain tomographic measurements of the bias-weighted mean SFR density,
$\langle b\rho_{\rm SFR}\rangle$, by fitting the phenomenological model in Eq.
29 to the power spectrum data described in the previous section. The model
consists on one galaxy bias parameter $b_{g}$ and one value of $\langle
b\rho_{\rm SFR}\rangle$ for each redshift bin, as well as a free white-noise
amplitude ($n_{xy}$) for each individual power spectrum. Although the white
noise levels for different CIB frequencies should be tightly correlated, we
choose to allow each of them to vary freely, in order to account for different
levels of point-source subtraction, and non-Poissonian stochastic
contributions in galaxy clustering. The full model therefore has a total of
$6\times N_{\rm bin}=36$ free parameters. Note, however, that since the white
noise amplitudes are linear parameters, it is possible to marginalise over
them analytically, significantly reducing the number of degrees of freedom we
have to explore to $2\times N_{\rm bin}=12$. This is achieved by modifying the
inverse covariance matrix used in the Gaussian likelihood of Eq. 35 as
(Rybicki & Press, 1992)

$\left({\sf C}_{\bf d}^{-1}\right)^{\prime}={\sf C}_{\bf d}^{-1}-{\sf C}_{\bf
d}^{-1}{\sf T}\left({\sf T}^{T}{\sf C}_{\bf d}^{-1}{\sf T}\right)^{-1}{\sf
T}^{T}{\sf C}_{\bf d}^{-1},$ (36)

where ${\sf T}$ is a matrix containing, as columns, the templates
corresponding to the different linear parameters of the model (in this case, a
constant vector for each of the different power spectra included in the data
vector). We used flat, uninformative priors on all parameters.

Figure 8 shows the resulting measurements of $\langle b\rho_{\rm SFR}\rangle$.
Our fiducial measurements, simultaneously using all redshift bins and CIB
frequencies, with a small-scale cut of $k_{\rm max}=0.15\,{\rm Mpc}^{-1}$ are
shown in black. The associated numerical values can be found in Table 2. We
obtain $\sim 10\%$ constraints on $\langle b\rho_{\rm SFR}\rangle$ at
$z\gtrsim 0.5$, and a $\sim 20\%$ measurement in the second redshift bin
($z\simeq 0.37$). The conservative scale cuts used here leave only 3 data
points in the cross-correlation with the first redshift bin, and we are only
able to provide an upper bound on $\langle b\rho_{\rm SFR}\rangle$, due to the
degeneracy between this parameter and the white noise amplitude. The model is
able to reproduce the data with a good $\chi^{2}$, with a probability-to-
exceed ${\rm PTE}=0.245$. The full posterior distribution of the $\langle
b\rho_{\rm SFR}\rangle$ parameters is shown in Fig. 9. Although the
measurements are largely statistically independent from one another, we can
appreciate a non-negligible correlation between the measurements in adjacent
DELS bins (e.g. bins 2-3 and 3-4). More importantly, we find that the
posterior distribution is highly Gaussian. The dot-dashed black lines in the
Figure show a multivariate Gaussian distribution with the same mean and
covariance as that obtained from the MCMC chains, which reproduces the true
posterior almost exactly. This will enable us to use the measured values of
$\langle b\rho_{\rm SFR}\rangle$ to constrain halo-based models of the star
formation history in Section 4.3. For convenience, we provide the correlation
matrix of these measurements:

${\sf r}=\left(\begin{array}[]{cccccc}1.00&0.14&-0.01&0.01&0.02&-0.03\\\
0.14&1.00&0.26&0.02&-0.02&-0.01\\\ -0.01&0.26&1.00&0.27&0.01&0.00\\\
0.01&0.02&0.27&1.00&0.02&-0.01\\\ 0.02&-0.02&0.01&0.02&1.00&0.06\\\
-0.03&-0.01&0.00&-0.01&0.06&1.00\end{array}\right).$ (37)

To validate these measurements, we repeat our analysis changing the number of
redshift bins and CIB frequencies used, as well as the small-scale cut. The
results for each of these alternative analyses are listed in Tables 3, 5, and
4 of Appendix A. First, in order to quantify the effect of the overlap between
different redshift bins, we obtain new measurements of $\langle b\rho_{\rm
SFR}\rangle$ for each bin separately. We recover constraints that are almost
equal to our fiducial measurements, and the model is able to reproduce the
data for each redshift bin individually, with $\chi^{2}$ $p$-values above
$0.1$. These results are shown as grey points in Fig. 8.

The model used here relies on having an accurate description of the mean
infrared spectral energy distribution ($S^{\rm eff}_{\nu}$ in Eq. 13), and
inaccuracies in it would lead to inconsistent measurements of $\langle
b\rho_{\rm SFR}\rangle$ for different frequencies. To test for this, we repeat
our measurement of $\langle b\rho_{\rm SFR}\rangle$ three more times, using
only one of the three CIB frequencies at a time. The results are presented in
Table 5, and as orange/red/brown points in Fig. 8 for the 353/545/857 GHz
channels respectively. Overall we find good agreement between the three
frequencies. The 353 GHz channel achieves consistently larger uncertainties
than the other two, and the recovered values of $\langle b\rho_{\rm
SFR}\rangle$ are also consistently higher, although compatible with the
fiducial measurements within $1\sigma$. Although this could be a hint for a
small inconsistency between the CIB data and the effective spectra of
Béthermin et al. (2013, 2015, 2017), it may also be an indication of
contamination by other extragalactic foregrounds at lower frequencies,
particularly the thermal Sunyaev-Zel’dovich effect. Since we do not observe a
definite trend in the difference between the 353 GHz constraints and the rest
as a function of redshift, and since the differences lie within the 68%
uncertainties, we did not investigate this further, and proceeded to combine
all frequencies in our fiducial analysis.

Finally we repeat the analysis relaxing the small-scale cut to $k_{\rm
max}=0.2\,{\rm Mpc}^{-1}$, which allows us to include one or two additional
bandpowers in each redshift bin (see Fig. 5). Overall, the recovered values of
$\langle b\rho_{\rm SFR}\rangle$ are again in good agreement with our fiducial
results, within $1\sigma$ uncertainties. With these additional data points,
the uncertainties on $\langle b\rho_{\rm SFR}\rangle$ reduce by $\sim
20-30\%$. More importantly, the extra data point for the first DELS redshift
bins allows us to break the degeneracy with the white noise amplitude and
obtain a measurement of $\langle b\rho_{\rm SFR}\rangle$ at $z\simeq 0.2$. The
value recovered is in reasonable agreement with the predictions of the best-
fit halo models derived in the next Section when extrapolated to lower
redshifts. The model is still able to provide a reasonable fit to the data
with this less conservative scale cut (${\rm PTE}=4\%$). Our measurements of
$\langle b\rho_{\rm SFR}\rangle$ are therefore robust to the choice of scale
cut made for our analysis.

###  4.3 Constraints on halo-based SFR models

Figure 10: Constraints on the two halo-based SFR models (Model A and Model B
in the left and right panels, respectively). The red bands in the left panel
show the $1\sigma$ and $2\sigma$ constraints on Model A obtained by Yan et al.
(2022). The red lines in the right panel show the best-fit parameters obtained
by Moster et al. (2018) for Model B. Figure 11: Star-formation rate density as
a function of redshift. The points with error bars correspond to the direct
measurements from the infrared luminosity function, extracted from Madau &
Dickinson (2014a); Davies et al. (2016b). The solid black shows the best-fit
prediction from the CIB analysis of Maniyar et al. (2021). The two bands show
the 1-$\sigma$ constraints obtained from our measurements of $\langle
b\rho_{\rm SFR}\rangle$ using the two halo-based SFR models described in
Section 2.1.3 (Model A in yellow, Model B in green).

We use the tomographic measurements of $\langle b\rho_{\rm SFR}\rangle$
presented in the previous section to derive constraints on the two halo models
for the SFR history described in Section 2.1.3. For this, we use again the
Gaussian likelihood of Eq. 35, where ${\bf d}$ is the set of six measured
values of $\langle b\rho_{\rm SFR}\rangle$, ${\sf C}_{\bf d}$ is the
covariance of these measurements estimated from the MCMC chains, and ${\bf m}$
is the theoretical model, given by Eq. 14. The model has 4 free parameters for
both models:

  * •

Model A has $\Theta=\\{\eta_{\rm max},\log_{10}M_{\rm
max},\sigma_{M,0},\tau\\}$ with the following top-hat priors:

$\displaystyle 0\leq\eta_{\rm max}\leq 1.2,$ (38) $\displaystyle
10.5\leq\log_{10}M_{\rm max}\leq 14,$ (39) $\displaystyle
0\leq\sigma_{M,0}\leq 4,$ (40) $\displaystyle 0\leq\tau\leq 3.$ (41)

  * •

Model B has $\Theta=\\{m_{0},m_{z},\eta_{0},\eta_{z}\\}$ with the following
top-hat priors:

$\displaystyle 9\leq m_{0}\leq 14,$ (42) $\displaystyle-6\leq m_{z}\leq 6,$
(43) $\displaystyle 0\leq\eta_{0}\leq 1,$ (44) $\displaystyle
0\leq\eta_{z}\leq 1.$ (45)

Since aiming to constrain 4 parameters with only 6 data points will likely
lead to overfitting, we do not expect to be able to obtain competitive
constraints on all of these parameters. Our main aim, however, is to use these
models to map our measurements of $\langle b\rho_{\rm SFR}\rangle$ onto
constraints on the evolution of the star-formation rate density $\rho_{\rm
SFR}(z)$.

The posterior distribution of the parameters of Model A is shown in Fig. 10
(left panel). Our constraints are in good agreement with those found by Yan et
al. (2022) using the same CIB maps and galaxy clustering data from KiDS. These
are shown as red vertical bands in the diagonal panels of the figure. Our
measurements seem to favour slightly lower values of $\tau$ than those of Yan
et al. (2022), albeit with a long tail towards large $\tau$ that makes them
broadly compatible. The peak efficiency $\eta_{\rm max}\sim 0.4$ is found to
be in good agreement with both Yan et al. (2022) and the CIB constraints from
Maniyar et al. (2021). On the other hand, our data favour a peak star
formation mass of $\log_{10}M_{\rm max}\sim 11.5$, significantly lower than
that found by Maniyar et al. (2021) ($\log_{10}M_{\rm max}=12.94\pm 0.02$),
but compatible with the constraints of Yan et al. (2022).

The posterior distribution for Model B is shown in Fig. 10 (right panel). Our
constraints are in good agreement with the best-fit parameters of Moster et
al. (2018)444The constraints found by Moster et al. (2018) using a large
dataset are much tighter than those found here, and thus we only show the
final best-fit values., shown as dashed red lines. The data are able to
constrain the peak mass parameters to

$m_{0}=11.50\pm 0.85,\hskip 12.0ptm_{z}=0.5\pm 1.3,$ (46)

although the peak efficiency parameters are largely unconstrained within their
prior due to the degeneracy with $m_{0}$ and $m_{z}$. The peak mass is roughly
constant, with $\log_{10}M_{\rm max}\simeq 11.6\pm 0.8$ within the redshift
range of our data. This is consistent with the results found with Model A.

In order to place constraints on the star-formation rate history from our
measurements of $\langle b\rho_{\rm SFR}\rangle$, we randomly select 1000
samples from the MCMC chains for both models and calculate $\langle b\rho_{\rm
SFR}\rangle$ and $\rho_{\rm SFR}$ as a function of redshift via Eqs. 14 and 15
respectively. The yellow and green bands in Fig. 8 show the $68\%$ confidence
interval for the predicted evolution of $\langle b\rho_{\rm SFR}\rangle$
extracted from Model A and Model B respectively. The corresponding constraints
on the star-formation rate density are shown in Fig. 11. The figure also shows
various direct measurements of $\rho_{\rm SFR}$ from the infrared luminosity
function (Shang et al., 2012) (points with error bars), as well as the
prediction from the best-fit model found by Maniyar et al. (2021) from the CIB
auto-correlation. The direct $\rho_{\rm SFR}$ measurements are from Sanders et
al. (2003); Takeuchi et al. (2003); Magnelli et al. (2011); Magnelli et al.
(2013); Gruppioni et al. (2013); Marchetti et al. (2016) and Davies et al.
(2016a), and were extracted from Madau & Dickinson (2014a) and Davies et al.
(2016b). Overall, our constraints are in rough agreement with these
independent direct measurements with two main caveats. At low redshifts, our
constraints on Model A predict a marginally higher SFR density compared with
the direct measurements, whereas the predictions of Model B agree well with
the data. This may be due to the fact that, unlike Model A, Model B allows for
both the peak efficiency and the maximum efficiency mass to depend on
redshift. The conversion from $\langle b\rho_{\rm SFR}\rangle$ to $\rho_{\rm
SFR}$ at low-redshift prediction is therefore somewhat sensitive to the
details of the underlying model. On the other hand, at high redshift both
models predict a marginally higher $\rho_{\rm SFR}$ compared with the direct
measurements. This is likely driven by the large amplitude of the cross-
correlation with the high-redshift eBOSS quasar sample, and therefore cross-
correlations with other high-redshift tracers would be valuable to confirm
this result.

These caveats aside, the fact that our measurements of $\langle b\rho_{\rm
SFR}\rangle$ are compatible with direct measurements of the SFR density is an
interesting and non-trivial result, beyond the fact that both constraints come
from independent datasets. The contribution of haloes of different masses to
$\langle b\rho_{\rm SFR}\rangle$ is weighted by their linear bias, and thus
$\langle b\rho_{\rm SFR}\rangle$ is more sensitive than $\rho_{\rm SFR}$ to
star formation in the most massive haloes. An agreement between measurements
of both quantities within a given model thus provides non-trivial support for
the physical interpretation of the dependence of star formation on halo mass.
The measurements of $\langle b\rho_{\rm SFR}\rangle$ presented here, and those
that will be achievable with future datasets, can therefore be used to
strengthen constraints on star formation and galaxy evolution models, two
areas of astrophysics of critical importance for cosmology in the next decade.

##  5 Conclusions

We have presented a tomographic analysis of Cosmic Infrared Background maps
through cross-correlation with a suite of six galaxy samples spanning the
redshift range $z\lesssim 2$ comprising galaxies from the DESI Legacy Survey
and quasars from eBOSS. The CIB is a good tracer of the obscured star
formation history, and thus the main motivation behind this exercise was to
place direct constraints on the evolution of star formation in the last $\sim
10$ billion years. Our analysis has been carried out in two steps.

First, we used a phenomenological parametrisation of the angular power
spectra, inspired on the halo model, which allows us to directly measure the
mean star-formation rate density weighted by halo mass, $\langle b\rho_{\rm
SFR}\rangle$, in a virtually model-independent manner (see Eq. 29). Although
this involves restricting our analysis to the largest scales, thus discarding
significant amounts of information, we have argued that the resulting
constraints are robust to a wide range of theoretical uncertainties that a
full modelling of the small-scale $C_{\ell}$s would be sensitive to. These
include accurately describing the scale dependence of the mean and variance of
the galaxy overdensity in haloes, the covariance between galaxy overdensity
and SFR density in haloes (or, equivalently, characterising star formation in
the target galaxy sample), the transition between the 1-halo and 2-halo-
dominated regimes, and the different contributions (Poisson or otherwise) to
the shot-noise component of all auto- and cross-correlations involved.

We find that this phenomenological model provides a good fit to the data on
the scales probed. We thus obtained measurements for $\langle b\rho_{\rm
SFR}\rangle$ at six different redshifts, reported in Table 2. The measurements
are robust with respect to the choice of scale cuts, as well as the
combination of datasets used (i.e. the choice of redshift bins and frequency
maps). Furthermore, we have shown that the posterior distribution for these
$\langle b\rho_{\rm SFR}\rangle$ measurements is well approximated by a
multivariate normal. It is therefore easy to use these measurements to derive
constraints on any halo-based parametrisation of star formation using a simple
Gaussian likelihood. Since $\langle b\rho_{\rm SFR}\rangle$ is more sensitive
to the contribution to the SFR density from massive haloes than direct
measurements of $\rho_{\rm SFR}$, a combination of both measurements can be
used to put constraints on the halo mass dependence of star formation.

As a second step, we use the measured values of $\langle b\rho_{\rm
SFR}\rangle$ as data to constrain the parameters of two halo-based models for
the SFR history, characterising the mass and redshift dependence of the
efficiency with which gas is converted into stars (see Fig. 1). Besides
allowing us to place constraints on the models themselves, this procedure
makes it possible to map our constraints on $\langle b\rho_{\rm SFR}\rangle$
onto constraints on $\rho_{\rm SFR}$ itself. The result is shown in Fig. 11.
Our constraints are in good qualitative agreement with independent
measurements of the SFR density. This supports the physical interpretation of
the halo-based SFR models of Moster et al. (2018); Maniyar et al. (2021), and
the role of the CIB as a cosmic tracer of star formation. We have shown,
however, that the predicted SFR density at low redshifts depends on the
details of the underlying model, whereas our data seem to consistently predict
a slightly higher $\rho_{\rm SFR}$ at $z\sim 2$. Thus our measurements of
$\langle b\rho_{\rm SFR}\rangle$ contain information that is complementary to
that of direct measurements of the SFR density, and which can be used to
improve our understanding of star formation in haloes.

A number of caveats associated with our analysis must be noted. First and
foremost, as Fig. 5 shows, the cross-correlation between the CIB and our
galaxy samples is detected at very high significance over all scales and
redshifts explored. Nevertheless, the simplified model used here requires us
to discard most of these data. A more sophisticated model would allow us to
derive significantly tighter constraints on both models of the star formation
history, and on the astrophysical properties of the galaxy samples used. As we
have noted above, several major modelling challenges must be overcome in order
to construct such a model given the high statistical significant of the
measurement. Ensuring that the model is able to recover unbiased constraints
on its parameters would likely require the use of simulations incorporating a
somewhat realistic description of star formation and galaxy evolution. The
rewards of developing more sophisticated models describing this and other
cross-correlations (Hurier et al., 2015; Hurier et al., 2017; Krolewski et
al., 2021; McCarthy & Madhavacheril, 2021; Krolewski & Ferraro, 2022; Tröster
et al., 2022) are nevertheless worth the effort. For now, we have limited
ourselves to providing these conservative but robust constraints, and leave a
more complex analysis, modelling the small-scale power spectra directly for
future work.

Our model also relies on a number of simplified assumptions. First, we have
assumed perfect knowledge of the mean infrared spectral energy distributions
as tabulated by Béthermin et al. (2013, 2015, 2017). Since our measurement is
based on the large-scale bias extracted from the 2-halo term, it is
insensitive to the scatter in this spectrum. However, any mis-modelling of the
redshift dependence of the mean SED would directly lead to a biased recovery
of the evolution of $\langle b\rho_{\rm SFR}\rangle$ (although the fact that
we recover compatible results from the three Planck frequencies is
reassuring). Furthermore, we have assumed perfect knowledge of the
relationship between infrared luminosity and SFR, through a simple
proportionality constant (Eq. 11) that assumes particular universal initial
mass function (Chabrier, 2003). This assumption would affect our results and
the direct measurements of $\rho_{\rm SFR}$ in the same manner. The model also
assumes that the SFR at a given time is proportional to the baryon accretion
rate at the same time, whereas the dependence is likely more complicated. Our
results are likely insensitive to most of these assumptions within the
uncertainties reported, but they may become relevant in the presence of more
data.

In terms of systematics, the presence of other sources of extragalactic
contamination in the CIB maps, which would correlate with the galaxy
distribution, could bias the measured values of $\langle b\rho_{\rm
SFR}\rangle$. The most harmful contaminants in this case are thermal Sunyaev-
Zel’dovich and emission from radio point sources. The fact that we recover
compatible results from the three Planck frequencies reassures that any
residual contamination is likely negligible within the range of scales and
frequencies explored here. This may not be the case if a more ambitious
analysis was carried out, using smaller-scale data or other frequencies, and a
more thorough quantification of the impact from possible sources of
contamination would have been necessary. On the galaxy side, Galactic
extinction may cause a non-cosmological modulation in the observed number
counts of galaxies. Since Galactic dust is also a contaminant of the CIB maps,
the cross-correlation between both tracers may be affected by this systematic.
In our case, the effect should be negligible given the conservative large-
scale cut we have imposed, but a more careful study of this cross-
contamination may be required for more sensitive data, or if the large-scale
correlation of both tracers were to be exploited (as would be the case if
targeting ultra-large scale cosmological science cases such as primordial non-
Gaussianity Tucci et al. (2016)). Finally, photometric redshift uncertainties
leading to an incorrect determination of the galaxy redshift distribution may
also impact the recovered constraints on $\langle b\rho_{\rm SFR}\rangle$
(Hang et al., 2021). In our case this is only a potential concern for the DELS
sample. Uncertainties in the sample’s $p(z)$ can arise qualitatively in two
ways: shifts in the mean of the distribution and changes in its width. Given
the broad kernel of the CIB, our constraints should be impervious to moderate
mean shifts. An incorrect determination of the width, however, would affect
the amplitude of the galaxy auto-correlation, leaving the cross-correlation
with the CIB almost untouched (Alonso et al., 2021). This would lead to an
incorrect determination of the galaxy bias, which would then translate into a
biased $\langle b\rho_{\rm SFR}\rangle$. However, by examining the cross-bin
correlations of the DELS sample used here, Hang et al. (2021) were able to
tightly constrain the overlap between bins (which would change with the $p(z)$
width), and therefore this effect should be small in our analysis. A careful
calibration of the redshift distribution would be needed, however, when
attempting this measurement with other photometric samples.

Further progress in the study of the CIB, and in its use as a cosmological
tool, will become possible with the advent of next-generation experiments. On
the CIB side, high-resolution data from experiments such as CCAT-prime (Stacey
et al., 2018) or the Simons Observatory (Ade et al., 2019) will improve the
quality, sensitivity, and frequency coverage of current CIB observations. On
the galaxy side, future deep photometric surveys such as the Rubin Observatory
Legacy Survey of Space and Time (LSST, LSST Science Collaboration et al.,
2009) will enable us to trace the large-scale structure tomographically to
redshifts $z\sim 3$, significantly improving our ability to trace the
evolution of star formation at earlier times. Likewise, spectroscopic surveys
such as DESI (Levi et al., 2019) will provide denser quasar samples, enabling
the same type of constraints with better control over redshift uncertainties.
With a better understanding of the star formation history, it will also be
possible to use the CIB itself as a cosmological tool, allowing us to map the
Universe’s structures on very large scales with a sensitivity higher than will
be achievable by future CMB lensing experiments. However, a large number of
modelling and data analysis challenges need to be tackled before we can fully
exploit the rich information encoded in the CIB.

## Acknowledgements

We would like to thank Jens Chluba, Martin Rey, Aprajita Verma, and Ziang Yan
for useful discussions. BJ is supported by the ENS Paris-Saclay. CGG is
supported by European Research Council Grant No: 693024 and the Beecroft
Trust. DA is supported by the Science and Technology Facilities Council
through an Ernest Rutherford Fellowship, grant reference ST/P004474. JRZ and
NK are supported by an STFC doctoral studentship. NK also acknowledges funding
from the University of Oxford, Department of Physics. We also made extensive
use of computational resources at the University of Oxford Department of
Physics, funded by the John Fell Oxford University Press Research Fund.

We made extensive use of the numpy (Oliphant, 2006; Van Der Walt et al.,
2011), scipy (Virtanen et al., 2020), astropy (Astropy Collaboration et al.,
2013, 2018), healpy (Zonca et al., 2019), GetDist Lewis (2019), and matplotlib
(Hunter, 2007) python packages.

The Legacy Surveys consist of three individual and complementary projects: the
Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David
Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID
#2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy
Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS
together include data obtained, respectively, at the Blanco telescope, Cerro
Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward
Observatory, University of Arizona; and the Mayall telescope, Kitt Peak
National Observatory, NOIRLab. The Legacy Surveys project is honored to be
permitted to conduct astronomical research on Iolkam Duág (Kitt Peak), a
mountain with particular significance to the Tohono O’odham Nation.

NOIRLab is operated by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the National Science
Foundation.

This project used data obtained with the Dark Energy Camera (DECam), which was
constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES
Projects has been provided by the U.S. Department of Energy, the U.S. National
Science Foundation, the Ministry of Science and Education of Spain, the
Science and Technology Facilities Council of the United Kingdom, the Higher
Education Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli
Institute of Cosmological Physics at the University of Chicago, Center for
Cosmology and Astro-Particle Physics at the Ohio State University, the
Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho
de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de
Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de
Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia,
Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the
Collaborating Institutions in the Dark Energy Survey. The Collaborating
Institutions are Argonne National Laboratory, the University of California at
Santa Cruz, the University of Cambridge, Centro de Investigaciones
Energeticas, Medioambientales y Tecnologicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Consortium, the University
of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi
National Accelerator Laboratory, the University of Illinois at Urbana-
Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de
Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig
Maximilians Universitat Munchen and the associated Excellence Cluster
Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, the Ohio State University, the University of Pennsylvania, the
University of Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, and Texas A&M University.

BASS is a key project of the Telescope Access Program (TAP), which has been
funded by the National Astronomical Observatories of China, the Chinese
Academy of Sciences (the Strategic Priority Research Program “The Emergence of
Cosmological Structures” Grant # XDB09000000), and the Special Fund for
Astronomy from the Ministry of Finance. The BASS is also supported by the
External Cooperation Program of Chinese Academy of Sciences (Grant #
114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant #
11433005).

The Legacy Survey team makes use of data products from the Near-Earth Object
Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet
Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by
the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director,
Office of Science, Office of High Energy Physics of the U.S. Department of
Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility under the
same contract; and by the U.S. National Science Foundation, Division of
Astronomical Sciences under Contract No. AST-0950945 to NOAO.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions.

SDSS-IV acknowledges support and resources from the Center for High
Performance Computing at the University of Utah. The SDSS website is
www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

## Data Availability

The code developed for this work as well as the derived data sets produced
(power spectra and covariances) are available upon request. The catalogues and
maps used were made publicly available by the authors, as described in the
text.

## References

  * Ade et al. (2019) Ade P., et al., 2019, J. Cosmology Astropart. Phys., 2019, 056
  * Alonso et al. (2019) Alonso D., Sanchez J., Slosar A., LSST Dark Energy Science Collaboration 2019, MNRAS, 484, 4127
  * Alonso et al. (2021) Alonso D., Bellini E., Hale C., Jarvis M. J., Schwarz D. J., 2021, MNRAS, 502, 876
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Bekhti et al. (2016) Bekhti N. B., et al., 2016, Astronomy & Astrophysics, 594, A116
  * Berlind & Weinberg (2002) Berlind A. A., Weinberg D. H., 2002, The Astrophysical Journal, 575, 587
  * Béthermin et al. (2013) Béthermin M., Wang L., Doré O., Lagache G., Sargent M., Daddi E., Cousin M., Aussel H., 2013, A&A, 557, A66
  * Béthermin et al. (2015) Béthermin M., et al., 2015, A&A, 573, A113
  * Béthermin et al. (2017) Béthermin M., et al., 2017, A&A, 607, A89
  * Blum et al. (2016) Blum R. D., et al., 2016, in American Astronomical Society Meeting Abstracts #228. p. 317.01
  * Cao et al. (2020) Cao Y., Gong Y., Feng C., Cooray A., Cheng G., Chen X., 2020, ApJ, 901, 34
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chen et al. (2016) Chen C.-C., et al., 2016, ApJ, 831, 91
  * Chiang et al. (2019) Chiang Y.-K., Ménard B., Schiminovich D., 2019, ApJ, 877, 150
  * Chiang et al. (2020) Chiang Y.-K., Makiya R., Ménard B., Komatsu E., 2020, ApJ, 902, 56
  * Chisari et al. (2019) Chisari N. E., et al., 2019, ApJS, 242, 2
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Physics reports, 372, 1
  * Darwish et al. (2021a) Darwish O., Sherwin B. D., Sailer N., Schaan E., Ferraro S., 2021a, arXiv e-prints, p. arXiv:2111.00462
  * Darwish et al. (2021b) Darwish O., et al., 2021b, MNRAS, 500, 2250
  * Davies et al. (2016a) Davies L. J. M., et al., 2016a, MNRAS, 461, 458
  * Davies et al. (2016b) Davies L. J., et al., 2016b, Monthly Notices of the Royal Astronomical Society, 461, 458
  * Dey et al. (2019) Dey A., et al., 2019, AJ, 157, 168
  * Dole et al. (2006) Dole H., et al., 2006, A&A, 451, 417
  * Dunkley et al. (2013) Dunkley J., et al., 2013, J. Cosmology Astropart. Phys., 2013, 025
  * Dwek et al. (1998) Dwek E., et al., 1998, ApJ, 508, 106
  * Fakhouri et al. (2010) Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267
  * Fang et al. (2022) Fang X., Eifler T., Schaan E., Huang H.-J., Krause E., Ferraro S., 2022, MNRAS, 509, 5721
  * Flaugher et al. (2015) Flaugher B., et al., 2015, AJ, 150, 150
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * García-García et al. (2019) García-García C., Alonso D., Bellini E., 2019, J. Cosmology Astropart. Phys., 2019, 043
  * García-García et al. (2021) García-García C., Ruiz-Zapatero J., Alonso D., Bellini E., Ferreira P. G., Mueller E.-M., Nicola A., Ruiz-Lapuente P., 2021, J. Cosmology Astropart. Phys., 2021, 030
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Gruppioni et al. (2013) Gruppioni C., et al., 2013, MNRAS, 432, 23
  * Hamimeche & Lewis (2008) Hamimeche S., Lewis A., 2008, Phys. Rev. D, 77, 103013
  * Hang et al. (2021) Hang Q., Alam S., Peacock J. A., Cai Y.-C., 2021, MNRAS, 501, 1481
  * Hivon et al. (2002) Hivon E., Górski K. M., Netterfield C. B., Crill B. P., Prunet S., Hansen F., 2002, ApJ, 567, 2
  * Hou et al. (2021) Hou J., et al., 2021, MNRAS, 500, 1201
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Hurier et al. (2015) Hurier G., Douspis M., Aghanim N., Pointecouteau E., Diego J. M., Macias-Perez J. F., 2015, A&A, 576, A90
  * Hurier et al. (2017) Hurier G., Singh P., Hernández-Monteagudo C., 2017, arXiv e-prints, p. arXiv:1711.10774
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ARA&A, 36, 189
  * Kennicutt & Evans (2012) Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531
  * Kereš et al. (2005) Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
  * Knox et al. (2001) Knox L., Cooray A., Eisenstein D., Haiman Z., 2001, ApJ, 550, 7
  * Koukoufilippas et al. (2020) Koukoufilippas N., Alonso D., Bilicki M., Peacock J. A., 2020, MNRAS, 491, 5464
  * Krolewski & Ferraro (2022) Krolewski A., Ferraro S., 2022, J. Cosmology Astropart. Phys., 2022, 033
  * Krolewski et al. (2021) Krolewski A., Ferraro S., White M., 2021, J. Cosmology Astropart. Phys., 2021, 028
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lenz et al. (2019) Lenz D., Doré O., Lagache G., 2019, ApJ, 883, 75
  * Levi et al. (2019) Levi M., et al., 2019, in Bulletin of the American Astronomical Society. p. 57 (arXiv:1907.10688)
  * Lewis (2019) Lewis A., 2019, arXiv e-prints, p. arXiv:1910.13970
  * Limber (1953) Limber D. N., 1953, ApJ, 117, 134
  * Madau & Dickinson (2014a) Madau P., Dickinson M., 2014a, arXiv preprint arXiv:1403.0007
  * Madau & Dickinson (2014b) Madau P., Dickinson M., 2014b, ARA&A, 52, 415
  * Magnelli et al. (2011) Magnelli Elbaz, D. Chary, R. R. Dickinson, M. Le Borgne, D. Frayer, D. T. Willmer, C. N. A. 2011, A&A, 528, A35
  * Magnelli et al. (2013) Magnelli B., et al., 2013, A&A, 553, A132
  * Maniyar et al. (2018) Maniyar A. S., Béthermin M., Lagache G., 2018, A&A, 614, A39
  * Maniyar et al. (2021) Maniyar A., Béthermin M., Lagache G., 2021, A&A, 645, A40
  * Marchetti et al. (2016) Marchetti L., et al., 2016, MNRAS, 456, 1999
  * McCarthy & Madhavacheril (2021) McCarthy F., Madhavacheril M. S., 2021, Phys. Rev. D, 103, 103515
  * Mead & Verde (2021) Mead A. J., Verde L., 2021, MNRAS, 503, 3095
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Ménard et al. (2013) Ménard B., Scranton R., Schmidt S., Morrison C., Jeong D., Budavari T., Rahman M., 2013, arXiv e-prints, p. arXiv:1303.4722
  * Morrison et al. (2017) Morrison C. B., Hildebrandt H., Schmidt S. J., Baldry I. K., Bilicki M., Choi A., Erben T., Schneider P., 2017, MNRAS, 467, 3576
  * Moster et al. (2018) Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822
  * Neveux et al. (2020) Neveux R., et al., 2020, MNRAS, 499, 210
  * Newman (2008) Newman J. A., 2008, ApJ, 684, 88
  * Nicola et al. (2021) Nicola A., García-García C., Alonso D., Dunkley J., Ferreira P. G., Slosar A., Spergel D. N., 2021, J. Cosmology Astropart. Phys., 2021, 067
  * Oliphant (2006) Oliphant T. E., 2006, A guide to NumPy. Vol. 1, Trelgol Publishing USA
  * Osborne et al. (2014) Osborne S. J., Hanson D., Doré O., 2014, J. Cosmology Astropart. Phys., 2014, 024
  * Pandey et al. (2019) Pandey S., et al., 2019, Phys. Rev. D, 100, 063519
  * Partridge & Peebles (1967) Partridge R. B., Peebles P. J. E., 1967, ApJ, 147, 868
  * Peacock & Smith (2000) Peacock J., Smith R., 2000, Monthly Notices of the Royal Astronomical Society, 318, 1144
  * Planck Collaboration et al. (2011) Planck Collaboration et al., 2011, A&A, 536, A18
  * Planck Collaboration et al. (2014a) Planck Collaboration et al., 2014a, A&A, 571, A18
  * Planck Collaboration et al. (2014b) Planck Collaboration et al., 2014b, A&A, 571, A30
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 596, A109
  * Puget et al. (1996) Puget J. L., Abergel A., Bernard J. P., Boulanger F., Burton W. B., Desert F. X., Hartmann D., 1996, A&A, 308, L5
  * Ross et al. (2020) Ross A. J., et al., 2020, MNRAS, 498, 2354
  * Rybicki & Press (1992) Rybicki G. B., Press W. H., 1992, ApJ, 398, 169
  * Sailer et al. (2021) Sailer N., Schaan E., Ferraro S., Darwish O., Sherwin B., 2021, Phys. Rev. D, 104, 123514
  * Sanders et al. (2003) Sanders D. B., Mazzarella J. M., Kim D.-C., Surace J. A., Soifer B. T., 2003, The Astronomical Journal, 126, 1607
  * Schmittfull & Seljak (2018) Schmittfull M., Seljak U., 2018, Phys. Rev. D, 97, 123540
  * Scranton et al. (2005) Scranton R., et al., 2005, ApJ, 633, 589
  * Seljak (2000) Seljak U., 2000, Monthly Notices of the Royal Astronomical Society, 318, 203
  * Serra et al. (2014) Serra P., Lagache G., Doré O., Pullen A., White M., 2014, A&A, 570, A98
  * Shang et al. (2012) Shang C., Haiman Z., Knox L., Oh S. P., 2012, MNRAS, 421, 2832
  * Silk (2003) Silk J., 2003, MNRAS, 343, 249
  * Silva et al. (2016) Silva D. R., et al., 2016, in American Astronomical Society Meeting Abstracts #228. p. 317.02
  * Stacey et al. (2018) Stacey G. J., et al., 2018, in Marshall H. K., Spyromilio J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 10700, Ground-based and Airborne Telescopes VII. p. 107001M (arXiv:1807.04354), doi:10.1117/12.2314031
  * Takada & Hu (2013) Takada M., Hu W., 2013, Phys. Rev. D, 87, 123504
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Takeuchi et al. (2003) Takeuchi T. T., Yoshikawa K., Ishii T. T., 2003, The Astrophysical Journal, 587, L89
  * Tinker & Wetzel (2010) Tinker J. L., Wetzel A. R., 2010, ApJ, 719, 88
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, The Astrophysical Journal, 688, 709
  * Tinker et al. (2010) Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878
  * Tinsley (1980) Tinsley B. M., 1980, Fundamentals Cosmic Phys., 5, 287
  * Tröster et al. (2022) Tröster T., et al., 2022, A&A, 660, A27
  * Tucci et al. (2016) Tucci M., Desjacques V., Kunz M., 2016, MNRAS, 463, 2046
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science & Engineering, 13, 22
  * Viero et al. (2013) Viero M. P., et al., 2013, ApJ, 772, 77
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2015) Wang L., et al., 2015, MNRAS, 449, 4476
  * Yan et al. (2021) Yan Z., et al., 2021, arXiv e-prints, p. arXiv:2102.07701
  * Yan et al. (2022) Yan Z., van Waerbeke L., Wright A. H., Bilicki M., Gu S., Hildebrandt H., Maniyar A. S., Tröster T., 2022, arXiv e-prints, p. arXiv:2204.01649
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298
  * Zou et al. (2019) Zou H., et al., 2019, ApJS, 245, 4
  * van Engelen et al. (2014) van Engelen A., Bhattacharya S., Sehgal N., Holder G. P., Zahn O., Nagai D., 2014, ApJ, 786, 13

##  Appendix A $\langle b\rho_{\rm SFR}\rangle$ constraints from alternative
analysis choices

Table 3 shows the constraints on $\langle b\rho_{\rm SFR}\rangle$ obtained by
treating each redshift bin independently. The constraints found with a less
conservative scale cut of $k_{\rm max}=0.2\,{\rm Mpc}^{-1}$ are shown in Table
4. Finally, the results found including only one CIB frequency map at a time
are shown in Table 5. The correlation matrix for the measurements presented in
the Table 4 is

${\sf r}=\left(\begin{array}[]{cccccc}1.00&0.11&0.02&0.03&0.00&0.01\\\
0.11&1.00&0.28&0.04&0.01&-0.01\\\ 0.02&0.28&1.00&0.25&0.00&-0.01\\\
0.03&0.04&0.25&1.00&0.01&0.00\\\ 0.00&0.01&0.00&0.01&1.00&0.05\\\
0.01&-0.01&-0.01&0.00&0.05&1.00\end{array}\right).$ (47)

Bin | $b_{g}$ | $\langle b\rho_{\rm SFR}\rangle$ | $\chi^{2}/N_{\rm d.o.f}$ | PTE
---|---|---|---|---
($M_{\odot}{\rm yr}^{-1}{\rm Mpc}^{-3}$)
1 | $1.073\pm 0.044$ | $-0.002\pm 0.008$ | $1.75/10=0.175$ | $0.998$
2 | $1.377\pm 0.026$ | $0.024\pm 0.006$ | $19.98/18=1.109$ | $0.334$
3 | $1.307\pm 0.017$ | $0.041\pm 0.005$ | $34.18/26=1.315$ | $0.131$
4 | $1.733\pm 0.019$ | $0.048\pm 0.005$ | $23.38/30=0.779$ | $0.799$
5 | $2.069\pm 0.146$ | $0.111\pm 0.013$ | $56.23/50=1.125$ | $0.253$
6 | $2.158\pm 0.192$ | $0.264\pm 0.032$ | $61.43/54=1.138$ | $0.227$
Table 3: Summary of the result values for the galaxy bias and the bias weighted star-formation rate within each redshift bin as presented in 2, when all parameters are sampled bin by bin. Bin | $b_{g}$ | $\langle b\rho_{\rm SFR}\rangle$
---|---|---
($M_{\odot}{\rm yr}^{-1}{\rm Mpc}^{-3}$)
1 | $1.125\pm 0.030$ | $0.012\pm 0.005$
2 | $1.366\pm 0.018$ | $0.023\pm 0.004$
3 | $1.307\pm 0.012$ | $0.042\pm 0.004$
4 | $1.739\pm 0.014$ | $0.053\pm 0.004$
5 | $1.970\pm 0.135$ | $0.118\pm 0.013$
6 | $2.321\pm 0.152$ | $0.271\pm 0.024$
$\chi^{2}/N_{\rm d.o.f}=270.1/232=1.16,\hskip 12.0pt{\rm PTE}=0.044$
Table 4: Summary of the result values for the galaxy bias and the bias weighted star-formation rate for the least restrictive value of $k_{\rm max}=0.20\,{\rm Mpc}^{-1}$. Bin | $b_{g}$ | $\langle b\rho_{\rm SFR}\rangle$
---|---|---
($M_{\odot}{\rm yr}^{-1}{\rm Mpc}^{-3}$)
$353$ GHz
1 | $1.071\pm 0.045$ | $0.008\pm 0.025$
2 | $1.382\pm 0.026$ | $0.044\pm 0.017$
3 | $1.312\pm 0.016$ | $0.045\pm 0.014$
4 | $1.736\pm 0.019$ | $0.063\pm 0.012$
5 | $2.064\pm 0.149$ | $0.133\pm 0.020$
6 | $2.244\pm 0.186$ | $0.263\pm 0.033$
$\chi^{2}/N_{\rm d.o.f}=84.6/88=0.96,\hskip 12.0pt{\rm PTE}=0.583$
$545$ GHz
1 | $1.073\pm 0.045$ | $0.000\pm 0.017$
2 | $1.382\pm 0.026$ | $0.041\pm 0.012$
3 | $1.311\pm 0.017$ | $0.036\pm 0.010$
4 | $1.737\pm 0.019$ | $0.054\pm 0.009$
5 | $2.062\pm 0.151$ | $0.128\pm 0.017$
6 | $2.244\pm 0.187$ | $0.254\pm 0.030$
$\chi^{2}/N_{\rm d.o.f}=89.5/88=1.02,\hskip 12.0pt{\rm PTE}=0.440$
$857$ GHz
1 | $1.073\pm 0.045$ | $-0.003\pm 0.011$
2 | $1.382\pm 0.026$ | $0.034\pm 0.008$
3 | $1.310\pm 0.017$ | $0.039\pm 0.007$
4 | $1.736\pm 0.019$ | $0.051\pm 0.007$
5 | $2.075\pm 0.148$ | $0.119\pm 0.014$
6 | $2.246\pm 0.185$ | $0.227\pm 0.030$
$\chi^{2}/N_{\rm d.o.f}=94.7/88=1.08,\hskip 12.0pt{\rm PTE}=0.296$
Table 5: Summary of the result values for the galaxy bias and the bias
weighted star-formation rate within each redshift bin as presented in 2, when
all parameters are sampled frequency per frequency.