SIDES: Simulated Infrared Dusty Extragalactic Sky

• SIDES is a simulation framework that models the dusty extragalactic sky using a joint treatment of galaxy populations, clustering, and observational biases.
• The framework uses subhalo abundance matching and evolving SED libraries to realistically reproduce far-infrared to millimeter number counts and redshift distributions.
• SIDES forward-models observational processes including beam convolution, confusion, and source extraction, making it critical for interpreting Herschel and JWST data.

Searching arXiv for SIDES-related papers to ground the article in current literature. arxiv_search(query="SIDES simulated infrared dusty extragalactic sky Bethermin 2017 (Bethermin et al., 2017)", max_results=5) SIDES denotes a semi-empirical lightcone framework for simulating the dusty extragalactic sky from the mid-infrared to the millimeter. In the literature, the acronym is expanded both as the “Simulated Infrared Dusty Extragalactic Sky” and the “Simulated Infrared Extragalactic Dusty Sky.” Its defining feature is the joint treatment of galaxy populations, large-scale structure, spectral energy distributions, beam convolution, and source extraction, so that intrinsic catalogs and observed maps can be compared on the same footing. The framework was introduced to explain how confusion, blending, clustering, and finite angular resolution reshape far-infrared and submillimeter observables, and it was later extended to sub-mm/mm spectroscopy, polarization forecasts, and JWST/MIRI source-population modeling (Bethermin et al., 2017, Bethermin et al., 2022, Béthermin et al., 2024, Vidal et al., 18 Sep 2025).

1. Historical setting and scientific scope

Before SIDES, a closely related infrared–radio sky simulation had already shown that an invariant mapping between radio luminosity, far-infrared luminosity, and infrared SEDs was insufficient to reproduce the full set of Spitzer and SCUBA number counts. In that earlier work, radio luminosities were converted into far-infrared luminosities through the Yun et al. $q$ relation,

$$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$

with a Gaussian distribution $p=2.34$, $\sigma=0.26$, and infrared templates were assigned from the Rieke et al. library. The model required stronger luminosity evolution for star-forming galaxies than in the original radio simulation, and it also required cooler templates with stronger PAH emission features at increasing redshift at fixed far-infrared luminosity, with a reversal from $z=1$ to 2 to avoid overproducing 850 $\mu$m counts (Wilman et al., 2010).

SIDES was constructed in a different but complementary context. Its immediate purpose was to reproduce the far-infrared to millimeter sky in a way that was consistent with both intrinsic galaxy populations and the observational biases introduced by Herschel-like and ground-based beams. The framework was then broadened in two directions: toward spectral-cube simulations for [CII], CO, and [CI] intensity mapping, and toward shorter-wavelength mid-infrared applications where quiescent galaxies, stellar photospheric emission, composites, and dust-obscured AGN become important (Bethermin et al., 2017, Bethermin et al., 2022, Vidal et al., 18 Sep 2025).

Paper	Scope	Distinctive addition
(Bethermin et al., 2017)	2 deg$^2$ far-IR to submm mock sky	Clustering, map-making, source extraction, resolution bias
(Bethermin et al., 2022)	Sub-mm/mm spectroscopy	CO, [CII], [CI], CONCERTO-like cubes, pySIDES
(Béthermin et al., 2024)	PRIMAger forecasting	Classical confusion in intensity and polarization
(Vidal et al., 18 Sep 2025)	JWST/MIRI regime	Quiescent galaxies, AGN/composites, probabilistic $f_{\rm AGN}$

2. Cosmological backbone and galaxy–halo connection

The 2017 SIDES realization is built on a Bolshoi-Planck lightcone of $1.4^\circ \times 1.4^\circ$, spanning $0$$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$0, with halo mass resolution down to $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$1. The simulation predicts continuum observables from $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$2 to $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$3, with particular emphasis on Herschel PACS/SPIRE wavelengths, $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$4, and 1.1–1.2 mm (Bethermin et al., 2017).

SIDES uses subhalo abundance matching to connect galaxies to the dark-matter lightcone. Rather than matching stellar mass directly to halo mass, it ranks halos and subhalos by the peak circular velocity $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$5, ranks galaxies by stellar mass, assigns them monotonically, and includes an intrinsic scatter of 0.2 dex in stellar mass at fixed $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$6. The observed stellar mass function is deconvolved for this scatter before matching. This procedure is explicitly used because the paper states that the 0.2 dex scatter is needed to reproduce observed galaxy clustering (Bethermin et al., 2017).

The stellar mass function is written as a double Schechter form,

$$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$7

