# Constraining the physics of star formation from CIB-cosmic shear cross-correlations

Baptiste Jego, David Alonso, Carlos García-García, Jaime Ruiz-Zapatero

► To cite this version:

Baptiste Jego, David Alonso, Carlos García-García, Jaime Ruiz-Zapatero. Constraining the physics of star formation from CIB-cosmic shear cross-correlations. Monthly Notices of the Royal Astronomical Society, 2023, 520 (1), pp.583-598. 10.1093/mnras/stad174. hal-04819770

HAL Id: hal-04819770

https://hal.science/hal-04819770v1

Submitted on 4 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining the physics of star formation from CIB-cosmic shear cross-correlations

Baptiste Jego,<sup>1,2</sup> David Alonso,<sup>2★</sup> Carlos García-García<sup>2</sup> and Jaime Ruiz-Zapatero<sup>1,2</sup>

Accepted 2023 January 8. Received 2022 December 29; in original form 2022 November 2

## ABSTRACT

Understanding the physics of star formation is one of the key problems facing modern astrophysics. The cosmic infrared background (CIB), sourced by the emission from all dusty star-forming galaxies since the epoch of reionization, is a complementary probe to study the star formation history, as well as an important extragalactic foreground for studies of the cosmic microwave background. In this paper, we make high signal-to-noise measurements of the cross-correlation between maps of the CIB from the Planck experiment, and cosmic shear measurements from the Dark Energy Survey and Kilo-Degree Survey. Cosmic shear is a direct tracer of the matter distribution and thus we can use its cross-correlation with the CIB to directly test our understanding of the link between the star formation rate (SFR) density and the matter density. We use our measurements to place constraints on a halo-based model of the SFR that parametrizes the efficiency with which gas is transformed into stars as a function of halo mass and redshift. These constraints are enhanced by using model-independent measurements of the bias-weighted SFR density extracted from the tomographic cross-correlation of galaxies and the CIB. We are able to place constraints on the peak efficiency at low redshifts, $\eta = 0.445^{+0.055}_{-0.11}$, and on the halo mass at which this peak efficiency is achieved today $\log_{10}(M_1/M_{\odot}) = 12.17 \pm 0.25$. Our constraints are in excellent agreement with direct measurements of the SFR density, as well as other CIB-based studies.

**Key words:** galaxies: star formation – large-scale structure of Universe.

## 1 INTRODUCTION

The study of the cosmic microwave background (CMB) has evolved rapidly during the last two decades (Spergel et al. 2003, 2007; Komatsu et al. 2009, 2011; Dunkley et al. 2011; Hinshaw et al. 2013; Sievers et al. 2013; Planck Collaboration XVI 2014a; Planck Collaboration XIII 2016; Calabrese et al. 2017; Planck Collaboration VI 2020; Bianchini et al. 2020; Aiola et al. 2020). Thanks to the swift advance in detector technology, and the construction of new ground-based facilities, we now have at our disposal high-sensitivity, wide-area maps of the radio and infrared sky at arc-minute resolutions. This situation will keep improving rapidly in the next few years with the advent of new facilities, such as the Simons Observatory (Ade et al. 2019), the LiteBIRD satellite (Hazumi et al. 2020), and CMB S4 (Abazajian et al. 2016). However, through this gain in sensitivity, the field has transitioned to a regime where the contribution from various non-cosmological contaminants can no longer be ignored or avoided via masking, and must instead be modelled and incorporated into the cosmological analysis. This is particularly true for total intensity observations, where extra-galactic emission from radio sources, the thermal and kinematic Sunyaev-Zel'dovich (SZ) effects, and the cosmic infrared background (CIB) (Dwek et al. 1998; Puget et al. 1996), dominate the emission on small scales at all frequencies

(Planck Collaboration XVIII 2011; Dunkley et al. 2013). Although the situation is less dire for polarization data, whose constraining power has greatly increased, total intensity maps are still vital, as they can be used to constrain, for example, the epoch of reionization, primordial non-Gaussianity, and the late-time growth of structure through the kinetic SZ effect (Smith & Ferraro 2016; Münchmeyer et al. 2019; Smith et al. 2018).

As a foreground contaminant, the CIB is particularly complex to treat. Its non-universal spectral energy distribution leads to significant decorrelation between observations at different frequencies, which complicates its removal through the simplest multifrequency component separation methods. As an indirect tracer of the large-scale matter fluctuations peaking at redshift $z \sim 2$, it has a markedly non-Gaussian structure that can contaminate the reconstruction of the CMB lensing potential, with which it correlates strongly (Osborne, Hanson & Doré 2014; van Engelen et al. 2014; Sailer et al. 2021; Darwish et al. 2021a). Understanding the physics of the CIB is therefore of vital importance for CMB cosmology, in order to devise more effective component separation techniques, and to incorporate the residual contamination in the model used to obtain cosmological constraints.

In spite of its role as a nuisance in CMB observations, the CIB is a remarkable tool for astrophysics. In cosmology, it can be used as a probe of structure on ultra-large scales (Tucci, Desjacques & Kunz 2016) and to revert the effects of gravitational lensing on maps of the CMB (Smith et al. 2012; Sherwin & Schmittfull 2015). More

<footnote n="1"><sup>1</sup>ENS Paris-Saclay, F-91 190 Gif-sur-Yvette, France</footnote>

<footnote n="2"><sup>2</sup>Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>* E-mail: david.alonso@physics.ox.ac.uk</footnote>
importantly, the CIB contains invaluable information to improve our understanding of the formation and evolution of galaxies. The widely accepted origin of the CIB is the combined infrared emission from dust in star-forming galaxies, heated by the absorption of ultraviolet (UV) light from massive short-lived stars (Partridge & Peebles 1967; Knox et al. 2001). As such, maps of the CIB contain information about the star formation rate (SFR) history from the epoch reionization until today (Dole et al. 2006). Understanding the formation of stars in galaxies of different types at different cosmic epochs is key to understand the formation and evolution of galaxies themselves (Tinsley 1980). Studies of the SFR history have evolved rapidly thanks to the observation of the UV and infrared luminosity function (Gruppioni et al. 2013; Magnelli et al. 2013; Marchetti et al. 2016; Davies et al. 2016). These studies have shown that the star formation rate density (SFRD) grows swiftly from the epoch of reionization, peaking at $z \sim 2$, and then decreasing as the gas fueling it is depleted (Madau & Dickinson 2014).

Although the global picture is qualitatively well understood, the details of the relation between SFR and galaxy properties, or the properties of the halos these galaxies reside in, is far murkier. Studies of the CIB are able to shed some light. First studying the CIB anisotropies through their auto-correlation can constrain the spatial distribution of infrared sources, although the projected nature of the CIB maps makes it difficult to disentangle the contributions from different redshifts (Shang et al. 2012; Viero et al. 2013; Planck Collaboration XXX 2014c; Maniyar, Béthermin & Lagache 2018, 2021). This can be remedied through cross-correlations. Jego et al. (2022) (J22 hereafter) recently showed that the cross-correlation between the CIB and a set of galaxy samples on large scales can be used to make a model-independent tomographic measurement of the bias-weighted SFR density $\langle b\rho_{SFR} \rangle$ which, when combined with direct measurements of the SFRD $\rho_{SFR}$, can shed light on the relation between SFR and halo mass as a function of redshift. The same cross-correlation on small scales is, in principle, sensitive to the relation between SFR and the properties of the target galaxies (Serra et al. 2014; Wang et al. 2015; Chen et al. 2016; Yan et al. 2022). In this regime the signal is sensitive to the fraction of star-forming galaxies in the target sample, which complicates its modelling and interpretation.

In this paper, we will turn instead to correlations between the CIB and tomographic galaxy weak lensing measurements. Since the cosmic shear signal from galaxies at different redshifts directly traces matter inhomogeneities, this cross-correlation is sensitive to the relation between the SFR and matter densities at different times. This will thus allow us to test the validity of different halo-based SFR models and, in general, to improve our understanding of the connection between star-forming galaxies and the underlying dark matter fluctuations. This approach is complementary to the study of the correlation between the CIB and the lensing convergence of the CMB (Planck Collaboration XVIII 2014b; Maniyar et al. 2018; Cao et al. 2020; Darwish et al. 2021b; McCarthy & Madhavacheril 2021) with two advantages. Firstly, cosmic shear data provide a handle on the redshift dependence of the signal via tomography. Secondly, while the CIB is a known contaminant for CMB lensing reconstruction, no such contamination exists for cosmic shear.

This paper is structured as follows. Section 2 presents the theoretical background and the methods used in the analysis. The data sets used are described in Section 3. The measured cross-correlations are presented and analysed in Section 4, where we also present the associate constraints on star formation models. We summarize and discuss our results in Section 5.

# 2 METHODS

## 2.1 Theory

Our theory prediction will follow the formalism described in (Moster, Naab & White 2018; García-García et al. 2021; Jego et al. 2022).

### 2.1.1 Angular power spectra and the halo model

The cosmic shear signal and CIB anisotropies can both be described in as projected tracers (i.e. sky maps) $u(\hat{\theta})$ of a three-dimensional field $U(\mathbf{x}, z)$ through a radial kernel $q_u(\chi)$:

$$u(\hat{\theta}) = \int d\chi q_u(\chi) U(\chi, \hat{\theta}, z), \quad (1)$$

where $z$ is the redshift corresponding to the comoving distance $\chi$. The angular power spectrum of two such quantities, $u$ and $v$ is then given by

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

where $P_{UV}(k, z)$ is the power spectrum of the corresponding three-dimensional (3D) quantities, and we have made use of the Limber approximation (Limber 1953), which is appropriate for the tracers analysed here. The multiplicative factors $f_{\ell}^{u/v}$ account for the potential angular derivatives relating the two-dimensional and 3D quantities.

We will model $P_{UV}(k, z)$ making use of the halo model (Seljak 2000; Peacock & Smith 2000; Cooray & Sheth 2002). In this formalism

