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

HAL
open science

# Identification of large equivalent width dusty galaxies at $4 < z < 6$ from sub-millimetre colours

D. Burgarella, P. Theulé, V. Buat, L. Gouiran, L. Turco, M. Boquien, T. J. L.

C. Bakx, A. K. Inoue, Y. Fudamoto, Y. Sugahara, et al.

► To cite this version:

D. Burgarella, P. Theulé, V. Buat, L. Gouiran, L. Turco, et al. Identification of large equivalent width dusty galaxies at $4 < z < 6$ from sub-millimetre colours. Astronomy & Astrophysics, 2023, 671, 10.1051/0004-6361/20224491. insu-04479103

HAL Id: insu-04479103

https://insu.hal.science/insu-04479103

Submitted on 3 Mar 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Identification of large equivalent width dusty galaxies at 4 < z < 6 from sub-millimetre colours

D. Burgarella<fnref n="1" />, P. Theulé<fnref n="1" />, V. Buat<fnref n="1" />, L. Gouiran<fnref n="1" />, L. Turco<fnref n="1" />, M. Boquien<fnref n="2" />, T. J. L. C. Bakx<fnref n="3" />,<fnref n="4" />, A. K. Inoue<fnref n="5" />,<fnref n="6" />, Y. Fudamoto<fnref n="4" />, Y. Sugahara<fnref n="4" />,<fnref n="6" />, and J. Zavala<fnref n="4" />

Received 13 July 2022 / Accepted 10 November 2022

## ABSTRACT

**Context.** Infrared (IR), sub-millimetre (sub-mm), and millimetre (mm) databases contain a huge quantity of high-quality data. However, a large part of these data are photometric, and they are thought not to be useful to derive quantitative information on the nebular emission of galaxies.

**Aims.** The aim of this project is first to identify galaxies at $z \ge 4-6$ and in the epoch of reionisation from their sub-millimetre colours. We also aim to show that the colours can be used to try and derive physical constraints from photometric bands when accounting for the contribution from the IR fine structure lines to these photometric bands.

**Methods.** We modelled the flux of IR fine structure lines with CLOUDY and added them to the dust continuum emission with CIGALE. Including (or not) emission lines in the simulated spectral energy distribution (SED) modifies the broad-band emission and colours.

**Results.** The introduction of the lines allows us to identify strong star forming galaxies at $z \ge 4-6$ from the $[\log_{10} \frac{\text{PSW}_{250\mu\text{m}}}{\text{PMW}_{350\mu\text{m}}}]$ versus $\log_{10} \frac{\text{LABOCA}_{870\mu\text{m}}}{\text{PLW}_{500\mu\text{m}}}]$ colour-colour diagram. By comparing the relevant models to each observed galaxy colour, we are able to roughly estimate the fluxes of the lines and the associated nebular parameters. This method allows us to identify a double sequence in a plot built from the ionisation parameter and the gas metallicity.

**Conclusions.** The HII and photodissociation region fine structure lines are an essential part of the SEDs. It is important to add them when modelling the spectra, especially at $z \ge 4-6$, where their equivalent widths can be large. Conversely, we show that we can extract some information on strong-IR fine structure lines and on the physical parameters related to the nebular emission from IR colour-colour diagrams.

**Key words.** galaxies: formation – galaxies: evolution – galaxies: high-redshift – galaxies: ISM – ISM: abundances – submillimeter: ISM

## 1. Introduction

Several papers have reported excesses of the flux densities of high-redshift galaxies in broad bands. For instance, a boost of the *Spitzer*/IRAC bands at $z \sim 7-8$ is observed when $H\alpha$, and [OIII]500.7 nm fall in the mid-infrared (mid-IR) filters (e.g., de Barros et al. 2013; Roberts-Borsani et al. 2020; Anders 2003). In the sub-millimetre (sub-mm) as well, Seaquist et al. (2004) suggested that about 25% of the 850 µm flux density could be due to the CO(3–2) molecular emission. More relevant to this paper, Smail et al. (2011) estimated that a [CII]158 µm fine structure line with 0.27% of the galaxy’s $L_{\text{dust}}$ would contribute 5–10% to the far-IR broadband flux densities at 850 µm. Because the line contribution scales linearly with $L_{\text{[CII]}}/L_{\text{dust}}$, the [CII]158 µm-to-dust-continuum luminosity ratio and sources with $L_{\text{[CII]}}/L_{\text{dust}} > 1\%$ will contribute more than four times this amount, reaching 20–40% of the 850 µm flux densities for galaxies at $z \sim 4-6$ (also see Seymour et al.

(2012) for the contribution to the Herschel/SPIRE 500 µm band).

On the contrary, $L_{\text{[CII]}}/L_{\text{dust}}$ presents a deficit for galaxies with large IR surface brightnesses or IR luminosities. Luhman et al. (2003) proposed that this deficit could be due to high values of the ionisation parameter<fnref n="1" /> ($\log_{10} U > -2.5$), for which a narrower photodissociation region (PDR) would lead to lower [CII]158 µm fluxes. This explanation is also supported by a number of other analyses (e.g., Abel et al. 2009; Díaz-Santos et al. 2013, 2017; Herrera-Camus et al. 2018).