with parameters interpolated smoothly with redshift from observed stellar mass functions. Galaxies are partitioned into passive and star-forming populations through

$$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$8

and star-forming systems are assigned SFRs using the Schreiber et al. main-sequence prescription. In the original SIDES far-infrared/submm implementation, passive galaxies are assumed to have negligible FIR/submm output (Bethermin et al., 2017).

Clustering is therefore not inserted as an imposed angular correlation function. It emerges from the halo distribution in the Bolshoi-Planck simulation after abundance matching. This choice is central to SIDES because the framework is designed not only to match number counts but also to model blending, beam contamination, and map-level observables (Bethermin et al., 2017).

3. Emission modeling and SED libraries

The continuum version of SIDES is based on an updated 2SFM, or two star-formation modes, model. In this framework, star-forming galaxies occupy a main sequence plus a starburst component, while passive galaxies are treated separately. Later descriptions of SIDES specify that the continuum model uses evolving SEDs for main-sequence and starburst galaxies, with scatter in dust temperature and UV radiation field intensity, and that these ingredients were calibrated to observed number counts, redshift distributions, pixel histograms, and CIB power spectra (Bethermin et al., 2017, Bethermin et al., 2022).

The JWST/MIRI extension modified this SED framework substantially. The original dust-only main-sequence and starburst templates were extended to include stellar photospheric emission. The update keeps the Magdis et al. dust templates, varying with radiation field intensity $$q = \log\left[\frac{L(\mathrm{FIR},\,\mathrm{W})}{3.75\times 10^{12}\,\mathrm{Hz}}\right] - \log\left[L_{1.4\mathrm{GHz}}(\mathrm{W\,Hz}^{-1})\right],$$9, and adds stellar light from SWIRE templates: the Sd template for main-sequence galaxies and the M82 template for starbursts. Quiescent galaxies, absent from the original FIR-focused SIDES treatment because their dust emission is weak, are explicitly added through the Ell5 SWIRE template and are normalized to the rest-frame $p=2.34$0-band using a redshift-dependent $p=2.34$1 relation derived from the Santa Cruz semi-analytic model lightcone (Vidal et al., 18 Sep 2025).

AGN and composites are modeled with the Kirkpatrick et al. template library, indexed by the mid-infrared AGN fraction $p=2.34$2 from 0 to 1 in 11 coarse bins. These templates are then modified by stitching in stellar emission at different wavelengths depending on AGN strength, and interpolated linearly to produce 2001 SEDs with $p=2.34$3. The classification used in the paper is star-forming for $p=2.34$4, composite for $p=2.34$5, and AGN for $p=2.34$6 (Vidal et al., 18 Sep 2025).

A further modification is the probabilistic assignment of the infrared AGN fraction. The model defines

$p=2.34$7

ties the AGN infrared luminosity to accretion-disk luminosity through $p=2.34$8 and $p=2.34$9, and then draws $\sigma=0.26$0 probabilistically from PDFs built from the Santa Cruz semi-analytic model in bins of redshift and $\sigma=0.26$1. The empirical cumulative distribution function is

$\sigma=0.26$2

The best overall match to observed AGN fractions is obtained for $\sigma=0.26$3 (Vidal et al., 18 Sep 2025).

4. Forward modeling, map generation, and resolution bias

A distinctive property of SIDES is that it forward-models the observational process rather than stopping at an intrinsic catalog. After stellar masses, SFRs, SEDs, and lensing magnifications are assigned, mock maps are generated with Gaussian beams corresponding to Herschel resolutions: 5.5″ at 70 $\sigma=0.26$4m, 6.5″ at 100 $\sigma=0.26$5m, 11″ at 160 $\sigma=0.26$6m, 18.2″ at 250 $\sigma=0.26$7m, 24.9″ at 350 $\sigma=0.26$8m, and 36.3″ at 500 $\sigma=0.26$9m. For 850 $z=1$0m and 1.1 mm, beam-smoothed maps are also produced at roughly 21″ and 26″. For pure count comparisons, no instrumental noise is included; confusion arises from unresolved faint galaxies. Herschel-like catalogs are extracted with FASTPHOT, typically using 24 $z=1$1m positional priors, whereas $z=1$2 blind extraction is used at 850 $z=1$3m and 1.2 mm because 24 $z=1$4m priors are no longer reliable there (Bethermin et al., 2017).