$$P_{UV}(k, z) = P_{UV}^{2h}(k, z) + P_{UV}^{1h}(k, z), \quad (3)$$

where the 1-halo and 2-halo contributions are given by

$$P_{UV}^{1h}(k) \equiv \int dM n(M) \langle U(k, M) V(k, M) \rangle, \quad (4)$$

$$P_{UV}^{2h}(k) \equiv \langle bU(k) \rangle \langle bV(k) \rangle P_{\text{lin}}(k), \quad (5)$$

$$\langle bU(k) \rangle \equiv \int dM n(M) b_h(M) \langle U(k, M) \rangle. \quad (6)$$

Here, $n(M)$ and $b_h(M)$ are the halo mass function and the halo bias respectively for halos of mass $M$, and $P_{\text{lin}}(k, z)$ is the linear matter power spectrum. $\langle U(k, M) \rangle$ is the Fourier transform of the mean halo profile of quantity $U$, i.e.

$$\langle U(k, M) \rangle = 4\pi \int_0^\infty dr r^2 \langle U(r, M) \rangle \frac{\sin kr}{kr}, \quad (7)$$

where $\langle U(r, M) \rangle$ is the mean value of $U$ at a distance $r$ from the center of a halo of mass $M$. Likewise, $\langle U(k, M)V(k, M) \rangle$ is the two-point cumulant of the two Fourier-space profiles.

A prediction for the angular cross-power spectrum between cosmic shear (labelled as $\gamma$ here) and CIB anisotropies at frequency $\nu$ thus requires the radial kernels of both probes, and a model for the statistics (scale-dependent mean and covariance) of the 3D quantities associated with them in halos of different masses. We describe these ingredients in the next two sections.

Although the halo model is able to describe the power spectrum on either large- or small-scales, where either the 2-halo or the 1-halo terms dominate, it is inaccurate at the ~10–20 per cent level on intermediate scales. The most likely cause for this is the oversimplified treatment of halo biasing used in the vanilla version of the halo model used here (Mead et al. 2021; Mead & Verde 2021).
We account for this following (Mead 2017; Koukoufilippas et al. 2020). We multiply the halo model prediction for $P_{\gamma\gamma}(k, z)$ in the range $0.05 \text{ Mpc}^{-1} < k < 2 \text{ Mpc}^{-1}$ by the factor

$$R(k, z) = \frac{P_{\text{halofit}}(k, z)}{P_{\text{HM}}(k, z)}, \qquad (8)$$

where $P_{\text{HM}}(k, z)$ is the halo model prediction for the matter-matter power spectrum, and $P_{\text{halofit}}(k, z)$ is the fit to the same quantity using the Halofit (Smith et al. 2003) parametrization of Takahashi et al. (2012).

### 2.1.2 Cosmic shear and matter fluctuations

Weak gravitational lensing distorts the shapes of background galaxies, correlating their ellipticities. This effect, known as ‘cosmic shear’ is quantified through a spin-2 projected field $\gamma$. At leading order, weak lensing contributes only to the parity-even ‘E-mode’ component of the field. This is a projected tracer of the matter overdensity $\Delta_m(\mathbf{x}, z)$ with radial kernel (Bartelmann & Schneider 2001)

$$q_{\gamma}(\chi) = \frac{3}{2} H_0^2 \Omega_m (1+z) \chi \int_z^\infty dz' p(z') \frac{\chi(z') - \chi}{\chi(z')}, \qquad (9)$$

where $H_0$ is the current value of the expansion rate, $\Omega_m$ is the fractional energy density of non-relativistic matter, and $p(z)$ is the redshift distribution of the source galaxies. We use natural units where the speed of light is $c = 1$.

The $\ell$-dependent prefactor, due to the relation between $\gamma$ and the angular Hessian of the Newtonian gravitational potential, is

$$f_{\ell}^{\gamma} = \sqrt{\frac{(\ell+2)!}{(\ell-2)!} \frac{1}{(\ell+1/2)^2}}, \qquad (10)$$

which is negligibly different from 1 on the scales used here.

The associated halo profile is simply the matter density profile normalized by the mean background matter density $\bar{\rho}_m$. For this we use the truncated Navarro–Frenk–White parametrization (Navarro, Frenk & White 1996). In this case

$$\langle u_{\gamma}(k, M) \rangle = \frac{M}{\bar{\rho}_m} u_{\text{NFW}}(k, M), \qquad (11)$$

where

$$u_{\text{NFW}}(k) = \left[ \ln(1+c) - \frac{c}{1+c} \right]^{-1} \left\{ \sin x [\text{Si}((1+c)x) - \text{Si}(x)] + \cos x [\text{Ci}((1+c)x) - \text{Ci}(x)] - \frac{\sin(cx)}{(1+c)x} \right\}. \qquad (12)$$

Here $x \equiv kR_{\Delta}/c(M)$, $R_{\Delta}(M)$ is the halo virial radius, $c(M)$ is the concentration-mass relation, and Si/Ci are the sine and cosine integrals.

The intrinsic alignment (IA) of the shapes of galaxies, caused by local gravitational tidal forces, is an important source of contamination for cosmic shear (Brown et al. 2002). Unlike the cosmic shear signal, IAs are a local effect and correlate with the structure at redshift of the source galaxies. Within the linear non-linear alignment model (LNLA, Hirata & Seljak (2004)), the effect of IAs can be taken into account by adding a local contribution to the cosmic shear kernel of the form:

$$q_{I}(\chi) = -A_{\text{IA}}(z) H(z) p(z), \qquad (13)$$

where $H(z)$ is the expansion rate at redshift $z$, and $A_{\text{IA}}$ is a linear amplitude parameter that parametrizes the strength with which local

tidal forces modify galaxy ellipticities. A common parametrization of this amplitude is (e.g. Troxel et al. 2018)

$$A_{\text{IA}}(z) = A_0 \left( \frac{1+z}{1+z_0} \right)^{\lambda} \frac{0.0139 \Omega_m}{D(z)}, \qquad (14)$$

where $D(z)$ is the linear growth factor, $z_0$ is a fixed pivot redshift, and $A_0$ and $\lambda$ are free parameters.

We must also propagate uncertainties due to the poor knowledge of the source redshift distribution $p(z)$ of the shear samples used here. We do so by marginalizing over a nuisance parameter $\Delta z$ in each redshift bin describing a shift in the mean of the redshift distribution:

$$p(z) \rightarrow p(z + \Delta z). \qquad (15)$$

This has been shown to encapsulate the main effect of redshift uncertainties on the weak lensing kernel (Tessore & Harrison 2020; Cordero et al. 2022).

### 2.1.3 Cosmic infrared background and star formation

To model the CIB signal we follow the prescription of Moster et al. (2018), Maniyar et al. (2018), Maniyar, Béthermin & Lagache (2021), Jego et al. (2022) (the reader is referred to these papers for further details to avoid repetition).

The CIB specific intensity at a given observed frequency $\nu$, $I_{\nu}(\hat{\theta})$, is a projected tracer of the SFRD, $\rho_{\text{SFR}}$:

$$I_{\nu}(\hat{\theta}) = \int d\chi q_{\nu}(\chi) \rho_{\text{SFR}}(\chi\hat{\theta}, z), \qquad (16)$$

with a radial kernel given by

$$q_{\nu}(\chi) = \frac{\chi^2 S_{\nu}^{\text{eff}}(z)}{K}. \qquad (17)$$

Here, $K = 10^{-10} \text{ M}_\odot \text{ yr}^{-1} \text{ L}_\odot^{-1}$ is the calibration constant relating the far infrared luminosity $L_{\text{IR}}$ and star formation rate ($L_{\text{IR}} = \text{SFR}/K$ Kennicutt 1998; Kennicutt & Evans 2012) for a Chabrier initial mass function (Chabrier 2003), and $S_{\nu}^{\text{eff}}(z) \equiv S_{\nu}(z)/L_{\text{IR}}$ is the normalized mean spectral energy distribution of sources as a function of redshift.<fnref n="1" /> The $\ell$-dependent factor is simply $f_{\ell}^{\nu} = 1$.

To model the SFRD profile, we add the contribution from central and satellite galaxies. The resulting mean Fourier-space profile is

$$\langle u_{\text{SFR}}(k, M) \rangle = \text{SFR}_c(M, z) + \text{SFR}_s(M, z) u_s(k, M), \qquad (18)$$

where $SFR_{c/s}(M, z)$ is the star formation from centrals and satellites in a parent halo of mass $M$ at redshift $z$, and $u_s(k, M)$ is the mean distribution of satellites in Fourier space [normalized to $u_s(k \rightarrow 0, M) = 1$]. For simplicity, we assume that satellites follow the dark matter distribution, and use $u_s(k, M) = u_{\text{NFW}}(k, M)$. Note that this is an assumption that can be explicitly tested through the CIB-cosmic shear cross-correlation on small, halo-sized scales, although we do have not attempted to do so here.

As in Maniyar et al. (2018), Maniyar et al. (2021), Moster et al. (2018), we parametrize the SFR of the central galaxy in terms of the efficiency $\eta$ to convert infalling gas into stars:

$$\text{SFR}_c(M, z) = \eta(M, z) \text{BAR}(M, z), \qquad (19)$$

<footnote n="1">¹We use the estimates of Béthermin et al. (2013), Béthermin et al. (2015), Béthermin et al. (2017) for the Planck channels, made available in https://github.com/abhimaniyar/halomodel.cib.tsz.cibxtsz.</footnote>
where BAR is the baryonic accretion rate, which we model after Fakhouri, Ma & Boylan-Kolchin (2010) as

$$ \text{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}, \quad (20) $$

with $\dot{M}_0 = 46.1 M_{\odot} \text{ yr}^{-1}$.

To calculate the satellite contribution, we will assume that all subhalos in a parent halo of mass M contain a satellite:

$$ \text{SFR}_s(M, z) = \int_{M_{\text{min}}}^{M} dM_{\text{sub}} \frac{dN}{dM_{\text{sub}}} \text{SFR}_{\text{sat}}(M_{\text{sub}}, M, z), \quad (21) $$

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

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

This guarantees that the SFR in any satellite never exceeds that of the central galaxy (see Maniyar et al. 2021).

Most parametrizations of $\eta(M, z)$ assume that the efficiency peaks at a particular halo mass $M_{\text{max}} \sim 10^{12-13} M_{\odot}$. This behaviour is physically motivated (Silk 2003; Kereš et al. 2005): at lower masses, the weaker gravitational potential and the impact of supernova feedback deplete the galaxy of gas. In turn, at higher masses the increased gas cooling time and the impact of AGN feedback have a similar effect. Here, we will choose the parametrization of Moster et al. (2018) (M18 hereafter), which assumes

$$ \eta(M, z) = \frac{2\eta_*}{(M_1/M)^\beta + (M/M_1)^\gamma}. \quad (22) $$

All parameters ($\eta_*$, $m \equiv \log_{10} M_1/M_{\odot}$, $\beta$, and $\gamma$) are allowed to vary with redshift as:

$$ x(z) = x_0 + x_z \frac{z}{1+z}. \quad (23) $$

As in J22, we will fix $\beta_0$, $\beta_z$, $\gamma_0$, and $\gamma_z$ to the best-fit values found by M18, and we will only consider four free parameters: $\eta_0$, $\eta_z$, $m_0$, and $m_z$. These parametrize the value and time dependence of the peak efficiency and the associated mass. Note that other models have been proposed in the literature. For instance, Maniyar et al. (2021), proposed using a log-normal form for the mass dependence of $\eta(M, z)$, with a constant peak efficiency and mass, but a redshift-dependent high-mass tail. We choose to follow the parametrization of M18 here, since it arguably allows for a more flexible, but simple, redshift evolution of the efficiency curve.

Given a halo-based model for the SFR, the mean star formation rate density can be calculated by simply integrating over the halo mass function:

$$ \rho_{\text{SFR}}(z) = \int dM n(M, z) \text{SFR}(M, z). \quad (24) $$

Likewise, the bias-weighted version of the same quantity is

$$ \langle b \rho_{\text{SFR}} \rangle = \int dM n(M, z) b_h(M, z) \text{SFR}(M, z). \quad (25) $$

## 2.2 Angular power spectra and covariances

We use the MASTER algorithm (Hivon et al. 2002) as implemented in NaMaster (Alonso et al. 2019) to estimate all power spectra used in this analysis. The method is based on estimating the statistical coupling between different power spectrum multipoles, caused by the survey mask, (the so-called mode-coupling matrix) using analytical

<figure><img src="image_2.png" /><figcaption>Figure 1. Bandpower window functions for the cross-correlation between the third KiDS redshift bin and the 545 GHz CIB map.</figcaption></figure>

approach enabled by the orthogonality of the Wigner-3j symbols. A thorough description of the method can be found in Alonso et al. (2019). Here, we will only describe some of the details specific to the use of cosmic shear data.

We follow the prescription of Nicola et al. (2021). The cosmic shear field $\gamma$ is only sampled at the positions of the observed galaxies. As such, the natural survey mask that optimizes the pixel-level signal-to-noise is given by the sum of shape measurement weights for all galaxies in each pixel. This leads to a complex mask with substantial structure, which in turn can cause a non-negligible statistical coupling between distant multipoles. When comparing the estimated power spectra with theoretical predictions we therefore convolve the predictions with the bandpower window functions, which fully describe this coupling:

$$ C_b = \sum_{\ell} \mathcal{F}_b^{\ell} C_{\ell}, \quad (26) $$

where $C_\ell$ is the per-$\ell$ theoretical power spectrum, $\mathcal{F}_b^\ell$ is the window function for the b-th bandpower, and $C_b$ is the corresponding prediction for that bandpower. To illustrate this, Fig. 1 shows the bandpower windows for the cross-correlation between the third KiDS bin and the 545 GHz CIB map.

We will only include angular scales $100 \le \ell \le 1500$ in our analysis. The large-scale cut is motivated by the loss of power on large scales in the CIB maps due to Galactic dust removal, while the small-scale cut ensures that our analysis is not significantly affected by modelling uncertainties on the impact of baryonic feedback processes on the matter power spectrum. Using the same $\ell$ bins as (García-García et al. 2021), this leaves 20 bandpower measurements in for each cross correlation.

We use the analytical approximation of García-García, Alonso & Bellini (2019), Nicola et al. (2021) (the so-called ‘Narrow-Kernel Approximation’) to estimate the Gaussian covariance matrix of the data, neglecting all non-Gaussian contributions. The approximation relies on an estimate of the angular power spectra of the different tracers involved. For this, we use the pseudo-$C_\ell$ estimate from the data for all pairs of maps, normalized by the mean of the product of their respective masks. This has the advantage of not relying on a particular astrophysical model for the signal, or a precise instrument model for the noise. We validate this estimate by recomputing the covariance for a subset of the spectra via jackknife resampling. As in J22, we find that the analytical approximation works well at the per-cent level, although a small correction of a few per cent must be applied to covariance matrix elements involving two different CIB
<table><caption>**Table 1.** Prior distributions for the nuisance parameters entering our analysis for each tracer. $U(a, b)$ and $\mathcal{N}(\mu, C)$ describe a uniform distribution with boundaries $(a, b)$ and a multivariate Gaussian distribution with mean $\mu$ and covariance $C$, respectively.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td>$m_0$</td><td>$U(9, 14)$</td><td>$m_z$</td><td>$U(-6, 6)$</td></tr><tr><td>$\eta_0$</td><td>$U(0, 1)$</td><td>$\eta_z$</td><td>$U(0, 1)$</td></tr><tr><td>$A_0^{\text{DES}}$</td><td>$\mathcal{N}(1, 0.5)$</td><td>$A_0^{\text{KiDS}}$</td><td>$\mathcal{N}(1, 0.5)$</td></tr><tr><td>$\Delta z_{\text{DES}}$ (x4)</td><td>$\mathcal{N}(\mu, C)$</td><td>$\Delta z_{\text{KiDS}}$ (x5)</td><td>$\mathcal{N}(\mu, C)$</td></tr></tbody></table>

maps. This is important in a multifrequency analysis given the very high correlation between the different frequencies. Since, as we will describe in Section 3.2.2, we will instead combine all frequencies at the map level for this analysis, this modification is less critical than in J22.

## 2.3 Likelihood

We will make use of a Gaussian likelihood,

$$ -2 \log p(\mathbf{d}|\Theta) = (\mathbf{d} - \mathbf{m}(\Theta))^T C_d^{-1} (\mathbf{d} - \mathbf{m}(\Theta)) + K, \quad (27) $$

to constrain the free parameters of the star-formation model. In this expression $\mathbf{d}$ is the data vector, containing all cross-correlations between tomographic cosmic shear bins and CIB maps, $C_d$ is the covariance of $\mathbf{d}$, $\mathbf{m}$ is the theoretical model for $\mathbf{d}$, and $\Theta$ denotes the free parameters of the model. The posterior distribution is then given by the product of this likelihood and the parameter priors.

Our fiducial model contains 4 parameters describing the halo model for the star formation history, described in Section 2.1.3: {$m_0, m_z, \eta_0, \eta_z$}. We also marginalize over two intrinsic alignment parameters $A_0^{\text{DES/KiDS}}$, describing the amplitude of IAs in the DES and KiDS samples. For both parameters we use a Gaussian prior centered at $A_0 = 1$ with standard deviation $\sigma(A_0) = 0.5$. This choice encompasses the values favoured in the cosmological analyses of (Troxel et al. 2018; Heymans et al. 2021). We fix the IA evolution parameter to $\lambda = 0$ in both cases. Finally, we marginalize over one redshift shift parameter $\Delta z^i$ for each redshift bin (five for KiDS, four for DES). For these, we use the Gaussian priors described in (Troxel et al. 2018; Heymans et al. 2021). Shape measurement uncertainties were propagated by marginalizing over one multiplicative bias parameter per redshift bin. This was done analytically, following the procedure of Joachimi et al. (2021) with the Gaussian priors of (Troxel et al. 2018; Heymans et al. 2021). In the most general case, our model therefore has 15 free parameters. The priors used for all of them are summarized in Table 1.

We sample this posterior distribution using the Metropolis-Hastings Markov-Chain Monte-Carlo (MCMC) algorithm implemented in `Cobaya` (Torrado & Lewis 2021, 2019). All cosmological theory predictions were computed using the Core Cosmology Library (CCL, Chisari et al. (2019)<fnref n="2" />). 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)$. We use the halo mass function and halo bias parametrization of Tinker et al. (2010), with halo masses defined for a spherical overdensity $\Delta = 200$ with respect to the critical density.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Radial kernels of the five KiDS (red) and four DES (blue) redshift bins used in this analyses. The vertical dashed line marks the redshift at which we evaluate the mean CIB spectra to define the coadding weights.</figcaption></figure>

# 3 DATA

## 3.1 Cosmic shear

We make use of two weak lensing catalogs released by the Dark Energy Survey and Kilo-Degree Survey collaborations. The methods used to process these catalogs are described in detail in García-García et al. (2021), and mostly follow the same procedures used in the official cosmological analyses of both collaborations (Zuntz et al. 2018; Troxel et al. 2018; Giblin et al. 2021; Asgari et al. 2021). We outline only the main features of these data, and refer the reader to García-García et al. (2021), as well as the original papers for further details.