Recent promising papers concerning the James Webb Space Telescope (JWST, e.g., Schaerer et al. 2022; Trump et al. 2022;

<footnote n="1">¹ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
e-mail: denis.burgarella@lam.fr</footnote>

<footnote n="2">² Centro de Astronomía (CITEVA), Universidad de Antofagasta, Avenida Angamos 601, Antofagasta, Chile</footnote>

<footnote n="3">³ Division of Particle and Astrophysical Science, Graduate School of Science, Nagoya University, Aichi 464-8602, Japan</footnote>

<footnote n="4">⁴ National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="5">⁵ Department of Physics, School of Advanced Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan</footnote>

<footnote n="6">⁶ Waseda Research Institute for Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan</footnote>

<footnote n="1">¹ The ionisation parameter is defined as the dimensionless ratio of the incident ionising photon density to the hydrogen density $-U = Q(H)/(4\pi R^2 n_H c)$— where $Q(H)$ is the number of hydrogen ionising photons per second, $c$ is the speed of light, $n_H$ is the hydrogen density, and $R$ is the distance of the ionising source to the illuminated face.</footnote>
<table><caption>**Table 1.** CIGALE modules and input parameters used for to create the 88 million models.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Fit w/o lines</th><th>Fit w/ lines</th></tr></thead><tbody><tr><td>Target sample</td><td></td><td>SPT DSFGs</td><td>SPT DSFGs</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Delayed SFH and recent burst</td></tr><tr><td>e-folding time scale of the delayed SFH</td><td>τ<sub>main</sub> [Myr]</td><td>500</td><td>500</td></tr><tr><td>Age of the main population</td><td>Age<sub>main</sub> [Myr]</td><td>10</td><td>10</td></tr><tr><td>Burst</td><td>f<sub>burst</sub></td><td>No burst</td><td>No burst</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">SSP</td></tr><tr><td>SSP</td><td></td><td>BC03</td><td>BC03</td></tr><tr><td>Initial mass function</td><td>IMF</td><td>Chabrier</td><td>Chabrier</td></tr><tr><td>Metallicity</td><td>Z</td><td>0.02</td><td>0.02</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Nebular emission</td></tr><tr><td>Ionization parameter</td><td>logU</td><td>—</td><td>31 values in [-4.0, -1.0]<br>with δ log U = 0.1</td></tr><tr><td>Gas metallicity</td><td>zgas</td><td>—</td><td>0.0001, 0.001, 0.0025, 0.005,<br>0.007, 0.008, 0.011, 0.014,<br>0.016, 0.019, 0.022, 0.025,<br>0.03, 0.033, 0.037, 0.041,<br>0.046, 0.051</td></tr><tr><td>Electron density</td><td>ne</td><td>—</td><td>10, 100, 1000</td></tr><tr><td>Line width [km s<sup>-1</sup>]</td><td>—</td><td>—</td><td>200</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Dust attenuation law (dustatt_modified_CF00)</td></tr><tr><td>V-band attenuation in the interstellar medium</td><td>Av_ISM</td><td>10</td><td>10</td></tr><tr><td>Av_ISM / (Av_BC+Av_ISM)</td><td>mu</td><td>0.44</td><td>0.44</td></tr><tr><td>Power law slope of the attenuation in the ISM</td><td>slope_ISM</td><td>-0.7</td><td>-0.7</td></tr><tr><td>Power law slope of the attenuation in the birth clouds</td><td>slope_BC</td><td>-1.3</td><td>-1.3</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Dust emission (casey2012)</td></tr><tr><td>Temperature of the dust in K</td><td>temperature</td><td>31 values in [30., 90.]</td><td>31 values in [30., 90.]</td></tr><tr><td>Emissivity index of the dust</td><td>beta</td><td>21 values in [1.0, 3.5]</td><td>21 values in [1.0, 3.5]</td></tr><tr><td>Mid-infrared powerlaw slope</td><td>alpha</td><td>2.0</td><td>2.0</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Redshifting</td></tr><tr><td>Redshift</td><td>redshift</td><td>81 values in [0.0, 8.0]</td><td>81 values in [0.0, 8.0]</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">No AGN emission</td></tr></tbody></table>

**Notes.** CF00 means Charlot & Fall (2000), BC03 means Bruzual & Charlot (2003), Casey2012 means Casey (2012), and the Chabrier IMF refers to Chabrier (2003).

Taylor et al. (2022) suggest that we are now able to spectroscopically measure the nebular parameters of galaxies in the epoch of reionisation (EoR) at $z \sim 5-8$. However, we still have large uncertainties when estimating the obscured star formation density at high redshift (e.g., Algera et al. 2023). Thus, we still need to identify and measure the luminosities of galaxies in the EoR to better understand how the total star formation rate density (SFRD) evolves in the early Universe to understand the formation and early evolution of galaxies at $z > 10$ (e.g., Finkelstein et al. 2022).

In this paper, we propose an original method to identify $4.5 < z < 6.0$ galaxies via a broad-band excess due to the [CII]158 µm line. This method is applicable to very large numbers of dusty galaxies that can be extracted from the already existing far-IR and sub-millimetre databases (Sect. 3). We present tests on a sample of galaxies with spectroscopic redshifts observed with the South Pole Telescope (SPT, Sect. 2) that seem to confirm the validity of the above method (Sect. 4). This method simultaneously provides a way to constrain the nebular parameters of galaxies at $2 < z < 7$ (Sect. 5). Moreover, the same method could be utilised for different redshift ranges and different photometric bands.

We assumed a Chabrier initial mass function (IMF, Chabrier 2003). We used WMAP7 cosmology (Komatsu et al. 2011). Finally, we assumed a solar metallicity $Z_{\odot} = 0.014$ from Asplund et al. (2009).

## 2. The SPT galaxy sample

Reuter et al. (2020) presented the final spectroscopic redshift analysis of a flux-limited ($S_{870\mu m} > 25$ mJy) sample of galaxies from the 1.4 mm SPT survey. In this 2500-square-degree survey observed at 1.4 mm and 2.0 mm, they identified 81 strongly lensed, dusty star-forming galaxies (DSFGs) at $1.9 \le z \le 6.9$. The spectroscopic observations were conducted with the Atacama Large Millimeter/submillimeter Array (ALMA) across the 3 mm spectral window, targeting carbon monoxide line emission. From them, spectroscopic redshifts have been estimated by combining ALMA data with ancillary data. They are used in the rest of this paper, and we did not estimate them when fitting the observed spectral energy distributions (SED) or when estimating the physical parameters.

With APEX/LABOCA, Strandet et al. (2016) obtained 870 µm flux densities in the period from 2010 September – 2012 November. The Herschel/SPIRE maps at 250 µm, 350 µm, and 500 µm were observed in two observing programmes, in the period from 2012 August – 2013 March. As described in Strandet et al. (2016), the flux densities were extracted by fitting a Gaussian to the source. The peak of the Gaussian is taken as the flux density. The noise was estimated by taking the RMS in the central few arcmins of the map, which was then added in
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Evolution of CIGALE models with nebular parameters. This sample of models created by CIGALE is an extract of the entire and much larger one used for the SED fitting (88 × 10⁶ models). Here, in addition to the nebular parameters given in the legend, we fix the SFH (delayed with τmain = 500 Myrs, and age main = 100 Myrs), A_V (ISM) = 0.3, Tdust = 40 K, the dust emissivity β = 2.0, and the redshift z = 0. However, to improve the visibility, we offset the spectra by δz = 0.1 on the X axis and 1 dex on the Y axis. The bottom spectra, with log₁₀ U = -2 are more representative of the IR spectra emitted by an HII-region dominated galaxy, with a strong [OIII]88.3 µm line (blue-shaded area) detected, for example, in Lyman break galaxies in the early Universe. The top spectra, with log₁₀ U = -4, resemble an IR spectra emitted by a PDR-dominated galaxy, with a strong [CII]157.6 µm line (red-shaded area), detected, for example, in DSFGs in the early Universe. A higher gas metallicity amplifies the strength of the metal lines.</figcaption></figure>

quadrature to the uncertainty due to pixelation. The data them-
selves are from Reuter et al. (2020).

This SPT galaxy sample contains bright IR galaxies that
might include an active galactic nucleus (AGN). So, a prelim-
inary check is needed to see if an AGN can substantially modify
the IR spectral energy distribution of galaxies. However, none
of the attempted observations were able to confirm the pres-
ence of an AGN in these objects, or at least not a strong AGN.
De Breuck et al. (2019) derived a low AGN fraction<fnref n="2" /> (frac<sub>AGN</sub> <
5%) for SPT 0418-47, and they claim that this is consistent with
the AGN fractions estimated in other DSFGs from the SPT sam-
ple. In the galaxy SPT 2132-58, Béthermin et al. (2016) identi-
fied an evolved interstellar medium (ISM) with 0.5 < Z/Z<sub>☉</sub>
1.5, dominated by PDRs. Its CO spectral line energy distri-
bution does not allow us to reach a conclusion on the pres-
ence or absence of an AGN. Gururajan et al. (2022) showed that
the presence of an AGN in SPT 0103-45 is unlikely (see also,
Ma et al. 2016; Spilker et al. 2018). Finally, Apostolovski et al.
(2019) found no signs for a signature pointing to the presence of
an AGN in SPT 0346-52. Thus, for the sake of the present paper,
we made the assumption that there are no major AGNs in the
studied sample.

3. A colour-colour approach for the selection of galaxies at z > 4

We modelled the Herschel and APEX/LABOCA colours
log₁₀ PSW₂₅₀µm / PMW₃₅₀µm and log₁₀ LABOCA₈₇₀µm / PLW₅₀₀µm with CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019). With the input