This end-to-end treatment is required because intrinsic and extracted observables differ strongly at long wavelengths. SIDES reproduces number counts from 70 to 500 $z=1$5m and at 850 $z=1$6m and 1.2 mm, redshift distributions across 100 $z=1$7m to 1.1 mm, number counts split by redshift bin, CIB anisotropy power spectra, and the $z=1$8 pixel histogram. The stacking bias is defined as

$z=1$9

and the paper emphasizes that $\mu$0, stacking, and CIB power spectra are all sensitive to clustering (Bethermin et al., 2017).

The main quantitative result is that observed Herschel number counts between 5 and 50 mJy at 350 and 500 $\mu$1m are biased high by about a factor of 2 because single-dish beams blend multiple galaxies. The same simulation simultaneously matches higher-resolution interferometric counts and Herschel single-dish counts once the observational process is modeled consistently. Multiplicity is also quantified directly: at 250 $\mu$2m the brightest real galaxy contributes about 85% of the measured flux, whereas at 500 $\mu$3m it contributes only about 58–60% on average. Extracted redshift distributions are biased low at 350 and 500 $\mu$4m because the brightest 24 $\mu$5m prior in a 500 $\mu$6m beam is not always the galaxy dominating the 500 $\mu$7m emission. SIDES further suggests that the reported excess of very red Herschel sources can be explained largely by noise, blending and resolution effects, steep flux and color distributions, and some contribution from lensing, rather than requiring a failure of the galaxy-evolution model (Bethermin et al., 2017).

This aspect of SIDES is frequently treated as its central methodological contribution: theory should not be compared directly to single-dish counts without forward-modeling the beam, noise, and extraction. The public 2 deg$\mu$8 products were released for that purpose at the SIDES website (Bethermin et al., 2017).

5. Spectroscopic extension and intensity mapping

SIDES was later updated for sub-mm/mm spectroscopy and [CII] intensity mapping. In that version, the original Bolshoi-Planck lightcone and continuum model are retained, but the framework is extended to include the main lines around 1 mm—CO, [CII], and [CI]—together with a spectral-cube generator. The code was rewritten in Python as pySIDES, and both code and products were released publicly (Bethermin et al., 2022).

CO is modeled with a relation between $\mu$9 and $^2$0,

$^2$1

with an intrinsic scatter of 0.2 dex, plus a $^2$2 dex offset for starbursts. [CII] is assigned through two empirical prescriptions, De Looze et al. 2014 and Lagache et al. 2018, again with 0.2 dex scatter in the SFR–[CII] relation. CI and CI are calibrated empirically from observed line compilations. CONCERTO-like cubes are generated over 125–305 GHz with 1 GHz spectral channels and 5 arcsec pixels, in line-only and full continuum-plus-lines forms (Bethermin et al., 2022).

The updated framework reproduces observed CO luminosity functions from COLDz, ASPECS pilot and large programs, and PdBI HDFN; [CI] luminosity functions within 1$^2$3 of ASPECS constraints; ALPINE-based [CII] luminosity functions at $^2$4 and $^2$5 with model-dependent differences; and the mmIME shot-noise measurement around 100 GHz, for which the predicted total CO+[CI] shot noise averaged over the mmIME frequency range is $^2$6 (Bethermin et al., 2022).

For [CII] intensity mapping, the major result is the scale of the astrophysical uncertainty. Depending on the assumed SFR–[CII] relation and the high-redshift star-formation history, the predicted $^2$7 [CII] power spectrum varies by about two orders of magnitude. The dust continuum is by far the brightest component, by a factor of 3 to 100 depending on frequency, while CO is the main line contaminant and [CI] is smaller but non-negligible. Masking known galaxies from deep surveys reduces CO+[CI] by more than an order of magnitude and can lower the residual foregrounds to about 20% of the [CII] power spectrum up to $^2$8, but the paper states that masking alone is not sufficient at higher redshift (Bethermin et al., 2022).

6. Confusion-limited forecasting in intensity and polarization

In forecasts for PRIMAger, SIDES is used as the astrophysical truth catalog behind simulated maps. The authors use the 2 deg$^2$9 SIDES realization and the integrated map maker make_maps.py to build noiseless intensity and polarization maps containing confusion from unresolved galaxies but no instrumental noise. Gaussian beams, rectangular filters, and beam-sampling pixel scales are adopted; no Galactic cirrus is included; maps are expressed in Jy/beam; and source detection is intentionally basic, using photutils without matched filtering, in order to provide a classical confusion benchmark (Béthermin et al., 2024).