The DES catalog corresponds to the official cosmic shear sample from the first-year (Y1) data release. We make use of the Metacalibration catalog, which uses the shape-measurement algorithm Sheldon & Huff (2017) of the same name to determine individual galaxy ellipticities. The catalog is split into four different tomographic redshift bins spanning the range $z \lesssim 1.5$. We make shear maps by averaging the two galaxy ellipticities in each pixel, corrected by the mean response tensor trace. All galaxies are given equal weights, and a mean ellipticity is subtracted from the maps in each redshift bin, following Troxel et al. (2018). We use the official redshift distributions released by DES, marginalizing over shifts in the mean redshift as described in Section 2.1.2. The catalog covers an area of ~1300 deg<fnref n="2" />, of which approximately 1000 deg<fnref n="2" /> overlap with the CIB maps used here.

We also use the cosmic shear catalog publicly released by the KiDS collaboration as part of the fourth data release (DR4), the so-called KiDS-1000 sample (Giblin et al. 2021). The sample is divided into the same five redshift bins used by the KiDS collaboration, and we use the official redshift distributions released by KiDS to carry out our analysis (again marginalizing over mean shifts). Galaxy shear is estimated using `lensfit` (Miller et al. 2013), and we generate shear maps taking the corresponding shape measurement weights into account. We correct for a mean multiplicative bias and subtract the residual mean ellipticity in each redshift bin. The catalog covers an area of ~1000 deg<fnref n="2" />, of which ~500 deg<fnref n="2" /> overlap with the CIB footprint.

Fig. 2 shows the weak lensing kernels, defined in equation (9), corresponding to each of the nine redshift bins (four DES bins in solid blue, five KiDS bins in dashed red). The four DES bins lead to lensing kernels that are remarkable similar to the four higher-

<footnote n="2">²The source code can be found at https://github.com/LSSTDESC/CCL</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. Sky footprint of the CIB and weak lensing surveys used in this analysis in Galactic coordinates.</figcaption></figure>

redshift KiDS bins. The weak lensing signal in the first KiDS bin is rather small and, as we will see, contributes negligibly to our cross-correlation analysis. The last two redshift bins, which dominate the cross-correlation signal, peak at around $z \sim 0.5$ (marked with a dashed black vertical line in the figure), and thus our measurements are mostly sensitive to the star formation history in this range of redshifts. The sky footprints of both samples are shown in Fig. 3, together with that of the CIB maps described in the next section.

## 3.2 Cosmic infrared background maps

### 3.2.1 Multifrequency maps

For our analysis, we use the CIB maps constructed by Lenz, Doré & Lagache (2019) from the Planck 353, 545 and 857 GHz temperature maps, corrected for contamination from Galactic dust using neutral hydrogen data from the HI4PI survey (HI4PI Collaboration et al. 2016). As in J22, we use the 20 per cent sky masks provided with these data, with an HI column density threshold of $N_{HI} > 2.5 \times 10^{20} \text{ cm}^{-2}$, apodized with a 15 arcmin FWHM Gaussian kernel. All power spectra were corrected for the effective beam of each of these maps, which is also provided in the public data release.

To avoid the loss of power on large scales reported by Lenz et al. (2019), due to the local removal of Galactic dust in overlapping sky patches, we exclude multipoles $l < 100$ from the analysis.

### 3.2.2 A coadded cosmic infrared background map

As shown in Yan et al. (2022), Jego et al. (2022), the cross-correlation signal of the three different frequency maps is highly (90–95 per cent) correlated. This has the advantage that the analysis can be carried out on each frequency map independently, and the results can then be compared to validate the internal consistency of the model used to connect the CIB to the SFR history. The disadvantage, however, is that care must be taken when combining power spectra from all frequencies, since a slight mis-estimation of the cross-covariance terms involving different frequency channels can lead to artificially high or low $\chi^2$ values and, in the worst case,

to a bias in the final parameter constraints. In J22 we showed that the analytical covariance matrix had to be corrected by a few per cent in these cross-terms in order to match the jackknife estimate. Such a small correction is normally irrelevant; however, in this case it can change the $\chi^2$ values by up to a factor 2, due to the tight correlation between frequencies. Since the current analysis covers a wider range of scales than that of J22, which might require a more careful modelling of this correction, we have opted to avoid it altogether.

To do so, we construct a single co-added map as an optimal linear combination of the three frequency maps. To determine the linear coefficients of this combination we model the three frequency maps as:

$$ \mathbf{m}(\hat{\theta}) = s c(\hat{\theta}) + \mathbf{n}(\hat{\theta}), \qquad (28) $$

where $\mathbf{m}(\hat{\theta})$ is a three-element vector containing the three frequency maps at sky position $\hat{\theta}$, $s$ is a vector containing the effective CIB spectrum, $c(\hat{\theta})$ is the coadded map we are trying to reconstruct, and $\mathbf{n}(\hat{\theta})$ is the contribution from instrumental noise. Assuming Gaussian noise, an optimal estimator for $c$ can then be found via least-squares minimization to be

$$ \hat{c}(\hat{\theta}) = \mathbf{w}^T \mathbf{m}(\hat{\theta}), \qquad \mathbf{w}^T = \frac{\mathbf{s}^T \mathbf{N}^{-1}}{\mathbf{s}^T \mathbf{N}^{-1} \mathbf{s}}, $$

where $N$ is the noise covariance matrix, and $\mathbf{w}$ are the linear weights we were seeking.

One must note that equation (28) is only correct in the limit where the CIB frequency maps are 100 per cent correlated at the signal level. Within the model used here (equation 16), this would only be strictly valid if the frequency and redshift dependence of the effective infrared spectra [$S_v^{\text{eff}}(z)$ in equation 17] were factorizable. Nevertheless, since correlation between the three maps is very high, assuming a perfect correlation in order to find $\mathbf{w}$ will only lead to a small loss of sensitivity when using the resulting coadded map instead of all the frequency maps. Since, as we have shown, the weak lensing kernel peaks in the range $z \sim 0.5$ for the cosmic shear maps, we construct the vector of spectra $\mathbf{s}$ by evaluating $S_v^{\text{eff}}$ at a pivot redshift $z_0 = 0.5$. In detail, in order to preserve the units of the
coadded CIB map (MJy/sr) we construct the elements of s as:

$$s_v = S_{\nu}^{\text{eff}}(z_0)/S_{857}^{\text{eff}}(z_0), \qquad (29)$$

where $S_{857}^{\text{eff}}$ is the spectrum in the 857 GHz channel (our highest signal-to-noise map). The coadded map will therefore have an amplitude comparable to that of the 857 map.

To estimate the noise covariance matrix N we first compute the noise power spectrum of the frequency maps $N_{\ell}^{\nu\nu'}$ (including all auto- and cross-correlations) using the difference between half-mission maps to null out the signal. Each entry of N is then calculated by averaging $(2\ell + 1)N_{\ell}^{\nu\nu'}$ in the range $100 \le \ell \le 1500$ used here. The resulting set of coadd weights is:

$$(w_{353}, w_{545}, w_{857}) = (0.06, 0.22, 0.92). \qquad (30)$$

We found that varying the pivot redshift in the range $0.2 \le z_0 \le 1$, as well as the range of multipoles used to determine N did not change the results significantly.

Finally, since the different frequency maps have different masks and effective beams, in practice we coadd the frequency maps at the level of the power spectra. The resulting cross-power spectrum of cosmic shear with the coadded map, and its covariance become:

$$C_{\ell}^{\gamma, \text{co}} = \sum_{\nu} w_{\nu} C_{\ell}^{\gamma, \nu}, \qquad (31)$$

$$\text{Cov}(C_{\ell}^{\gamma, \text{co}}, C_{\ell'}^{\gamma', \text{co}}) = \sum_{\nu\nu'} w_{\nu} w_{\nu'} \text{Cov}(C_{\ell}^{\gamma, \nu}, C_{\ell'}^{\gamma', \nu}). \qquad (32)$$

Since the weights are scale-independent, this is equivalent to a linear combination at the map level.

# 4 RESULTS

## 4.1 Power spectra, validation, and goodness of fit

The data vector that will form the basis of our fiducial analysis is the set of nine cross-correlations between the coadded CIB map described in Section 3.2.2 and the cosmic shear E-mode in the four DES and five KiDS tomographic redshift bins. These measurements and their statistical uncertainties are shown in Fig. 4. We obtain clear detections of the cross-correlation especially at the higher redshifts, where the amplitude of both the cosmic shear and CIB signals grows. As a first model-independent estimate of the total signal-to-noise ratio (S/N) of this cross-correlation, we find

$$S/N = \sqrt{\chi_0^2 - N_d} = 20.3, \qquad (33)$$

where $\chi_0^2$ is the $\chi^2$ statistic for a null hypothesis (i.e. $\chi_0^2 \equiv d^T C_d^{-1} d$ in the notation of Section 2.3), and $N_d = 180$ is the size of the data vector. This is, therefore one of the highest-significance detections of this cross-correlation (see also Tröster et al. 2022).

The orange solid lines in Fig. 4 show the best-fit theory prediction for the SFR model of M18 (see Section 2.1.3), found by maximizing the posterior distribution of the full set of nine power spectra. The dashed blue lines, in turn, show the best-fit model obtained by fitting the power spectra corresponding to DES and KiDS separately. The similarity between both best-fit models in all cases, compared with the size of the error bars, therefore implies that the measurements made with both surveys are consistent (we will show this over the full parameter space in the next section). The overall best-fit model has a $\chi^2 = 187.6$, which corresponds to an acceptable probability-to-exceed (PTE) of $p = 0.28$ assuming $N_d - 4 = 176$ degrees of

freedom.<fnref n="3" /> Having a best-fit model at hand allows us to compute the detection significance in a different way, as

$$S/N = \sqrt{\chi_0^2 - \chi_{\text{BF}}^2} = 20.1, \qquad (34)$$