CIGALE parameters listed in Table 1, CIGALE creates 88 mil-
lion models. This particular set of colours is selected to try and
identify line-boosted galaxies at z ≥ 4–6 as outliers in the above
colour-colour diagram.

A grid of nebular emission lines that includes HII regions and PDRs are pre-computed with CLOUDY (Ferland et al. 2017) by varying the nebular parameters N<sub>e</sub>, Z<sub>gas</sub>, and U. These emission lines are included into CIGALE when building the modelled spectra (Fig. 1). The photo-ionising field shape is generated with the single stellar population (SSP) model library (Bruzual & Charlot 2003) using a constant star formation history (SFH) over 10 Myr and accounting for a range of metallicities, ionisation parameters, and number densities of hydrogen. The radiation field intensity is given by the dimensionless ionisation parameter U ≡ n<sub>γ</sub>/n<sub>H</sub>, where n<sub>γ</sub> is the number density of photons capable of ionising hydrogen, and n<sub>H</sub> the number density of hydrogen, which equals the number density of electrons in a fully ionised medium. The limit of the effective HII region is set by an ionisation fraction ≤10<sup>-3</sup>, while the end of the effective PDR is set by the visual attenuation A<sub>V</sub> ≤ 10 (Röllig et al. 2007), and as such includes part of the molecular region. The size and abundance distribution of grains, typical for the ISM of our galaxy, is used to account for the extinction in the PDR; it includes both a graphitic and a silicate component with R<sub>V</sub> ≡ A<sub>V</sub>/E(B − V) = 3.1. The grain density scales with the hydrogen density n<sub>H</sub>. The line fluxes are re-scaled with the number of Lyman continuum photons, extracted from the stellar emission of the modelled galaxies. These line models will be detailed in Theulé et al. (in prep.).

For the dust emission, we selected the modified black-
body module in CIGALE, with a power law in the mid-IR

<footnote n="2">² frac<sub>AGN</sub>, the AGN fraction, is defined as L<sub>AGN</sub> / (L<sub>AGN</sub> + L<sub>dust</sub>).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** How the dust continuum was taken into account in the analysis. Top panel shows the distribution of $T_{\text{dust}}$ and $\beta_{\text{RL}}$ for the SPT sample. In the bottom panel, we used the CIGALE mock analysis to check whether $T_{\text{dust}}$ and $\beta_{\text{RL}}$ can be correctly estimated using the available set of data for this SPT sample. In this mock analysis, we use the best models and each of the exact parameters used to build the models for each of the objects. We add the observed noise to the models, and refit the simulated data to re-estimate the same dust and $\beta_{\text{RL}}$ parameters. A good correlation between the 'estimated' and the 'exact' $T_{\text{dust}}$ and $\beta_{\text{RL}}$ suggests that we are able to derive them correctly.</figcaption></figure>

(Casey 2012). In Fig. 2, we present the dust temperature and the Rayleigh-Jeans slope ($T_{\text{dust}}$ and $\beta_{\text{RL}}$) controlling the shape of the IR continuum SEDs that are derived by fitting the observed SEDs with CIGALE. $T_{\text{dust}}$ and $\beta_{\text{RL}}$ derived by fitting the data are only meant to check their measured range. They are not used in the rest of the analysis. The CIGALE mock analysis, also presented in Fig. 2, suggests that $T_{\text{dust}}$ and $\beta_{\text{RL}}$ can be well estimated, with a coefficient of correlation of $r^2 = 0.92$ for $T_{\text{dust}}$ and $r^2 = 0.80$ for $\beta_{\text{RL}}$.

In Fig. 3, the modelled colours are compared to the observed ones in the colour-colour diagram with and without adding emission lines to the modelled continuum. As above, we assume a mid-IR power law and a modified blackbody with a dense sampling in dust temperature, $T_{\text{dust}}$, and emissivity, $\beta_{\text{RL}}$, to best reproduce the rest-frame far-IR dust continuum emission.

The main apparent trend observed in the colour-colour diagram without emission lines, is a redshift-related sequence from the bottom right to the top left of the plot. This prime sequence is due to the peak of the IR dust emission, passing in the broad bands (see e.g., Amblard et al. 2010). When no fine structure