The confusion limit is defined as the faintest flux density at which sources can be reliably extracted in the limit of zero instrumental noise. In intensity, the classical confusion limit rises from 21 $f_{\rm AGN}$0Jy at 25 $f_{\rm AGN}$1m to 46 mJy at 235 $f_{\rm AGN}$2m. The measured flux density is dominated by the brightest galaxy in the beam, but the brightest-source contribution decreases from more than 0.91 at the classical confusion limit in the PHI bands to 0.72 in the longest-wavelength PPI band, implying that at 235 $f_{\rm AGN}$3m other objects contribute about 30% of the measured flux. Physically associated sources contribute $f_{\rm AGN}$4; the additional flux is mostly due to unrelated galaxies at other redshifts (Béthermin et al., 2024).

The same study finds that clustering has a mild effect on confusion in intensity, up to 25% in the summary, but negligible impact on polarization confusion and 50% completeness. In polarization, each galaxy is assigned a polarization fraction $f_{\rm AGN}$5 and angle $f_{\rm AGN}$6, and the Stokes parameters are defined by

$f_{\rm AGN}$7

Polarization confusion limits are more than two orders of magnitude lower than in intensity, reaching 0.03, 0.078, 0.17, and 0.25 mJy in the four PPI bands. The paper forecasts roughly 8,000 detections up to $f_{\rm AGN}$8 for the conservative sensitivity case, and emphasizes a strong synergy in which intensity surveys at short wavelength and polarization surveys at long wavelength tend to reach confusion at similar depth (Béthermin et al., 2024).

These PRIMAger results also clarify a recurrent ambiguity in the SIDES literature. Clustering can matter strongly for extracted counts, $f_{\rm AGN}$9, and beam contamination, yet still have only a mild effect on the formal confusion limit in a particular survey configuration. The distinction is observational rather than conceptual: different summary statistics respond differently to the same clustered sky (Bethermin et al., 2017, Béthermin et al., 2024).

7. JWST/MIRI extension and current interpretation

The MIRI update extends SIDES into the 5–25 $1.4^\circ \times 1.4^\circ$0m regime, where stellar photospheric emission, PAH features, hot dust, quiescent galaxies, composites, and dust-obscured AGN all contribute materially. Earlier SIDES versions did not explicitly model quiescent galaxies or AGN and had SEDs that were too limited for the MIRI regime. The revised framework was therefore built to reproduce observed MIRI source number counts, redshift distributions, and AGN population fractions while preserving the earlier SIDES success in the far-infrared (Vidal et al., 18 Sep 2025).

The resulting population mix differs from some earlier theoretical expectations. AGN dominate at bright flux densities, approximately $1.4^\circ \times 1.4^\circ$1, while main-sequence galaxies dominate at the faint end. Composites are the second most important population and are especially significant at intermediate fluxes and at cosmic noon redshifts. Quiescent galaxies contribute mainly at low redshift and never dominate the MIRI counts. The paper also notes that many MIRI-selected “AGN” are in fact AGN+star-formation composites, because AGN-only counts fall below the observed fraction at intermediate fluxes and agreement is restored when composites are included (Vidal et al., 18 Sep 2025).

The MIRI redshift distributions are often multi-modal because the filters repeatedly sample the $1.4^\circ \times 1.4^\circ$2 stellar bump at shorter wavelengths and PAH features at longer wavelengths. The framework therefore provides not only counts but also diagnostic MIRI color–color diagrams and joint NIRCam/MIRI flux distributions. A practical conclusion is that essentially all MIRI sources should have NIRCam counterparts in overlapping fields, but many MIRI-selected AGN and composites remain difficult to identify using NIRCam alone because host-galaxy stellar light dominates in the optical and near-infrared (Vidal et al., 18 Sep 2025).

The update also quantifies cosmic variance explicitly. Surveys smaller than about 20–25 arcmin$1.4^\circ \times 1.4^\circ$3 suffer strong cosmic variance; the abstract states that surveys with areas below $1.4^\circ \times 1.4^\circ$4 suffer from $1.4^\circ \times 1.4^\circ$5 uncertainty in bright AGN counts, and the detailed Monte Carlo sub-sampling shows that a field like SMILES, at $1.4^\circ \times 1.4^\circ$6, still has about 30% fractional uncertainty at the bright end. This places survey area on the same practical footing as depth for bright AGN statistics (Vidal et al., 18 Sep 2025).

Taken together, the SIDES literature defines a simulation framework whose scientific role is dual. It is a galaxy-population model anchored to large-scale dark-matter structure, and it is an observation model that propagates those galaxies through beam smearing, confusion, extraction, and survey selection. Its recurring conclusion is that apparent tensions between theory and data at dusty wavelengths are often inseparable from how the sky is observed.