where $\chi_{\text{BF}}^2$ is the $\chi^2$ of the best-fit model, and $\chi_0^2$ was defined after equation (33). The different panels in Fig. 4 include the signal-to-noise of each cross-correlation calculated in this manner. We detect the cross-correlation above $4\sigma$ significance in all the DES redshift bins and in the 3 highest redshift KiDS bins, whereas the two lowest redshift KiDS bin contribute negligibly to the signal. The difference in the detectability of the signal in the first DES bin and the second KiDS bin (which have similar kernels, as shown in Fig. 2) is due to the higher area overlap of the DES footprint.

Since weak lensing should lead to an undetectable B-mode signal (at least given current sensitivities), ensuring that the measured B-mode power spectra are compatible with zero is a good test for potential systematics in the data. As shown in Appendix A, we find that the total B-mode signal is compatible with zero, and that this is also the case for most of the individual power spectra. However, we find that the cross-correlation of the second KiDS bin with the coadded CIB map yields an unacceptably low $\chi^2$ probability ($p = 0.001$) for a null signal. We find this to be the case for the cross-correlation of this bin with the individual CIB frequency maps as well. This may be a statistical fluke, or a sign of systematics in this sample, which was also flagged specifically in the official cosmic shear analysis by the KiDS collaboration (Asgari et al. 2021). Since this particular bin does not contribute significantly to the model constraints obtained here, we have included it in our analysis, but verified that our results are unaffected by this choice (see Section 4.3).

Fig. 5 shows the correlation matrix of the power spectra shown in Fig. 4. There are significant off-diagonal correlations (at the level of 30–40 per cent) between different redshift bins in the same survey, which extend to almost all pairs of bins due to the cumulative nature of weak lensing. The covariance between the DES and KiDS measurements is zero since their footprints do not overlap.

We have also studied the cross-correlations with the individual CIB frequency maps. The 12 cross-correlations between the three frequency maps and the four DES redshift bins are shown in Fig. 6, with qualitatively equivalent results found for KiDS. As before, the solid orange lines show the best-fit model found using all the cross-correlations with the coadded CIB map, and the dotted blue line shows the best-fits using only cross-correlations with DES. In turn, the red dashed lines show the best-fit prediction using only the cross-correlations with one of the three frequency maps. The fiducial best-fit model is largely consistent with the predictions from the 545 and 857 GHz channels. It, however, consistently under-predicts the 353 GHz power spectra in both DES and KiDS. A similar trend, where the amplitude of the CIB signal is somewhat consistently higher in the 353 GHz channel was also found in cross-correlation with galaxies by J22. This could be due to both systematics in the data (e.g. contamination from CIB or other extragalactic foregrounds), or a mis-modeling of the CIB spectrum in this channel. Nevertheless, as we will see in Section 4.3, the 353 GHz data lead to constraints on the model parameters that are largely compatible with those of the other two frequencies, and with our fiducial constraints. Thus, and since the coadded CIB map is dominated by the 857 and 545 GHz channels,

<footnote n="3">³We subtract only the 4 SFR parameters from the total number of datapoints when computing the number of degrees of freedom, since all other model parameters have tight priors. In any case, subtracting all model parameters yields a lower, but still acceptable probability $p = 0.12$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top panels: cross-power spectra between the coadded CIB map and the four DES tomographic samples. Bottom panels: as above for the five KiDS tomographic bins. The orange solid lines show the best-fit model found by fitting all power spectra simultaneously, while blue dotted lines show the best-fit model for each survey (DES or KiDS) independently. We find consistent results between both surveys.</figcaption></figure>

we have kept the 353 GHz data in our fiducial analysis. Furthermore, we find that the best-fit model from the coadded CIB map is a good fit to the full multifrequency data vector, with a $\chi^2 = 547.6$ for 536 degrees of freedom ($p = 0.35$). Calculating the significance of the detection for the multifrequency data, as described above in the case of the coadded power spectra, yields

$$S/N = \sqrt{951.3 - 547.6} = 20.1. \qquad (35)$$

This is in striking agreement with the detection significance using only the coadded CIB map, which demonstrates that virtually no information is lost by coadding the three CIB maps, due to the tight correlation between them. This has the advantage of simplifying the analysis significantly, since each likelihood evaluation is three times faster and, more importantly, we are not sensitive to the impact on the final parameter constraints of a potential mis-modelling of the cross-frequency covariance matrix.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Correlation matrix $r_{ij}$ of the total data vector. The first 4 × 4 and the last 5 × 5 blocks correspond to the cross-correlations with the coadded CIB map of the DES and KiDS tomographic bins respectively. The off-diagonal elements between DES and KiDS $C_{\ell s}$ are zero, since both surveys do not overlap on the sky.</figcaption></figure>

## 4.2 Constraints on star formation models

Combining the nine cross-power spectra with the coadded CIB map, we are able to place constraints on the four free parameters of the halo-based SFR model of Section 2.1.3, {$m_0, m_z, \eta_0, \eta_z$}. The resulting 1 and 2-σ contours are shown in red in the left panel of Fig. 7. The data are able to constrain both the current value of the peak mass, $m_0$, as well as its evolution with redshift $m_z$. In turn, the two parameters describing the redshift evolution of the peak efficiency $\eta_*$, {$\eta_0, \eta_z$}, are significantly degenerate and, whereas we are able to constrain $\eta_0, \eta_z$ is completely unconstrained. In particular the data are most sensitive to a linear combination of both parameters that is approximately equal to the value of $\eta_*(z)$ at the peak of the lensing kernel $z = 0.4$ of the galaxy redshift bins used in the analysis (see Fig. 2). We will label this parameter $\eta_{0.4}$

$$ \eta_{0.4} = \eta_0 + \eta_z \frac{0.4}{1 + 0.4}. \qquad (36) $$

The constraints on these three parameters, {$m_0, m_z, \eta_{0.4}$}, from the shear-CIB power spectra, are listed in the first row of Table 2. The right panel in Fig. 7 shows, in red, the 1 and 2σ bounds on the evolution of the two SFR functions, $M_1$ (top) and $\eta_*$ (bottom), derived from our constraints on these three parameters. The figure also shows the best-fit model of M18 as a black solid line. Overall, our constraints on $M_1$ are in good agreement with M18, and recover the predicted $\eta_*$ within ~2σ. The data, however, seem to prefer a consistently higher peak star formation efficiency, particularly at low redshifts.

It is interesting to explore to what extent our constraints are in agreement with other measurements of the star formation history from CIB data. In particular, recently J22 made model-independent measurements of the ‘bias-weighted SFR density’ $\langle b\rho_{SFR} \rangle$ (see equation 25) through cross-correlations of the CIB with galaxy samples in the range $0 \lesssim z \lesssim 2$. Fig. 8 shows the measurements of $\langle b\rho_{SFR} \rangle$ of J22 (black points) together with the 1σ confidence interval on this quantity derived from our constraints on the SFR model parameters (red band). Our constraints are in reasonable agreement

with the $\langle b\rho_{SFR} \rangle$ measurements, with a marginal preference for lower values at high redshifts.

Given the compatibility of our measurements of $C_{\ell}^{\gamma, \text{CIB}}$ with the $\langle b\rho_{SFR} \rangle$ data, we can improve our constraints on the SFR model parameters by combining both data sets. This can be done easily by importance-sampling the MCMC chains obtained from the shear power spectrum measurements. Since, as shown in J22, the uncertainties on the measured $\langle b\rho_{SFR} \rangle$ are Gaussianly distributed, we assign a multiplicative weight to each sample i given by

$$ w_i \propto \exp \left[ -\frac{1}{2} (\hat{\mathbf{b}} - \mathbf{b}_i)^T C_{\mathbf{b}}^{-1} (\hat{\mathbf{b}} - \mathbf{b}_i) \right]. \qquad (37) $$

Here $\hat{\mathbf{b}}$ is a vector with the six measurements of $\langle b\rho_{SFR} \rangle$, $C_{\mathbf{b}}^{-1}$ is the covariance matrix of those measurements (also provided in J22), and $\mathbf{b}_i$ is the prediction for $\langle b\rho_{SFR} \rangle$ using the model parameters of the i-th sample. It is worth noting that this procedure assumes that the $\langle b\rho_{SFR} \rangle$ measurements are uncorrelated with the angular power spectra used in this analysis. This is not entirely true. At $z \lesssim 0.8$, the $\langle b\rho_{SFR} \rangle$ data were measured using data from the DESI Legacy Survey (DELS Dey et al. 2019), which partially overlaps with the KiDS sample used here. However, the $\langle b\rho_{SFR} \rangle$ measurements were carried out using the full area overlap between DELS and the CIB maps used here, of which only a small fraction coincides with the KiDS footprint. Furthermore, the uncertainties in the power spectra measured here receive a significant contribution from shape noise in the KiDS sample, which is uncorrelated with the DELS data. Thus, the actual correlation between the $\langle b\rho_{SFR} \rangle$ measurements and the cross-correlation data used here is minimal. We estimate that the correlation between the CIB-shear power spectra used here, and the DELS × CIB cross-correlations used in Jego et al. (2022) is below 10 per cent for all scales used here. We therefore neglect this correlation in our analysis.

The grey contours in the left panel of Fig. 7 show the constraints on the SFR parameters obtained solely from the $\langle b\rho_{SFR} \rangle$ measurements (which coincide with the results published by J22). The blue contours, in turn, show the constraints obtained from the combination of the shear-CIB $C_{\ell s}$ and the $\langle b\rho_{SFR} \rangle$ measurements. These are reported in the second row of Table 2. The constraints on $m_0, m_z$, and $\eta_{0.4}$ improve by ~22 per cent, 52 per cent, and 43 per cent respectively after including the $\langle b\rho_{SFR} \rangle$ data. The corresponding constraints on the evolution of $M_1$ and $\eta_*$ from this combination are shown in blue in the right panel of Fig. 7. Although $M_1$ is still in good agreement with the best-fit model of M18 (black solid line), adding the $\langle b\rho_{SFR} \rangle$ data accentuates the preference for a higher $\eta_*$, particularly at high redshifts. This, however, is mostly driven by the $\langle b\rho_{SFR} \rangle$ data, since the cosmic shear data used here is not very sensitive to redshifts $z \gtrsim 1.0$. The resulting constraints on $\langle b\rho_{SFR} \rangle$ are shown in blue in Fig. 8. Unsurprisingly, the results are in good agreement with the direct measurements of J22.