emission lines are added in CIGALE models, all the models are located in, or just below this prime sequence. However, we can see in Fig. 3 that a clump of objects is located above the prime sequence, at $[\log_{10} \frac{\text{PSW}_{250\mu\text{m}}}{\text{PMW}_{350\mu\text{m}}}, \log_{10} \frac{\text{LABOCA}_{870\mu\text{m}}}{\text{PLW}_{500\mu\text{m}}}] \approx [-0.1, -0.1]$, with a total offset from the main sequence of the order of $\Delta_{\text{colour}} \sim 0.1-0.2$. This clump contains most of the highest redshift galaxies at $z \ge 4-6$. They cannot be explained by changes in the dust continuum only. In Appendix A, we quantitatively explore the possibility that emission lines are the most important piece of this puzzle. We find that galaxies exhibiting large [CII]157.6 µm equivalent width (EW) reaching EW([CII]157.6 µm) ~ 10–20 µm could explain such outliers. This is in agreement with the estimates from Smail et al. (2011), which found that the [CII]158 µm line could provide as much as 40% of the 850 µm broadband flux density for $z \sim 4-6$ galaxies, when $L_{[CII]}/L_{\text{dust}} > 1\%$.

Figure 4 compares the evolution in redshift of a model with emission lines, when the PDR dominates the nebular emission ($\log_{10} U = -4.0$), and that of a model also with emission
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Models computed by CIGALE for the $[\log_{10}(\frac{\text{PSW}_{250\mu m}}{\text{PMW}_{350\mu m}}) \text{ versus } \log_{10}(\frac{\text{LABOCA}_{870\mu m}}{\text{PLW}_{500\mu m}})]$ colours. Left: when no fine structure emission lines are added to the modelled spectra, the modelled galaxies are mainly located on a prime sequence extending from the bottom right to the top left. However, a clump of high-redshift observed galaxies, with colours $[\log_{10}(\frac{\text{PSW}_{250\mu m}}{\text{PMW}_{350\mu m}}) \text{ versus } \log_{10}(\frac{\text{LABOCA}_{870\mu m}}{\text{PLW}_{500\mu m}})] \approx [-0.1, -0.1]$, and with a total offset from the main sequence of the order of $\Delta_{\text{colour}} \sim 0.1-0.2$ cannot be reached by these models without emission lines. Right: emission lines are added to the continuum. When large equivalent width fine structure lines (and more specifically, EW([CII]157.6 µm) ~ 10–20 µm) are added (Appendix A), the CIGALE models cover this clump of high-redshift galaxies. The three highest points are only partially covered by the models. Even though the two top ones (namely SPT0243-49 and SPT0245-63) are clearly high-redshift objects, respectively at $z = 5.702$ and $z = 5.626$, we assumed that the partial coverage is not sufficient and do not use these two objects in the rest of the paper. For the third one (SPT0348-62 at $z = 5.654$), the coverage by the models allows us to keep it in the analysis. We also identify SPT0418-47, for which we have five emission lines in the considered wavelength range, and an amplification factor of $\mu = 32.70$ (Reuter et al. 2020).</figcaption></figure>

lines, but when the HII regions dominate the nebular emission ($\log_{10} U = -2.0$). Both are also compared to the same models (i.e. for the same dust temperature, $T_{\text{dust}}$, and dust emissivity on the Rayleigh-Jeans side, $\beta_{\text{RL}}$) without lines. The [CII]157.6 µm has a strong impact on the colours of PDR-dominated models. The models with large EWs ([CII]157.6 µm) entering the LABOCA$_{870\mu m}$ band at $z \sim 4-6$ can explain the clump of galaxies above the prime sequence. At lower redshifts and for HII models, the combination of the [OIII]51.7 µm, and [OIII]88.3 µm lines (and others at lower levels) also impacts on the broad band colours, and offsets the models to below the prime sequence on the $\log_{10}(\frac{\text{LABOCA}_{870\mu m}}{\text{PLW}_{500\mu m}})$ axis. However, some of the models without lines overlap with this region (Fig. 3), which makes the identification of galaxies in this redshift range less conclusive for galaxies dominated by HII regions.

## 4. A colour-colour approach to estimate the nebular parameters of SPT galaxies

We saw the impact of emission lines on the colour distributions. We now try to constrain the physical parameters driving the intensity of the emission lines. For each of the SPT objects, we compute the mean of the models that lie inside ellipses delimited by the uncertainties in both colours. We stress that for each object we only keep the models within $\delta z = \pm 0.1$ from the spectroscopic redshift. This provides us with the mean and standard deviation of the physical parameters used to compute the models, as well as those for the line and continuum fluxes. As already mentioned, the top two objects are very partially covered by the

models in Fig. 3: SPT0243-49 and SPT0245-63 at $z = 5.702$ and $z = 5.626$, respectively. They are not used hereafter.