Finally, we compare the evolution of the star formation rate density predicted from our model constraints with independent direct measurements of $\rho_{SFR}$ based on measurements of the infrared luminosity function. The results are shown in Fig. 9, where the different points with error bars show the measurements of Sanders et al. (2003), Takeuchi, Yoshikawa & Ishii (2003), Magnelli et al. (2011), Magnelli et al. (2013), Gruppioni et al. (2013), Marchetti et al. (2016), Davies et al. (2016), some collected by Madau & Dickinson (2014). The red band shows the 68 per cent constraints derived from our measurements of the CIB-weak lensing cross-correlation, while the blue band shows the combination with the $\langle b\rho_{SFR} \rangle$ measurements. Our constraints are generally in good agreement with all independent measurements of the SFR density covering the last ~10 billion years
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cross-power spectra between the four DES tomographic bins (from top to bottom) and the three CIB frequency maps (353, 545, and 857 GHz from left to right). The points with error bars show our measurements. As in Fig. 4, the orange solid line shows the best-fit model prediction obtained by analysing all the cross-spectra with the coadded CIB map, while the dotted blue line shows the best fit obtained using only the cross-correlations with the DES data. The red dashed lines show the best-fit model obtained from the cross-correlations of all tomographic bins with each frequency map separately. The best-fit model obtained from the coadded map marginally under-estimates the correlation although, as described in Section 4.3, the resulting model constraints are compatible within the statistical uncertainties. Similar results are found for the KiDS data, which we do not show here for brevity.</figcaption></figure>

of cosmic history. This reinforces the interpretation of the CIB as being sourced primarily by the re-emission of UV light absorbed by dust in star-forming galaxies. The solid black line in the same figure shows the best-fit star formation model found by Maniyar et al. (2021) from the analysis of the CIB auto-correlation. Although our measurements are broadly compatible with this preferred model, the combination of $C_{\ell}^{\gamma, \text{CIB}}$ and $\langle b\rho_{\text{SFR}} \rangle$ data predicts a higher SFRD at redshift ~2. As pointed out by J22, this is mostly driven by the high-redshift $\langle b\rho_{\text{SFR}} \rangle$ data, based on the large-scale cross-correlation between the CIB and quasars. An alternative tomographic cross-correlation with the CIB at high redshifts would thus be useful to confirm or rule out this trend. This could be achieved through the cross-correlation with maps of the CMB lensing convergence, assuming that the systematics arising from cross-contamination with the CIB can be kept under control.

## 4.3 Robustness tests

We have quantified the robustness of the results presented in the previous sections to various choices made in our analysis, as well as various potential sources of systematic uncertainty. First, in order

to test for the potential impact of systematics in the cosmic shear measurements, which would affect the DES and KiDS data sets differently, we have re-run MCMC chains using only one of the two data sets. The resulting constraints are shown in the left panel of Fig. 10 in red and blue for DES and KiDS, respectively, together with our fiducial combined constraints as black solid lines. We show results for the three parameters that our data are able to constrain: $\{m_0, m_z, \eta_{0.4}\}$, and the corresponding $1\sigma$ bounds are listed in Table 2. The constraints found with either data set are compatible within ~1$\sigma$. However, we find that the KiDS data set has a preference for lower peak efficiency mass values. This is in line with the results from the cosmological analysis of the KiDS collaboration (Asgari et al. 2021; Heymans et al. 2021), where the weak lensing amplitude, parametrized by $S_8 = \sigma_8\sqrt{\Omega_m/0.3}$, was found to be lower than the preferred Planck value at the ~3$\sigma$ level. Since we use the best-fit Planck cosmology as fiducial in our analysis, a lower halo mass value (and hence a lower halo bias) is required in order to match the KiDS data. Cosmological constraints from weak lensing data are sensitive to systematics in the calibrated redshift distributions used to interpret them. Although we have made use of the official redshift distributions, and marginalized over systematic uncertainty
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left: constraints on the parameters of the M18 halo-based SFR model, described in Section 2.1.3. The parameter $\eta_0$ describes the maximum efficiency at which gas is transformed into stars at low redshifts, while $\eta_z$ parametrizes its redshift evolution. In turn, $m_0$ and $m_z$ parametrize the low-redshift value and evolution of the logarithmic halo mass at which this peak efficiency is attained. Results are shown for our fiducial measurements of the cross-correlation between CIB and cosmic shear data (red), and for its combination with the measurements of $\langle b\rho_{SFR} \rangle$ made by J22 (blue). The constraints from the $\langle b\rho_{SFR} \rangle$ data alone are shown in grey. We find that, while we are able to constrain $m_0$ and $m_z$, we can only accurately constrain a linear combination of the two efficiency parameters, corresponding to the value of the peak efficiency at $z \simeq 0.4$. Right: corresponding constraints on the evolution of the peak efficiency and corresponding peak mass (lower and upper panels) from the shear cross-correlation data (red), and from their combination with the $\langle b\rho_{SFR} \rangle$ measurements. The solid black line shows the best-fit model of M18.</figcaption></figure>

<table><caption>**Table 2.** Constraints on the SFR halo model parameters, characterizing the mass and redshift dependence of the efficiency with which gas is transformed into star (see Section 2.1.3), for various analysis choices.</caption><thead><tr><th>Data</th><th>m<sub>0</sub></th><th>m<sub>z</sub></th><th>η<sub>0.4</sub></th></tr></thead><tbody><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub></td><td>11.90 ± 0.32</td><td>0.0 ± 1.6</td><td>0.519<sup>+0.087</sup><sub>-0.18</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub> + &lt;bρ<sub>SFR</sub>&gt;</td><td>12.17 ± 0.25</td><td>-0.55<sup>+0.66</sup><sub>-0.87</sub></td><td>0.445<sup>+0.055</sup><sub>-0.11</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, DES</td><td>12.19 ± 0.39</td><td>-1.5<sup>+1.9</sup><sub>-2.3</sub></td><td>0.53<sup>+0.12</sup><sub>-0.26</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, KiDS</td><td>11.25<sup>+0.59</sup><sub>-0.44</sub></td><td>1.9 ± 2.1</td><td>0.65<sup>+0.16</sup><sub>-0.25</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, KiDS (no bin 2)</td><td>11.31<sup>+0.66</sup><sub>-0.50</sub></td><td>1.5<sup>+2.7</sup><sub>-2.3</sub></td><td>0.66<sup>+0.17</sup><sub>-0.25</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, 353 GHz</td><td>12.27<sup>+0.37</sup><sub>-0.31</sub></td><td>-1.2<sup>+1.6</sup><sub>-2.0</sub></td><td>0.80 ± 0.20</td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, 545 GHz</td><td>11.90<sup>+0.42</sup><sub>-0.33</sub></td><td>-0.3<sup>+1.8</sup><sub>-2.3</sub></td><td>0.64<sup>+0.15</sup><sub>-0.27</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, 857 GHz</td><td>11.88 ± 0.31</td><td>0.1 ± 1.6</td><td>0.515<sup>+0.088</sup><sub>-0.17</sub></td></tr><tr><td>C<sup>γ,CIB</sup><sub>ℓ</sub>, no IA</td><td>11.62 ± 0.42</td><td>1.7 ± 1.7</td><td>0.398<sup>+0.057</sup><sub>-0.10</sub></td></tr></tbody></table>

parameters as done by both collaborations, there is some evidence of potential mis-calibration of the DES $p(z)$s that could lower $S_8$ to values that are more in line with the KiDS data (Joudaki et al. 2020). Regardless of this, the differences between both data sets observed here are small in relation with the statistical uncertainties and are compatible with a statistical fluctuation. The main conclusion, nevertheless, is that the constraints on SFR model parameters depend on the background cosmological model assumed. A joint cosmology-SFR analysis making use of all auto- and cross-correlations of CIB and cosmic shear data would therefore be able to robustly account for cosmological uncertainties.

As an additional test for systematics in the data, we have quantified the significance of the correlations of the CIB intensity with the

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Bias-weighted SFR density as a function of redshift. The measurements of J22 are shown as black dots with error bars. The $1\sigma$ constraints derived from our measurements of the shear-CIB cross-correlation are shown as a shaded red band, while the constraints combining both data sets are shown in blue.</figcaption></figure>

weak lensing B-modes. The details are described in Appendix A. Although we find that the B-mode correlations are compatible with zero, the second KiDS redshift bin yields consistently low $p$-values in all its correlations with the different CIB maps. Although this might be a statistical fluke, the KiDS collaboration also pointed out a similar discrepancy for this bin in their cosmic shear analysis (Asgari
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Star formation rate density as a function of redshift. The points with error bars correspond to the direct measurements from the infrared luminosity function by various authors (see legend and main text). The solid black line shows the best-fit prediction from the CIB analysis of Maniyar et al. (2021). Our constraints from the cross-correlation between the CIB and cosmic shear data is shown as a red shaded band. Adding these to the measurements of the bias-weighted SFR density $\langle b\rho_{SFR} \rangle$ of J22 yields the blue shaded band (both bands show the 68 per cent confidence intervals).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Robustness tests of our main results. The left panel shows the constraints on the SFR halo model parameters found using our full data vector of cross-correlations (solid black lines), while the red and blue contours show the results of considering only the cross correlations with DES or KiDS respectively. The right panel shows the dependence of the final results on the choice of CIB map. Our fiducial results from the coadded map are shown as solid black lines, while the results using the 353, 545, and 857 GHz maps are shown in grey, red, and blue, respectively.</figcaption></figure>

et al. 2021). Since the E-mode cross-correlation with this bin is not detected with any significance, and it contributes negligibly to the final constraints, we have included these correlations in our final analysis. However, as an a posteriori test, we have re-derived our

constraints after removing these from the data vector. As shown in the fifth row of Table 2, removing these data leads to shift in the final parameter values of ~0.5σ or less compared to the constraints used with the full KiDS complement (but without DES). Since our
constraints are in fact driven by the DES data, given its larger overlap area, we conclude that including this bin in the final analysis does not affect our final results significantly.

Our parametrization of the CIB intensity in terms of $\rho_{SFR}$ requires a good description of the mean infrared spectrum at the relevant frequencies. For this we used the measurements of Béthermin et al. (2013), Béthermin et al. (2015), and Béthermin et al. (2017). In order test the sensitivity of our results to potential inaccuracies in this model, as well as the impact of contamination in the CIB maps from other astrophysical sources (e.g. Sunyaev-Zel'dovich, radio point sources), we have repeated our analysis replacing the coadded CIB map with one of the three single-frequency maps released by Lenz et al. (2019). The results are shown in the right panel of Fig. 10 for the 353, 545, and 857 GHz channels in grey, red, and blue, respectively, together with our fiducial constraints, shown as black lines. Since, as described in Section 3.2, our coadded map receives its main contribution from the 857 GHz map, it is not surprising that the constraints obtained with either of these two maps are very similar. It is reassuring, however, to find that these constraints are also highly compatible with those derived from the 545 or 353 GHz maps alone, albeit with larger error bars. In more detail, the data from the 353 GHz map seem to prefer a model with a larger peak efficiency mass: $m_0 = 12.27_{-0.31}^{+0.37}$ as opposed to $11.90 \pm 0.32$ in our fiducial constraints, corresponding to a $\sim 1\sigma$ shift. Although this is a small shift, in principle, the strong correlation between the three frequency maps makes the difference more significant. This may be a signature of a miscalibration of the effective CIB spectra used in our model, contamination from other astrophysical foregrounds in the 353 map, or that a more sophisticated model is needed to link CIB intensity and SFR density (e.g. a mass dependence in the shape of the mean infrared spectrum). Nevertheless, these differences are small compared to the final statistical uncertainties for current data.

Finally, to quantify the impact of intrinsic alignments, we have rederived constraints on the SFR parameters setting the contribution from IAs to zero. The results are shown in the last row of Table 2. Ignoring intrinsic alignments leads to a $\sim 1\sigma$ downwards shift in both the low-redshift peak efficiency $\eta_{0.4}$, and the corresponding halo mass $m_0$. This is in line with what one would expect: assuming a positive IA amplitude, corresponding to the case where galaxy shapes are aligned with the principal eigenvectors of the tidal field, intrinsic alignments have a negative contribution to the cosmic shear signal. Without accounting for this contribution, the only way to fit the amplitude of the measured correlation with the CIB is to lower the peak halo mass (associated with a lower halo bias), or the peak efficiency itself. Given the sensitivity of current data, the effect is relatively small, but it will become more relevant as better weak lensing data become available. The simple LNLA model used here to describe IAs is therefore good enough for our purposes, but more sophisticated frameworks may be required in the future.

## 5 CONCLUSIONS

In this paper we have studied the cross-correlation between maps of the CIB and tomographic measurements of the matter inhomogeneities from cosmic shear data. This data combination is complementary to the study of the CIB auto-correlation, as well as cross-correlations with galaxies, or with the gravitational lensing of the CMB. Since weak lensing probes directly the matter overdensities, this cross-correlation is able to directly test the connection between matter and star formation rate density, in a way that is independent of the complicated relationship between galaxies and matter that plague the cross-correlation with galaxies on small scales. Unlike in the case

of CMB lensing, the availability of galaxy redshifts allows cosmic shear cross-correlations to probe the redshift dependence of the CIB, although the cumulative nature of gravitational lensing allows for less flexibility than galaxy cross-correlations. Finally, since cosmic shear catalogs are independent of the CMB data that often forms the basis of CMB lensing and CIB maps, this cross-correlation is less sensitive to instrumental systematics in the CMB data, or to cross-contamination from other astrophysical foregrounds. This cross-correlation is therefore a useful and versatile tool to validate physical models of the CIB, its connection with star formation, and its relation to the underlying matter density fluctuations.

We have presented measurements of the shear-CIB cross-correlation using the CIB maps of Lenz et al. (2019), and cosmic shear catalogs from the first data release of the Dark Energy Survey (Zuntz et al. 2018), and from the fourth data release of the Kilometer Survey (Giblin et al. 2021). The cross-correlation is detected with a significance of $\sim 20\sigma$, which constitutes, to our knowledge, the first high-sensitivity measurement of this signal. Furthermore, we have shown that the measurements can be accurately described by a halo-based model that links the CIB intensity and the star formation rate density within the standard $\Lambda$CDM cosmological model, and which is compatible with other studies of the CIB anisotropies.

Our measurements have allowed us to place constraints on free parameters of the model that describe the efficiency with which gas is converted into stars as a function of halo mass. In particular, we are able to constrain the halo mass at which this efficiency peaks today, as well as its time evolution, and the value of the peak efficiency around $z \sim 0.4$, $\eta_{0.4}$. By combining our cross-correlation data with the measurements of the bias-weighted SFR density of J22, we find that the peak efficiency is $\eta_{0.4} = 0.445_{-0.11}^{+0.055}$ and that the corresponding peak mass is

$$ \log_{10}(M_1/M_{\odot}) = 12.17 \pm 0.25 + (-0.55_{-0.87}^{+0.66}) \frac{z}{1+z}. \quad (38) $$

These results, shown in Fig. 7, are in reasonable agreement with the best-fit model of M18, although with a marginal preference for a higher peak efficiency. Qualitatively, the results agree also with the peak efficiency and corresponding mass found by Maniyar et al. (2021), Yan et al. (2022), and confirms the overall picture obtained by previous studies, by which about 40 per cent of all infalling gas is transformed into stars in halos with mass $M \sim 10^{12.5} M_{\odot}$. We have shown that these results are robust with respect to the choice of CIB map and cosmic shear catalog, and are therefore unlikely to be significantly affected by systematics in either data set. Most importantly, as shown in Fig. 9, using our parameter constraints to predict the star formation rate density as a function of time from $z \sim 2$, we find excellent agreement with independent measurements of $\rho_{SFR}$ from the luminosity function of infrared galaxies, further confirming the physical origin of the CIB.

The analysis presented here has made use of a number of assumptions, that should be more thoroughly tested in the future. First, we have ignored the effect of baryonic feedback on the matter power spectrum when predicting the cosmic shear signal. Although this is likely a subdominant effect on the scales $\ell < 1500$ used here, and the theoretical uncertainties in the simple SFR model we use are also probably larger, a more thorough characterization of the impact of baryons would be useful, particularly in the presence of more constraining cosmic shear data. It would also be interesting to explore whether CIB-derived constraints on star formation could indirectly help constrain feedback effects (e.g. by quantifying the diffuse baryon fraction not accreted into stars, or the rate at which massive, short-lived stars are formed). Secondly,
our model assumes a very simple scale dependence of the density of satellite galaxies, following the dark matter distribution through the same NFW profile. This is likely not accurate in detail, and may affect future, more sensitive observations using smaller-scale data, deeper in the 1-halo regime. We leave such a study for future work. Thirdly, we have shown that intrinsic alignments can affect our constraints at the $1\sigma$ level. Here we have used a rather simple model to describe IA's, assuming a direct proportionality with the local tidal field. Although this should suffice for the data studied here, future, more sensitive data sets will likely require the use of more sophisticated models. Furthermore, on the CIB side, we have assumed a simple linear relation to link SFR and infrared luminosity (Chabrier 2003). The value of the proportionality constant, and thus the link between the amplitude of CIB anisotropies and SFR density, crucially depends on the form of the initial mass function at large masses. Although this systematic is common to any constraints on the SFR history, and not just those based on the CIB, the accuracy of this assumption should be more thoroughly tested if a precise understanding of the physics of star formation, and its repercussions on other astrophysical processes, is pursued from these data. Finally, our analysis has assumed perfect knowledge of the background cosmological parameters. This is justified, since these parameters are known to better precision than the SFR model parameters constrained here, thanks to observations of the CMB, large-scale structure, and supernovae. However, given the existing tension between CMB and low-redshift in the values of some of these parameters, particularly in the case of the weak lensing amplitude, a joint analysis of the cosmic shear and CIB data, including all auto- and cross-correlations, and targeting both cosmological and SFR parameters, would allow us to account for cosmological model uncertainties in a self-consistent manner. Furthermore, since photometric redshift systematics may play a significant role in the context of these tensions, future studies of the CIB-shear correlation will also require thorough characterization and marginalization of redshift distribution uncertainties.

In spite of these caveats, our work has shown that important insight on the star formation history, and its dependence on halo mass, can be gained by combining CIB and cosmic shear data. This type of analysis has a promising future, with the advent of new ground-based, high-resolution facilities targeting the far infrared, such as CCAT-prime (Stacey et al. 2018) or the Simons Observatory (Ade et al. 2019), and wide and deep optical weak lensing surveys such as the Rubin Observatory Legacy Survey of Space and Time (LSST, LSST Science Collaboration et al. 2009). The growth in sensitivity and area enabled by these experiments, will vastly increase the significance of the CIB-shear signal, and the range of scales over which it can be effectively employed. In combination with other cross-correlations, including galaxy, CMB lensing, and Sunyaev-Zel'dovich data, this will allow us to improve our understanding of the physical relation between star formation and the matter distribution while simultaneously throwing light onto some of the astrophysical modelling uncertainties described above.

## ACKNOWLEDGEMENTS

The authors would like to thank Martin Rey, Aprajita Verma, and Ziang Yan for useful discussions. BJ is supported by the ENS Paris-Saclay. DA is supported by the Science and Technology Facilities Council through an Ernest Rutherford Fellowship, grant reference ST/P004474. CGG is supported by European Research Council Grant No: 693024 and the Beecroft Trust. JRZ is supported by an STFC doctoral studentship. The authors made extensive use of computational resources at the University of Oxford Department of

Physics, funded by the John Fell Oxford University Press Research Fund.

The authors made extensive use of the numpy (Oliphant 2006; Van Der Walt, Colbert & Varoquaux 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.

## 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 of the relevant papers, as described in the text.

## REFERENCES

* Abazajian K. N. et al., 2016, preprint (arXiv:1610.02743)

* Ade P. et al., 2019, J. Cosmology Astropart. Phys., 2019, 056

* Aiola S. et al., 2020, J. Cosmology Astropart. Phys., 2020, 047

* Alonso D., Sanchez J., Slosar A., LSST Dark Energy Science Collaboration, 2019, MNRAS, 484, 4127

* Asgari M. et al., 2021, A&A, 645, A104

* Astropy Collaboration et al., 2013, A&A, 558, A33

* Astropy Collaboration et al., 2018, AJ, 156, 123

* Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

* Béthermin M., Wang L., Doré O., Lagache G., Sargent M., Daddi E., Cousin M., Aussel H., 2013, A&A, 557, A66

* Béthermin M. et al., 2015, A&A, 573, A113

* Béthermin M. et al., 2017, A&A, 607, A89

* Bianchini F. et al., 2020, ApJ, 888, 119

* Brown M. L., Taylor A. N., Hambly N. C., Dye S., 2002, MNRAS, 333, 501

* Calabrese E. et al., 2017, Phys. Rev. D, 95, 063525

* Cao Y., Gong Y., Feng C., Cooray A., Cheng G., Chen X., 2020, ApJ, 901, 34

* Chabrier G., 2003, PASP, 115, 763

* Chen C.-C. et al., 2016, ApJ, 831, 91

* Chisari N. E. et al., 2019, ApJS, 242, 2

* Cooray A., Sheth R., 2002, Phys. Rep., 372, 1

* Cordero J. P. et al., 2022, MNRAS, 511, 2170

* Darwish O., Sherwin B. D., Sailer N., Schaan E., Ferraro S., 2021a, preprint (arXiv:2111.00462)

* Darwish O. et al., 2021b, MNRAS, 500, 2250

* Davies L. J. M. et al., 2016, MNRAS, 461, 458

* Dey A. et al., 2019, AJ, 157, 168

* Dole H. et al., 2006, A&A, 451, 417

* Dunkley J. et al., 2011, ApJ, 739, 52

* Dunkley J. et al., 2013, J. Cosmology Astropart. Phys., 2013, 025

* Dwek E. et al., 1998, ApJ, 508, 106

* Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

* García-García C., Alonso D., Bellini E., 2019, J. Cosmology Astropart. Phys., 2019, 043

* 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

* Giblin B. et al., 2021, A&A, 645, A105

* Gruppioni C. et al., 2013, MNRAS, 432, 23

* Hazumi M. et al., 2020, Proc. SPIE, 11443, 114432F

* Heymans C. et al., 2021, A&A, 646, A140

* H14PI Collaboration et al., 2016, A&A, 594, A116

* Hinshaw G. et al., 2013, ApJS, 208, 19

* Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526

* Hivon E., Górski K. M., Netterfield C. B., Crill B. P., Prunet S., Hansen F., 2002, ApJ, 567, 2

* Hunter J. D., 2007, Comput.Sci. Eng., 9, 90
Jego B., Ruiz-Zapatero J., García-García C., Koukoufilippas N., Alonso D., 2022, preprint (arXiv:2206.15394)

Joachimi B. et al., 2021, A&A, 646, A129

Joudaki S. et al., 2020, A&A, 638, L1

Kennicutt Robert C. J., 1998, ARA&A, 36, 189

Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531

Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2

Knox L., Cooray A., Eisenstein D., Haiman Z., 2001, ApJ, 550, 7

Komatsu E. et al., 2009, ApJS, 180, 330

Komatsu E. et al., 2011, ApJS, 192, 18

Koukoufilippas N., Alonso D., Bilicki M., Peacock J. A., 2020, MNRAS, 491, 5464

Lenz D., Doré O., Lagache G., 2019, ApJ, 883, 75

Lewis A., 2019, preprint (arXiv:1910.13970)

Limber D. N., 1953, ApJ, 117, 134

LSST Science Collaboration et al., 2009, preprint (arXiv:0912.0201)

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D., Frayer D. T., Willmer C. N. A., 2011, A&A, 528, A35

Magnelli B. et al., 2013, A&A, 553, A132

Maniyar A. S., Béthermin M., Lagache G., 2018, A&A, 614, A39

Maniyar A., Béthermin M., Lagache G., 2021, A&A, 645, A40

Marchetti L. et al., 2016, MNRAS, 456, 1999

McCarthy F., Madhavacheril M. S., 2021, Phys. Rev. D, 103, 103515

Mead A. J., 2017, MNRAS, 464, 1282

Mead A. J., Verde L., 2021, MNRAS, 503, 3095

Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401

Miller L. et al., 2013, MNRAS, 429, 2858

Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822

Münchmeyer M., Madhavacheril M. S., Ferraro S., Johnson M. C., Smith K. M., 2019, Phys. Rev. D, 100, 083508

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

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 T. E., 2006, A guide to NumPy. Vol. 1. Trelgol Publishing, USA

Osborne S. J., Hanson D., Doré O., 2014, J. Cosmology Astropart. Phys., 2014, 024

Partridge R. B., Peebles P. J. E., 1967, ApJ, 147, 868

Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144

Planck Collaboration XVIII, 2011, A&A, 536, A18

Planck Collaboration XVI, 2014a, A&A, 571, A16

Planck Collaboration XVIII, 2014b, A&A, 571, A18

Planck Collaboration XXX, 2014c, A&A, 571, A30

Planck Collaboration XIII, 2016, A&A, 594, A13

Planck Collaboration VI, 2020, A&A, 641, A6

Puget J. L., Abergel A., Bernard J. P., Boulanger F., Burton W. B., Desert F. X., Hartmann D., 1996, A&A, 308, L5

Sailer N., Schaan E., Ferraro S., Darwish O., Sherwin B., 2021, Phys. Rev. D, 104, 123514

Sanders D. B., Mazzarella J. M., Kim D. C., Surace J. A., Soifer B. T., 2003, AJ, 126, 1607

Seljak U., 2000, MNRAS, 318, 203

Serra P., Lagache G., Doré O., Pullen A., White M., 2014, A&A, 570, A98

Shang C., Haiman Z., Knox L., Oh S. P., 2012, MNRAS, 421, 2832

Sheldon E. S., Huff E. M., 2017, ApJ, 841, 24

Sherwin B. D., Schmittfull M., 2015, Phys. Rev. D, 92, 043005

Sievers J. L. et al., 2013, J. Cosmology Astropart. Phys., 2013, 060

Silk J., 2003, MNRAS, 343, 249

Smith K. M., Ferraro S., 2017, Phys. Rev. Lett., 119, 021301

Smith R. E. et al., 2003, MNRAS, 341, 1311

Smith K. M., Hanson D., LoVerde M., Hirata C. M., Zahn O., 2012, J. Cosmology Astropart. Phys., 2012, 014

Smith K. M., Madhavacheril M. S., Münchmeyer M., Ferraro S., Giri U., Johnson M. C., 2018, preprint (arXiv:1810.13423)

Spergel D. N. et al., 2003, ApJS, 148, 175

Spergel D. N. et al., 2007, ApJS, 170, 377

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. SPIE, Bellingham, Washington USA, p. 107001M

Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152

Takeuchi T. T., Yoshikawa K., Ishii T. T., 2003, ApJ, 587, L89

Tessore N., Harrison I., 2020, Open J. Astrophys., 3, 6

Tinker J. L., Wetzel A. R., 2010, ApJ, 719, 88

Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878

Tinsley B. M., 1980, Fund. Cosmic Phys., 5, 287

Torrado J., Lewis A., 2019, Astrophysics Source Code Library, record ascl:1910.019

Torrado J., Lewis A., 2021, J. Cosmology Astropart. Phys., 2021, 057

Tröster T. et al., 2022, A&A, 660, A27

Troxel M. A. et al., 2018, Phys. Rev. D, 98, 043528

Tucci M., Desjacques V., Kunz M., 2016, MNRAS, 463, 2046

Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng., 13, 22

van Engelen A., Bhattacharya S., Sehgal N., Holder G. P., Zahn O., Nagai D., 2014, ApJ, 786, 13

Viero M. P. et al., 2013, ApJ, 772, 77

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wang L. et al., 2015, MNRAS, 449, 4476

Yan Z., van Waerbeke L., Wright A. H., Bilicki M., Gu S., Hildebrandt H., Maniyar A. S., Tröster T., 2022, A&A, 665, A52

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Software, 4, 1298

Zuntz J. et al., 2018, MNRAS, 481, 1149

## APPENDIX: B-MODE POWER SPECTRA

As an additional test for systematics in the data used here, we have examined the correlations between the shear B-modes in each redshift bin of the DES and KiDS catalogs, and the coadded CIB map. To quantify the evidence for a potential systematic, we calculate the $\chi^2$ value of each power spectrum with respect to a null signal, and its associated p-value. The results are shown in Fig. A1. We find overall acceptable p-values ($p > 2$ per cent), with the exception of the second KiDS bin, for which $p = 0.001$. This is in agreement with the findings of the KiDS team, which identified evidence of systematics in this bin. As we have shown in the main text, the inclusion of this bin in the analysis has a negligible effect in our final results, since the shear signal at these low redshifts is small.
<figure><img src="image_12.png" /><figcaption>**Figure A1.** B-mode power spectra for the DES and KiDS data sets (top and bottom panels). The probability-to-exceed of each power spectrum, corresponding to its $\chi^2$ value with respect to a null signal is included in each panel.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.