We compare the modelled and observed colours in Fig. 5. The regressions provide correlation coefficient of 0.997 for the LABOCA$_{870\mu m}$/PLW colour and 0.957 for the PSW/PMW colour. The median and standard deviations of the modelled-to-observed colour differences ($\Delta[\log_{10}(\text{Colour}_{\text{Modelled}}) - \log_{10}(\text{Colour}_{\text{Observed}})] = 0.001 \pm 0.036$ for the PSW/PMW colour and $0.003 \pm 0.011$ for the LABOCA$_{870\mu m}$/PLW colour.

In order to check our results, we compare the observed emission line fluxes of the SPT sample to our modelled ones in Fig. 6. However, this comparison is limited because only a small sample of SPT galaxies has been spectroscopically observed, and most of them only with one line.

For a sub-sample of the SPT galaxies, Lagache et al. (2018) agreed that [CII]157.6 µm at high redshift mainly originates from the PDR. For the objects in common with the present sample, their [CII]157.6 µm luminosities are compared to our derived [CII]157.6 µm luminosities. Cunningham et al. (2020) showed that 57% of the SPT sample presents an $L_{[CII]157.6\mu m}/L_{[NII]205.2\mu m}$ luminosity ratio (or lower limit) in agreement with those expected from PDR (or shock regions). However, they suggest that a sub-sample (~27%) of the $3 < z < 6$ SPT galaxies would be consistent (within uncertainties) with a hybrid regime between the model predictions of PDR emission and HII regions. For the objects in common with the present sample, their [NII]205.2 µm luminosities are compared to our derived [NII]205.2 µm luminosities.
<figure><img src="image_6.png" /></figure>

**Fig. 4.** Evolution of the predicted colours with redshift. When compared to models without lines (crosses), at $log_{10}U = -4.0$ (upward triangles), the strong [CII]157.6 µm line induces an upward move of $log_{10}(\frac{\text{LABOCA}_{870\mu m}}{\text{PLW}_{500\mu m}})$ that corresponds to the clump of galaxies at $z \sim 4-6$. However, even though the move of the [OIII]51.8 µm and [OIII]88.3 µm lines of galaxies with a strong emission from HII regions at $log_{10}U = -2.0$ (downward triangles) could induce specific colours, the effect is less clear as models with no lines can also lie in this same region of the plot (Fig. 3). All the symbols are colour-coded by redshift.

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

**Fig. 5.** Comparison of compared and observed colours. This figure shows that we can reproduce the observed colours with the CIGALE models. Top: modelled colours (green boxes for $log_{10}(\frac{\text{LABOCA}_{870\mu m}}{\text{PLW}_{500\mu m}})$ and red dots for $log_{10}(\frac{\text{PSW}_{250\mu m}}{\text{PMW}_{350\mu m}})$) are in excellent agreement with the observed ones within the observed uncertainties. The black line shows the one-to-one line, and the grey shaded area presents the boundaries at $\pm 5\%$ ($y = 0.95x$ and $y = 1.05x$). Bottom: This conclusion is confirmed in this panel, where we present the absolute differences in the logarithm of the two colours used in Fig. 3. In both panels, the red circled point corresponds to SPT0348-62, which is already identified in Fig. 3 as one of the objects only partially covered by the models.

Statistical tests are performed with the Python LINMIX library (Fig. 6). This LINMIX method presents the advantage of using a hierarchical Bayesian approach for the linear regression, which takes errors in both X and Y into account (Kelly

2007). The tests are significant (see Fig. 6), given the number of points and the linear correlation coefficient between the observed and modelled fluxes: $r = 0.62$. We check whether the different locations of the lines, and most notably [CII]158 µm and [NII]205 µm for the present sample, are due to the physical mechanisms producing these two lines. However, we could not pinpoint any differences in the critical densities or ionisation potentials that could explain the differences in the figure (see e.g., Fig. 2 in Spinoglio et al. 2015). We tentatively conclude that the main reason for the increased distance to the one-to-one line is very likely the difference in line intensity for [CII]158 µm and [NII]205 µm. In other words, this method is certainly more sensitive to strong lines (CII and OIII) than to faint lines (NII). The fact that there is about one order of magnitude offset between the observed and modelled line fluxes, for the weak lines, is a limitation for the method presented in this paper because the nebular parameters are based on line ratios. We note, however, another horizontal structure that does not seem random. The redshift does not explain this structure. The galaxies with an intrinsic (i.e. corrected for the amplification using Reuter et al. 2020) dust luminosity ($log_{10}L_{\text{dust}}$) in the $12.7 < log_{10}L_{\text{dust}} < 13.5$ range are generally found at larger distances, above the one-to-one line. On the contrary, galaxies with $log L_{\text{dust}} < 12.7$ and galaxies with $log L_{\text{dust}} > 13.5$ are significantly closer to this one-to-one line. No clear physical origin is identified for this differential effect, though. It could be related to the evolution of these very exotic high redshift objects. More data, and especially rest-frame UV (from JWST) and far-IR (from ALMA or NOEMA) morphologies are fundamental clues to deciphering the structure of this diagram.

From the present SPT sample, only for SPT0418-47 do we have several emission lines (identified with large open circles in the right panel of Fig. 6) that allow us to check how well we model the line ratios. Three of the $L_{\text{line}}/L_{\text{[OIII]88.3\mu m}}$ modelled ratios are in agreement, within a factor of three at most, to the observed ones ($L_{\text{[NII]122\mu m}}/L_{\text{[OIII]88.3\mu m}} = 0.048 \pm 0.014$, $L_{\text{[CII]158\mu m}}/L_{\text{[OIII]88.3\mu m}} = 0.667 \pm 0.091$,
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Comparison of the CIGALE-derived line-to-dust luminosities to the observed ones. Left: comparison of observations and models for the spectroscopic sample of SPT galaxies. Most data are from [CII]157.6 µm (Lagache et al. 2018) and [NII]205.2 µm (Cunningham et al. 2020), except for SPT 0418-47 (De Breuck et al. 2019), for which we have several lines. Generally speaking, the brightest [CII]157.6 µm and [OIII]88.3 µm lines are better modelled than the [NII]122. µm and [NII]205.2 µm. The [OI]145. µm is an outlier in this frame. It could be because the former lines are stronger than the latter, or because of physical differences between them. Right: Symbols are colour-coded by intrinsic dust luminosity. The fits suggest that the objects with $log_{10} L_{dust} < 12.7$ and $log_{10} L_{dust} > 13.5$ are found closer to the one-to-one line, or vice versa, meaning that the line predictions for those with $12.7 < log_{10} L_{dust} < 13.5$ are worse. The objects with a large green circle correspond to the sole object with different lines: SPT0418-47 from De Breuck et al. (2019) (see left panel to identify the lines).</figcaption></figure>

and $L_{[NII]205\mu m}/L_{[OIII]88.3\mu m} = 0.032 \pm 0.004$), while our estimates are 0.143, 0.619, and 0.085, respectively. The observed $L_{[OI]145\mu m}/L_{[OIII]88.3\mu m} = 0.189 \pm 0.050$, while we find 0.011. This is a factor of almost 17 for [OI]145 µm, which is significantly larger than the other line ratios. However, De Breuck et al. (2019) also had the same problem: their predicted range for this line is $0.02-1.9 \times 10^{10} L_{\odot}$, while their detection amounts to $2.1 \times 10^{10} L_{\odot}$ above the maximum predicted value. The conclusion again is that we need to collect more objects with these line ratios to be able to use a statistical approach and clarify the situation. For SPT0418-47, De Breuck et al. (2019) derived a gas metallicity of $0.3 < Z/Z_{\odot} < 1.3$ and an ionisation parameter of $-3.2 < \log_{10} U < -2.0$. These values are in reasonable agreement (within the uncertainties) with ours: $Z/Z_{\odot} \sim 1.47 \pm 1.14$ and $\log_{10} U \sim -3.02 \pm 0.80$.

We compute the mean of the model parameters inside ellipses delimited by the observational uncertainties. The metallicities are in the $0.5 \le Z/Z_{\odot} \le 2.5$ range and the ionisation parameters cover from $-4.0 \le \log_{10} U \le -1.5$. (Fig. 7). Very high metallicity objects are rare but observed in the local Universe (Gallazzi et al. 2005; Peeples et al. 2008; Maiolino & Mannucci 2019) with metallicities that extend to $>3Z_{\odot}$. Such high values for the ionisation parameter were also measured in the overlap region of the Antennae galaxies by, for example, Snijders et al. (2007), Kewley et al. (2019). Yeh & Matzner (2012) found that radiation pressure confinement sets an upper limit to $\log_{10} U = -1$ in individual regions. However, when observing unresolved starbursts, the mean values are of the order of $\log_{10} U = -2.3$ due to the variety of regions inside a galaxy. As noted earlier in this paper, the two top most objects in Fig. 3 are either very extreme cases among the diversity of galaxies, or the observed colour uncertainties are underestimated. We do not keep them the analysis. However, they probably deserve targeted studies.

**5. A structure in the log₁₀ U versus log₁₀ Z/Zₐ diagram**

From the parameters derived in the previous section, we build the $\log_{10} U$ versus $Z/Z_{\odot}$ diagram (Fig. 7), where we identify two sequences: the top one with a strong $L_{[OIII]88.3\mu m}$ and the bottom one with a strong $L_{[CII]157.6\mu m}$ emission. Even though the models are created on a dense regular grid, with a flat probability (bottom panel of Fig. 7), we stress that the metallicities and ionisation parameters derived from the method presented in this paper cannot be precise enough to define the two sequences as clearly as they appear in Fig. 7. The well-defined sequences might be due to the fact that the mean values are estimated from wide probability distribution functions, as confirmed by the uncertainties shown in Fig. 7. The mean values of these wide distributions regularly evolve in the plot, that could suggest the impression of well-defined sequences.

The bottom sequence presents $\log_{10} U$ values that are similar to PDR-dominated galaxies, while they are more similar to HII region-dominated galaxies for the top sequence. This type of objects would be galaxies that present a so-called [CII] deficit. These [CII]-deficit galaxies show a ratio of $L_{[OIII]88.3\mu m}/L_{[CII]157.6\mu m} \approx 3-20$ that is about ten times higher than $z \sim 0$ galaxies. Harikane et al. (2020) identified nine $z = 6-9$ galaxies whose observed properties are in agreement with being such [CII] deficit galaxies. Numerous explanations have been proposed: differences in C and O abundance ratios, observational biases, and differences in ISM properties. Carniani et al. (2020) suggested that a surface brightness dimming of the extended [CII] emission would be responsible for the [CII] deficit. Harikane et al. (2020) explained these high $L_{[OIII]88.3\mu m}/L_{[CII]157.6\mu m}$ ratios by high ionisation parameters or low PDR covering fractions, both of which are consistent with their [NII] observations. This scenario could be reproduced
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Sequence where we find LBGs at $z \sim 3.3$ sample (Onodera et al. 2016) and fit (dashed line) to the sequence from Kojima et al. (2017) are shown in grey. Two different branches are identified for which the [CII]157.6 µm (centre) and [OIII]88.3 µm (top) lines are strong. The top branch contains galaxies that present a [CII] deficit and have a stronger ionisation parameter, $log_{10} U \lesssim -2.5$ while the bottom branch corresponds to galaxies with a normal-to-extreme PDR emission. The bottom panel shows the same information colour-coded with the electronic density, with the grid of models superimposed showing that we should be able to derive metallicities and ionisation parameters in between the two branches.</figcaption></figure>

by a density-bound nebula with a PDR deficit. In radiation-hydrodynamics simulations, Abel et al. (2009) concluded that the effects of high ratios of impinging ionising radiation density

to particle density (i.e. again, high ionisation parameters) can reproduce the observational characteristics of ultra luminous IR galaxies (ULIRG). When $U$ increases, the fraction of UV photons absorbed by dust increases, and fewer photons are available to photoionise and heat the gas. This leads to a dust-bound nebula, which can explain the observed [CII] deficit (see also Fischer et al. 2014, for a slightly more complex but consistent explanation).

## 6. Conclusions

We performed an analysis of the SEDs of SPT galaxies with spectroscopic redshifts. We built more than $88 \times 10^6$ models with CIGALE and compared the observed objects to the models in a [$log_{10}$ (PSW$_{250\mu m}$/PMW$_{350\mu m}$) versus $log_{10}$(LABOCA$_{870\mu m}$/PLW$_{500\mu m}$)] colour-colour diagram. This set of colours was selected to identify galaxies at $z \gtrsim 4-6$ from the influence of the [CII]158 µm fine structure lines on broad bands. This method also allows us to roughly estimate the nebular parameters for this SPT sample.

From this analysis, we find the following results.

- The position of the SPT $z \sim 4-6$ galaxies in the colour-colour diagram can only be explained when adding the contribution of fine structure far-IR emission lines to the dust continuum.

- By averaging the models and their associated physical parameters in ellipses delimited by the observed uncertainties in both colours, we can estimate the flux of the fine structure far-IR emission lines, the gas metallicity ($Z_{gas}$), and the ionisation parameter ($log_{10} U$). We find that all the SPT galaxies have a high gas metallicity of $0.6 < Z_{gas} < 2.5$ and they cover a wide range of ionisation parameters in the $log_{10} U = -4.0$ to -1.5 range. However, we add a word of caution because the faintest lines are overestimated. Thus, line ratios involving faint lines could bias the estimation of nebular parameters.

- In the $log_{10} U$ versus $Z_{gas}$ diagram, we identify two branches with high [OIII]88 µm/[CII]158 µm ratios for the top branch and low [OIII]88 µm/[CII]158 µm ratios for the bottom branch. The top branch, with high $log_{10} U \lesssim -2.5$, presents a deficit in $L_{CII}/L_{dust}$ with respect to the bulk of the galaxy sample, while the bottom branch is the extension of a sequence that continues to LBGs at lower metallicities.

- Without emission lines, outliers are offset from the prime sequence by $|\Delta(log_{10}(PSW_{250\mu m}/PMW_{350\mu m}))| \approx 0.05$ and $|\Delta(log_{10}(LABOCA_{870\mu m}/PLW_{500\mu m}))| \approx 0.09$. That is a total offset of the order of $\Delta_{total} \sim 0.10$. For SPT galaxies at $z \gtrsim 4-6$, the main effect is due to [CII]158 µm with $0.6 \lesssim EW([CII]157.6\mu m) \lesssim 25.0\mu m$. For the most extreme cases, this line could be at the origin of almost half of the flux density at 850 µm, for galaxies at $z \sim 4-6$.

- In order to make the most efficient use of this method, a set of medium and broad bands in the mid- and far-IR would be ideal because the effect of the lines would be stronger in narrower bands. One of the caveats of this method is the need to collect an SED as complete as possible to correctly estimate the line fluxes: the more complete the SED, the better the line fluxes can be estimated. Thus, to be efficient, the utility of the method presented here relies on large photometric samples, which are cheaper to obtain than spectroscopy. In other words, our method would benefit from large photometric surveys on a large galaxy sample. Otherwise, spectroscopic observations present the advantage of providing much better estimates. With this configuration, a project like the PRobe far-Infrared Mission for Astrophysics
(PRIMA<fnref n="3" />) allows a statistical approach that would permit us to understand the cosmic rise of metals up to the reionisation. The results found in this paper will be useful for two reasons: the first one is the easier identification of high-redshift galaxies by making use of well-thought colours (as a function of the redshift) in where strong emission lines will produce a measurable excess. The second one by providing a rough but statistical information on nebular parameters for large photometric samples.

*Acknowledgements.* The authors thank Charles M. (Matt) Bradford for a very useful discussion. DB acknowledges support from the Centre National d'Etudes Spatiales (CNES) for this effort of simulating galaxies and their properties derived from space facilities like Herschel and predict new space observations. Médéric Boquien gratefully acknowledges support by the ANID BASAL project FB210003 and from the FONDECYT regular grant 1211000. A.K.I., Y.S., and Y.F. are supported by NAOJ ALMA Scientific Research grant No. 2020-16B.

## References

Abel, N. P., Dudley, C., Fischer, J., Satyapal, S., & van Hoof, P. A. M. 2009, ApJ, 701, 1147

Algera, H. S. B., Inami, H., Oesch, P. A., et al. 2023, MNRAS, 518, 6142

Amblard, A., Cooray, A., Serra, P., et al. 2010, A&A, 518, L9

Anders, P., & Fritze-v. Alvensleben, U., 2003, A&A, 401, 1063

Apostolovski, Y., Aravena, M., Anguita, T., et al. 2019, A&A, 628, A23

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Béthermin, M., De Breuck, C., Gullberg, B., et al. 2016, A&A, 586, L7

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413

Carniani, S., Ferrara, A., Maiolino, R., et al. 2020, MNRAS, 499, 5136

Casey, C. M. 2012, MNRAS, 425, 3094

Chabrier, G. 2003, PASP, 115, 763

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Cunningham, D. J. M., Chapman, S. C., Aravena, M., et al. 2020, MNRAS, 494, 4090

de Barros, S., Nayyeri, H., Reddy, N., & Mobasher, B. 2013, in SF2A-2013: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. L. Cambresy, F. Martins, E. Nuss, & A. Palacios, 531

De Breuck, C., Weiß, A., Béthermin, M., et al. 2019, A&A, 631, A167

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. A&A, 53, 385

Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2022, ApJ, 940, L55

Fischer, J., Abel, N. P., González-Alfonso, E., et al. 2014, ApJ, 795, 117

Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41

Gururajan, G., Béthermin, M., Theulé, P., et al. 2022, A&A, 663, A22

Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93

Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2018, ApJ, 861, 95

Kelly, B. C. 2007, ApJ, 665, 1489

Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511

Kojima, T., Ouchi, M., Nakajima, K., et al. 2017, PASJ, 69, 44

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Lagache, G., Cousin, M., & Chatzikos, M. 2018, A&A, 609, A130

Luhman, M. L., Satyapal, S., Fischer, J., et al. 2003, ApJ, 594, 758

Ma, J., Gonzalez, A. H., Vieira, J. D., et al. 2016, ApJ, 832, 114

Maiolino, R., & Mannucci, F. 2019, A&A Rev., 27, 3

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Onodera, M., Carollo, C. M., Lilly, S., et al. 2016, ApJ, 822, 42

Peeples, M. S., Pogge, R. W., & Stanek, K. Z. 2008, ApJ, 685, 904

Reuter, C., Vieira, J. D., Spilker, J. S., et al. 2020, ApJ, 902, 78

Roberts-Borsani, G. W., Ellis, R. S., & Laporte, N. 2020, MNRAS, 497, 3440

Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187

Schärer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4

Seaquist, E., Yao, L., Dunne, L., & Cameron, H. 2004, MNRAS, 349, 1428

Seymour, N., Altieri, B., De Breuck, C., et al. 2012, ApJ, 755, 146

Smail, I., Swinbank, A. M., Ivison, R. J., & Ibar, E. 2011, MNRAS, 414, L95

Snijders, L., Kewley, L. J., & van der Werf, P. P. 2007, ApJ, 669, 269

Spilker, J. S., Aravena, M., Béthermin, M., et al. 2018, Science, 361, 1016

Spinoglio, L., Pereira-Santaella, M., Dasyra, K. M., et al. 2015, ApJ, 799, 21

Strandet, M. L., Weiss, A., Vieira, J. D., et al. 2016, ApJ, 822, 80

Taylor, A. J., Barger, A. J., & Cowie, L. L. 2022, ApJ, 939, L3

Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2022, AAS J., submitted, [arXiv:2207.12388]

Yeh, S. C. C., & Matzner, C. D. 2012, ApJ, 757, 108

<footnote n="3">³ See https://prima.ipac.caltech.edu/ and https://agora.lam.fr/</footnote>
# Appendix A: Contribution of emission lines to broad bands

In this appendix, we quantitatively estimate the contribution of the far-IR emission lines to the broad bands. We also show that only the emission lines with the largest equivalent widths, namely [OIII]88.3 µm and [CII]157.6 µm for the considered broad bands and redshift range, could have an impact on the observed colours.

The offset in the colour-colour diagram (Fig. 3) could be due to changes of the PSW<sub>250µm</sub> / PMW<sub>350µm</sub> and/or the LABOCA<sub>870µm</sub> / PLW<sub>500µm</sub> colours. Even though this is an over simplification (Fig. A.1), for the sake of clarity, we assume that the main effect is due to [CII]157.6 µm entering or exiting the LABOCA<sub>870µm</sub>, and nothing else is modified. The flux density in band PLW<sub>500µm</sub> does not change over the redshift range $4 \le z \le 6$. In this case, we have

$$F_{\text{PLW}_{500\mu\text{m}}}^{\text{line}} = F_{\text{PLW}_{500\mu\text{m}}}^{\text{noline}} \& F_{\text{LABOCA}_{870\mu\text{m}}}^{\text{line}} \ge F_{\text{PLW}_{500\mu\text{m}}}^{\text{noline}} \quad (\text{A.1})$$

We define the band ratio as follows:

$$R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{line}} = \frac{F_{\text{LABOCA}_{870\mu\text{m}}}^{\text{line}} \Delta\lambda_{870\mu\text{m}}}{F_{\text{PLW}_{500\mu\text{m}}}^{\text{line}} \Delta\lambda_{500\mu\text{m}}}, \quad (\text{A.2})$$

where $\Delta\lambda_{500\mu\text{m}} = 143$ µm, $\Delta\lambda_{500\mu\text{m}} = 186$ µm, $\Delta\lambda_{250\mu\text{m}} = 67$ µm, and $\Delta\lambda_{350\mu\text{m}} = 95$ µm (SVO Filter Service). Using the right end of Eq. C.1, we have

$$R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{line}} \ge R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{noline}} \quad (\text{A.3})$$

From Eq. C.2, where $F_{line}$ is the line flux contributing to the LABOCA<sub>870µm</sub> band,

$$R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{line}} = \frac{F_{\text{LABOCA}_{870\mu\text{m}}}^{\text{noline}} \Delta\lambda_{870\mu\text{m}} + F_{\text{line}}}{F_{\text{PLW}_{500\mu\text{m}}}^{\text{noline}} \Delta\lambda_{500\mu\text{m}}} \quad (\text{A.4})$$

From the left end of Eq. C.1, and the definition of the equivalent width ($EW_{line}$), we have $F_{line} = F_{\text{LABOCA}_{870\mu\text{m}}}^{\text{noline}} \Delta\lambda_{870\mu\text{m}}$. This gives:

$$R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{line}} = \frac{F_{\text{LABOCA}_{870\mu\text{m}}}^{\text{noline}}}{F_{\text{PLW}_{500\mu\text{m}}}^{\text{noline}}} (\Delta\lambda_{870\mu\text{m}} + EW_{\text{line}}) \quad (\text{A.5})$$

and, finally,

$$R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{line}} = R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{noline}} \left(1 + \frac{EW_{\text{line}}}{\Delta\lambda_{870\mu\text{m}}}\right). \quad (\text{A.6})$$

The mean colours of the outliers are $\langle\log_{10}(\text{PSW}_{250\mu\text{m}} / \text{PMW}_{350\mu\text{m}})\rangle = -0.205 \pm 0.113$ and $\langle\log_{10}(\text{LABOCA}_{870\mu\text{m}} / \text{PLW}_{500\mu\text{m}})\rangle = -0.187 \pm 0.113$. For the same redshift range, CIGALE models without lines give $\langle\log_{10}(\text{PSW}_{250\mu\text{m}} / \text{PMW}_{350\mu\text{m}})\rangle = -0.254 \pm 0.134$ and $\langle\log_{10}(\text{LABOCA}_{870\mu\text{m}} / \text{PLW}_{500\mu\text{m}})\rangle = -0.101 \pm 0.194$. The average move needed to reach these galaxies amounts to $|\Delta(\log_{10}(\text{PSW}_{250\mu\text{m}} / \text{PMW}_{350\mu\text{m}}))| \approx 0.05$ and $|\Delta(\log_{10}(\text{LABOCA}_{870\mu\text{m}} / \text{PLW}_{500\mu\text{m}}))| \approx 0.09$. That is a total offset of the order of $\Delta_{\text{total}} \sim 0.10$. Because we cannot know where the galaxy would be without accounting for the emission lines, both colours could contribute to this offset, but it is likely that [CII]157.6 µm is dominant.

We find that the models inside the ellipses, accounting for the observed uncertainties, have $0.6 \le EW([CII]157.6 \mu\text{m}) \le 25.0 \mu\text{m}$. From this, we obtain

$$1.0 \le \frac{R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{line}}}{R_{870\mu\text{m}_{500\mu\text{m}}}^{\text{noline}}} \le 1.2, \quad (\text{A.7})$$

$$0.0 \le \Delta(\log_{10}(\text{LABOCA}_{870\mu\text{m}} / \text{PLW}_{500\mu\text{m}})) \le 0.1, \quad (\text{A.8})$$

which is an offset in the colour-colour diagram of $\le 0.1$, that is, about the same order or the maximum one observed.

The [SIII]33.47 µm line enters the PSW<sub>250µm</sub> at about the same redshift range and modifies $\log_{10}(\text{PSW}_{250\mu\text{m}} / \text{PMW}_{350\mu\text{m}})$. However, with $0.05 \le EW([SIII]33.5 \mu\text{m}) \le 0.26$, we estimate that the [SII]33.47 µm line should not significantly contribute to the PSW<sub>250µm</sub> band.

No CIGALE models at $\delta z = \pm 0.5$ match the two highest objects in Fig. 3 (SPT0245-63 at $z = 5.626$ and SPT0243-49 at $z = 5.702$), and we cannot estimate any of the physical parameters. The reasons why we cannot reproduce the colours of these two objects is uncertain. Two plausible hypotheses could be made, though. First, the flux densities of these objects might not be correct, or, at least, the uncertainties could be underestimated. The other explanation might be that our present grid of models might not cover the entire possible range of data. These two objects deserve a closer analysis.
<figure><img src="image_10.png" /></figure>

**Fig. A.1.** This figure identifies the emission lines that produce excesses in the broad-bands when the redshift increases. Top: At $log_{10}$ U = -4.0, the strong [CII]157.6 µm lines enter the LABOCA$_{870\mu m}$ filters at z ~ 4 and exit at z ~ 5.2. The line induces an upward move of the corresponding $log_{10} \frac{\text{LABOCA}_{870\mu m}}{\text{PLW}_{500\mu m}}$ colour. [SIII]34.8 µm enters the PSW$_{250\mu m}$ at z ~ 5.5, which boosts the $log_{10} \frac{\text{PSW}_{250\mu m}}{\text{PMW}_{350\mu m}}$ colour. We find a total move of the galaxies by about $\Delta_{total} \sim 0.10$ in the colour-colour diagram. These two combined effects explain that PDR-dominated galaxies at 4 $\le$ z $\le$ 6 are offset from the other galaxies. Bottom: Movement of the [OIII]51.8 µm and [OIII]88.3 µm lines of galaxies with a strong emission from HII regions at $log_{10}$ U = -2.0 could induce specific colours that could help identify such galaxies from this colour-colour diagram. However, the effect is less clear as models with no lines can also lie here (Fig. 3).