# The ALPINE-CRISTAL-JWST survey: Spatially resolved star formation relations at z ~ 5

C. Accard, M. Béthermin, M. Boquien, V. Buat, L. Vallini, F. Renaud, K.

Kraljic, M. Aravena, P. Cassata, E. da Cunha, et al.

► To cite this version:

C. Accard, M. Béthermin, M. Boquien, V. Buat, L. Vallini, et al. The ALPINE-CRISTAL-JWST survey: Spatially resolved star formation relations at z ~ 5. Astronomy & Astrophysics - A&A, 2025, 702, pp.A206. 10.1051/0004-6361/202556140. hal-05331056

HAL Id: hal-05331056

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

Submitted on 25 Oct 2025

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
# The ALPINE-CRISTAL-JWST survey: Spatially resolved star formation relations at z ~ 5

C. Accard<sup>1,*</sup> <u>D</u>, M. Béthermin<sup>1</sup> <u>D</u>, M. Boquien<sup>2</sup> <u>D</u>, V. Buat<sup>3</sup>, L. Vallini<sup>4</sup> <u>D</u>, F. Renaud<sup>1</sup> <u>D</u>, K. Kraljic<sup>1</sup> <u>D</u>, M. Aravena<sup>5,6</sup>, P. Cassata<sup>8,9</sup> <u>D</u>, E. da Cunha<sup>10,11</sup>, P. Dam<sup>8</sup>, I. de Looze<sup>12</sup> <u>D</u>, M. Dessauges-Zavadsky<sup>13</sup> <u>D</u>, Y. Dubois<sup>14</sup>, A. Faisst<sup>15</sup> <u>D</u>, Y. Fudamoto<sup>16</sup>, M. Ginolfi<sup>17,18</sup>, C. Gruppioni<sup>4</sup> <u>D</u>, S. Han<sup>14</sup> <u>D</u>, R. Herrera-Camus<sup>6,7</sup>, H. Inami<sup>19</sup> <u>D</u>, A. M. Koekemoer<sup>20</sup> <u>D</u>, B. C. Lemaux<sup>21,22</sup>, J. Li<sup>10,11</sup> <u>D</u>, Y. Li<sup>23</sup> <u>D</u>, B. Mobasher<sup>24</sup> <u>D</u>, J. Molina<sup>25,6</sup> <u>D</u>, A. Nanni<sup>26,27</sup>, M. Palla<sup>4,28</sup> <u>D</u>, F. Pozzi<sup>4</sup> <u>D</u>, M. Relaño<sup>29</sup> <u>D</u>, M. Romano<sup>30,9</sup> <u>D</u>, P. Sawant<sup>27</sup> <u>D</u>, J. Spilker<sup>23</sup>, A. Tsujita<sup>34</sup> <u>D</u>, E. Veraldi<sup>31</sup> <u>D</u>, V. Villanueva<sup>32</sup> <u>D</u>, W. Wang<sup>15</sup> <u>D</u>, S. K. Yi<sup>33</sup>, and G. Zamorani<sup>4</sup>

*(Affiliations can be found after the references)*

Received 27 June 2025 / Accepted 15 August 2025

## ABSTRACT

**Context.** Star formation governs galaxy evolution, shaping stellar mass assembly and gas consumption across cosmic time. The Kennicutt-Schmidt (KS) relation, linking the star formation rate (SFR) and gas surface densities, is fundamental to understand star formation regulation, yet remains poorly constrained at $z > 2$ due to observational limitations and uncertainties in locally calibrated gas tracers. The [CII] 158 µm line has recently emerged as a key probe of the cold ISM and star formation in the early Universe.

**Aims.** We investigate whether the resolved [CII]–SFR and KS relations established at low redshift remain valid at $4 < z < 6$ by analysing 13 main-sequence galaxies from the ALPINE and CRISTAL surveys, using multi-wavelength data (HST, JWST, ALMA) at ~2 kpc resolution.

**Methods.** We performed pixel-by-pixel spectral energy distribution (SED) modelling with CIGALE on resolution-homogenised images. We developed a statistical framework to fit the [CII]–SFR relation that accounts for pixel covariance and compare our results to classical fitting methods. We tested two [CII]-to-gas conversion prescriptions to assess their impact on inferred gas surface densities and depletion times.

**Results.** We find a resolved [CII]–SFR relation with a slope of $0.87 \pm 0.15$ and intrinsic scatter of $0.19 \pm 0.03$ dex, which is shallower and tighter than previous studies at $z \sim 5$. The resolved KS relation is highly sensitive to the [CII]-to-gas conversion factor: using a fixed global $\alpha_{[CII]}$ yields depletion times of 0.5–1 Gyr, while a surface brightness-dependent $W_{[CII]}$, accounting for local ISM conditions, places some galaxies with high gas density in the starburst regime (<0.1 Gyr).

**Conclusions.** Future inputs from both simulations and observations are required to better understand how the [CII]-to-gas conversion factor depends on local ISM properties. We need to break this fundamental limit to properly study the KS relation at $z \ge 4$.

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

## 1. Introduction

The study of star formation in the early Universe ($z > 2$) is a cornerstone of modern astrophysics, crucial for understanding galaxy evolution, mass assembly, and cosmic star formation history (SFH; e.g. Madau & Dickinson 2014; Freundlich 2024). The Kennicutt-Schmidt (KS) relation, linking star formation rate (SFR) surface density ($\Sigma_{SFR}$) to gas surface density ($\Sigma_{gas}$) through a power law, has been fundamental for understanding the conversion of the gas into stars across cosmic time (e.g. Kennicutt 1998; de los Reyes & Kennicutt 2019). Initially proposed by Schmidt (1959) in terms of volume densities within the Milky Way, this relation was later reformulated for extragalactic studies using measurements of surface densities (Kennicutt 1998). While extensively studied locally (e.g. Bigiel et al. 2008; Kennicutt & Evans 2012; Leroy et al. 2013; Pessa et al. 2021; Jiménez-Donaire et al. 2023), its high-redshift behaviour ($z \ge 2$) remains poorly constrained, with significant uncertainties in both the slope and intrinsic scatter of the relation. These uncertainties arise primarily from low spatial resolution, limited sensitivity (e.g. Freundlich et al. 2013), and sample representativeness issues (e.g. Vallini et al. 2024), hindering progress towards an understanding of early-Universe star

formation. Thanks to existing observational facilities such as Hubble Space Telescope (HST) and Atacama large millimetre/submillimetre array (ALMA), and the recent advancements made possible with the James Webb Space Telescope (JWST), we can now perform multi-wavelength analysis to probe galaxies at unprecedented redshifts and resolutions (e.g. Nagy et al. 2023; Hashimoto et al. 2023; Wang et al. 2023; Zavala et al. 2024). Those observations have opened new avenues for investigating the physical processes governing star formation in the first few billion years after the Big Bang.

Traditional molecular gas tracers, such as CO, become increasingly difficult to detect at higher redshift ($z > 2$, e.g. Molina et al. 2019; Tacconi et al. 2020), creating a critical gap in our understanding of early-Universe star formation (e.g. Genzel et al. 2010; Scoville et al. 2016). This limitation has been partially overcome by the use of the [CII] 158 µm emission line, now regularly used as an alternative tracer of molecular gas (Zanella et al. 2018). That fine-structure line is typically one of the most luminous far-infrared line in star-forming galaxies and is redshifted into ALMA’s atmospheric windows at $z \ge 4$. This provides unique insights into cold gas kinematics and is found to correlate with star formation (e.g. De Looze et al. 2014; Herrera-Camus et al. 2015). In particular, Herrera-Camus et al. (2015) discuss how [CII], as the main

<footnote>* Corresponding author: cedric.accard@astro.unistra.fr</footnote>
coolant of photo-dissociation regions (PDRs; see Wolfire et al. 2022), is physically linked to the SFR through the photoelectric effect, which couples the heating of the gas to the presence of young, massive stars. Several theoretical and observational studies have demonstrated that [CII] emission could also be strongly linked to molecular gas mass, particularly in the local Universe (e.g. Zanella et al. 2018; Madden et al. 2020; D'Eugenio et al. 2023). However, recent work at high redshift suggests a more complex scenario. For instance, Dessauges-Zavadsky et al. (2020) find that, in main-sequence galaxies at $z \sim 5$, the gas masses inferred from [CII] luminosity are comparable to those derived from dynamical mass and far-infrared continuum estimates (Scoville et al. 2014). This implies that either the total gas reservoir in these galaxies is predominantly molecular, or that [CII] traces the total gas mass, including both molecular and atomic components, rather than exclusively the molecular phase. The complex origin of [CII] emission, arising from both PDRs and diffuse ionised gas, further complicates its interpretation, necessitating careful analysis of observational data (e.g. Vallini et al. 2015; Wolfire et al. 2022).

The ALMA large programme to investigate C+ at early times (ALPINE; Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2022) has significantly expanded our sample of [CII]-detected galaxies at $4 < z < 6$, offering a first opportunity to study star formation relations in main-sequence galaxies during this epoch. This survey, comprising 118 star-forming galaxies, has already yielded important insights into the interstellar medium (ISM) properties (Vanderhoof et al. 2022; Veraldi et al. 2025) and star formation characteristics of galaxies in the early Universe (Schaerer et al. 2020; Faisst et al. 2022; Romano et al. 2022). Follow-up studies combining the higher resolution of four ALPINE galaxies with HST observations have allowed a more detailed view to be built of individual galaxy properties. The new ALMA observations, allowing one to resolve the spatial extent of [CII] emission on kiloparsec scales, provide information about their morpho-kinematics (Devereaux et al. 2024) and the distribution of star formation and cold gas linked through the KS relation (Béthermin et al. 2023, B23 in what follows). The latter shows that, despite their disturbed morphologies, likely due to mergers, the studied galaxies appeared to follow an extension of the local KS relation at higher gas surface densities. These studies have so far been limited to small sample sizes, which limits the statistical significance of their conclusions. To overcome this limitation, the [CII] resolved ISM in star-forming galaxies with ALMA (CRISTAL) large programme (2021.1.00280.L, Herrera-Camus et al. 2025) was initiated as an ALMA Cycle 8 follow-up to the ALPINE survey. The programme aims to obtain spatially resolved images of 19 [CII]-detected galaxies at $4.430 \le z \le 5.689$ with existing HST rest-frame UV emission images. These galaxies were selected based on their specific star formation rate (sSFR) within a factor of three of the star-forming main sequence at their respective redshifts, and with stellar masses of $M_{\star} \ge 10^{9.5} M_{\odot}$.

Leveraging the expanded capabilities of the ALPINE and CRISTAL surveys, we present an analysis of 13 main-sequence galaxies at $z \sim 5$, combining high-resolution data from JWST, HST, and ALMA. This multi-wavelength approach probes both dust-obscured and unobscured star formation, as well as the cold gas component traced by [CII]. The inclusion of JWST data significantly enhances our spectral coverage, enabling panchromatic rest frame UV-to-radio spectral energy distribution (SED) modelling for each line of sight, which was previously unattainable with HST and ALMA data alone. This was performed using the Python-based code investigating galaxy emission (CIGALE;

Boquien et al. 2019). With the production of spatially resolved maps below the kiloparsec scale of inferred physical properties, we can investigate the slope and scatter of $\Sigma_{[CII]} - \Sigma_{SFR}$ relation in the spatially resolved manner. We adopt a comprehensive statistical approach to handle uncertainties of the different measurements and estimates, which allows us to derive the slope and intrinsic scatter of the [CII]-SFR relation. We explore the impact of different [CII]-to-gas conversion factors on estimates of gas masses and star formation efficiencies via the KS relation.

The goal of this work is to determine whether the resolved [CII]–SFR and KS relations established at low redshift remain valid at $z \ge 4$, or if evolving ISM conditions and gas density regimes at high redshift significantly impact the link between cold gas and star formation. The structure of this paper is as follows. In Sect. 2, we describe the observations, sample selection, and data homogenisation procedures. Sect. 3 details our methodology, including the resolved SED fitting and the statistical treatment of uncertainties. Our main results, including the resolved [CII]-SFR relation and the impact of [CII]-to-gas conversion factors on the KS relation, are presented in Sect. 4. We discuss the implications of our findings in Sect. 5, and summarise our conclusions in Sect. 6. To ensure consistency with past works on the ALPINE survey, we assume a flat $\Lambda$ cold dark matter cosmology ($h = 0.7$, $\Omega_{\Lambda} = 0.7$, $\Omega_{m} = 0.3$) and a Chabrier initial mass function (IMF, Chabrier 2003).

## 2. Observations and data homogenisation

### 2.1. Sample selection

For this study, we selected a subset of 11 galaxies from the CRISTAL survey with homogeneous ALMA, HST, and JWST observations, for which there is the same filter coverage across all sources: HST/ACS F814W, JWST/NIRCam F115W, F150W, F277W, and F444W, along with ALMA [CII] line and associated dust continuum observations. As ALMA provides the coarsest spatial resolution among these datasets, the effective resolution of our analysis is set by the ALMA beam (see Sect. 2.2). This consistent multi-wavelength coverage enables the application of a uniform data homogenisation pipeline across the entire sample. In addition to the CRISTAL objects, we considered another galaxy from the ALPINE survey re-observed as part of the ALMA project 2019.1.00226.S (PI: E. Ibar). Although not part of the original CRISTAL survey, this galaxy (DEIMOS_COSMOS_873756) was assigned a CRISTAL identifier (CRISTAL_24). Finally, we included vuds_cosmos_5110377875 (hereafter VC875) from the ALPINE survey, as it meets our selection criteria (ALMA project 2022.1.01118.S, PI: M. Béthermin).

The final sample contains 13 star-forming main-sequence galaxies between $z = 4.439$ and $z = 5.689$, all located in the Cosmic Evolution Survey field (COSMOS; Scoville et al. 2007). Accurate redshifts for these galaxies have been computed from extensive optical spectroscopy (Le Fèvre et al. 2015; Hasinger et al. 2018) and then confirmed and refined with [CII] ALPINE data (Béthermin et al. 2020). A summary of the final sample is found in Table 1.

### 2.2. Observations

ALMA observations were conducted in Band 7, targeting the 158 µm rest-frame [CII] line ($v_{\text{rest}} = 1900.54$ GHz) using various antenna configurations. For the CRISTAL survey objects, the final data combined extended (C43-5 or C43-6; typical beam sizes 0.1"–0.3" and spatial scales of 1.2–2.0 kpc) and
<table><caption>**Table 1.** Summary of the ALPINE-CRISTAL sample used in this work.</caption><thead><tr><th>Original name</th><th>CRISTAL ID</th><th>z<sub>[CII]</sub></th><th>α<sub>J2000</sub><br>deg</th><th>δ<sub>J2000</sub><br>deg</th><th>log<sub>10</sub>(M<sub>★</sub>)<br>M<sub>☉</sub></th><th>log<sub>10</sub>(SFR)<br>M<sub>☉</sub> yr<sup>-1</sup></th><th>Beam size<br>arcsec × arcsec</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">CRISTAL large programme</td></tr><tr><td>DEIMOS_COSMOS_842313</td><td>CRISTAL-01</td><td>4.530</td><td>150.2271</td><td>+2.5762</td><td>10.65 ± 0.50</td><td>2.31 ± 0.74</td><td>0.18 × 0.18</td></tr><tr><td>DEIMOS_COSMOS_848185</td><td>CRISTAL-02</td><td>5.294</td><td>150.0896</td><td>+2.5864</td><td>10.30 ± 0.23</td><td>1.71 ± 0.29</td><td>0.33 × 0.28</td></tr><tr><td>DEIMOS_COSMOS_536534</td><td>CRISTAL-03</td><td>5.689</td><td>149.9719</td><td>+2.1182</td><td>10.40 ± 0.29</td><td>1.79 ± 0.31</td><td>0.51 × 0.40</td></tr><tr><td>vuds_cosmos_5100822662</td><td>CRISTAL-04a</td><td>4.520</td><td>149.7413</td><td>+2.0809</td><td>10.15 ± 0.29</td><td>1.89 ± 0.21</td><td>0.42 × 0.35</td></tr><tr><td>vuds_cosmos_5100541407</td><td>CRISTAL-06a</td><td>4.562</td><td>150.2538</td><td>+1.8094</td><td>10.09 ± 0.30</td><td>1.62 ± 0.34</td><td>0.37 × 0.28</td></tr><tr><td>DEIMOS_COSMOS_873321</td><td>CRISTAL-07ab</td><td>5.154</td><td>150.0166</td><td>+2.6266</td><td>10.00 ± 0.33</td><td>1.89 ± 0.26</td><td>0.45 × 0.28</td></tr><tr><td>DEIMOS_COSMOS_519281</td><td>CRISTAL-09a</td><td>5.575</td><td>149.7537</td><td>+2.0910</td><td>9.85 ± 0.39</td><td>1.51 ± 0.32</td><td>0.18 × 0.16</td></tr><tr><td>DEIMOS_COSMOS_630594</td><td>CRISTAL-11</td><td>4.439</td><td>150.1358</td><td>+2.2579</td><td>9.68 ± 0.33</td><td>1.57 ± 0.31</td><td>0.25 × 0.21</td></tr><tr><td>vuds_cosmos_5100994794</td><td>CRISTAL-13</td><td>4.579</td><td>150.1715</td><td>+2.2873</td><td>9.65 ± 0.34</td><td>1.51 ± 0.41</td><td>0.18 × 0.14</td></tr><tr><td>vuds_cosmos_5101244930</td><td>CRISTAL-15</td><td>4.580</td><td>150.1986</td><td>+2.3006</td><td>9.69 ± 0.33</td><td>1.44 ± 0.24</td><td>0.23 × 0.21</td></tr><tr><td>DEIMOS_COSMOS_494763</td><td>CRISTAL-19</td><td>5.233</td><td>150.0213</td><td>+2.0534</td><td>9.51 ± 0.36</td><td>1.45 ± 0.36</td><td>0.39 × 0.29</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">ALPINE follow-up</td></tr><tr><td>DEIMOS_COSMOS_873756</td><td>CRISTAL-24</td><td>4.546</td><td>150.0113</td><td>+2.6278</td><td>10.53 ± 0.08</td><td>2.06 ± 0.22</td><td>0.22 × 0.15</td></tr><tr><td>vuds_cosmos_5110377875</td><td>None<sup>a</sup></td><td>4.551</td><td>150.3847</td><td>+2.4084</td><td>10.16 ± 0.21</td><td>1.99 ± 0.26</td><td>0.40 × 0.35</td></tr></tbody></table>

**Notes.** The first column presents the original COSMOS survey names, while the second column shows CRISTAL IDs following the naming convention in Mitsuhashi et al. (2024). The $z_{[CII]}$ column lists redshifts derived from ALMA [C II] line observations of the CRISTAL survey (Herrera-Camus et al. 2025), for which formal uncertainty estimates are not provided. Nonetheless, these values are consistent with those from Béthermin et al. (2020), who report typical uncertainties of approximately 0.0005. Co-ordinates are given in J2000 epoch. The integrated mass and SFR are those reported in Herrera-Camus et al. (2025), apart from for VC875, whose values are from the Béthermin et al. (2020) catalogue. Beam sizes of the ALMA observations are reported as major axis × minor axis in arcseconds. (a) VC875 has not been used in previous CRISTAL work, therefore it has no CRISTAL ID.

compact array (C43-1 or C43-2; 0.7"–1.0", 4.5–6 kpc) configurations, incorporating available ancillary data. CRISTAL_24 was observed using C43-3 (~0.4", ~2.7 kpc, medium-resolution) and C43-5 (~0.2", ~1.5 kpc, high-resolution) configurations. The CRISTAL team consistently re-processed these data to produce [CII] moment and dust continuum emission maps using Briggs weighting (see Sect. 5 of Herrera-Camus et al. 2025). VC875 was observed in the C43-4 configuration (~0.3", ~2.0 kpc) and reduced independently (Béthermin et al. 2023).

The HST observations are cut-outs from the COSMOS survey (Koekemoer et al. 2007), obtained with the HST ACS camera and the F814W filter, at a FWHM resolution of 0.13". At the considered redshifts, this corresponds to rest-frame emission of 120–147 nm and spatial scales of 0.76–0.86 kpc.

All sources were observed with four JWST/NIRCam filters: F115W, F150W, F277W, and F444W, probing rest-frame emissions at 183–211 nm, 238–276 nm, 440–509 nm, and 705–816 nm, respectively. The achieved FWHM resolutions are 0.037", 0.049", 0.088", and 0.140", giving spatial scales of, respectively, 0.21–0.24 kpc, 0.29–0.33 kpc, 0.52–0.59 kpc, and 0.82–0.93 kpc, for our redshift range. These observations were part of the COSMOS-Web programme (PID: 1727; co-PIs: Casey & Kartaltepe; Casey et al. 2023).

An example of the complete multi-wavelength coverage available for our sample is presented in Fig. 1, for VC875. VC875 shows two UV-bright components in HST F814W and $\lambda < 3$ µm JWST bands, but the JWST/F444W and ALMA bands, as well as the stellar mass map from CIGALE (see Fig. 3), reveal a single component, supporting its treatment as a single galaxy with UV-bright clumps.

## 2.3. Homogenisation of the resolution

To perform accurate pixel-by-pixel comparisons across all available bands for our SED modelling process, we need to ensure that the information contained within each pixel corresponds to

the same region of the sky, for all observational bands. This consistency is crucial for deriving reliable physical properties from multi-wavelength data. Therefore, we performed point spread function (PSF) homogenisation across all bands to match their resolutions. We used the ALMA beam as the target resolution since it is the coarsest in our sample. The HST and JWST PSFs central lobes were modelled with 2D Gaussian functions to retrieve their full width at half maximum and maximum intensity.

Cleaned synthesised ALMA beams models were created based on dust continuum map header parameters, assuming a 2D elliptical Gaussian shape. Convolution kernels were generated to transform HST and JWST PSFs to match the target ALMA beam taking into account different grid orientations, pixel scales, and normalisation conventions. We re-projected HST and JWST convolved images onto the ALMA astrometric grid using the `reproject_interp` module from the Python package `reproject`, ensuring consistent pixel sampling and astrometric alignment across all bands. Applying the procedure to a point-like source with unit flux in HST/JWST native resolution, we find the reconstructed average flux to be 2.6% lower than the expected peak value of unity across the different filters. The missing flux likely resides in the secondary lobes not taken into account, and we consider this effect to be negligible compared to our overall flux calibration uncertainties.

Figure 2 shows the effect of the resolution homogenisation of the JWST/NIRCam F150W filter in the case of VC875. The complete set of homogenised observations for this source is found in Fig. A.1. As a result of this resolution homogenisation, we achieve an average beam size of 0.3" (~2 kpc) across the sample. Although the pixel scale is ~0.2 kpc, the physical information is correlated on the scale of the beam. We take into account these correlations induced by our resolution homogenisation in our analysis of the [CII]-SFR relation (see Sect. 3.3).

<footnote n="1">¹ https://reproject.readthedocs.io</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Multi-wavelength data for the galaxy VC875, shown at their native resolutions. From left to right: The first row presents HST/ACS F814W (rest-frame 145 nm), JWST/NIRCam F115W (207 nm), F150W (270 nm), and F277W (499 nm). The second row shows JWST/F444W (800 nm), ALMA dust continuum emission around 157 µm, and ALMA [CII] 158 µm emission. The contours are those of the ALMA [CII] line emission maps starting from 3σ increasing by steps of 2σ. Ellipses in the bottom left corners represent the central lobe of the PSF for HST and JWST images, and the synthesised beam for ALMA images.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Native resolution image of VC875 (left panel) and the homogenised version matching the ALMA resolution (right panel) for the JWST/NIRCam F150W filter. The contours and ellipses are the same as in Fig. 1.</figcaption></figure>

## 2.4. Photometric and line flux unit conversion

To prepare the data for SED modelling, we converted all photometric data to micro jansky per beam, and ALMA dust continuum emission flux to dust luminosity. Here ‘beam’ refers to the synthesised ALMA beam, not the PSF of each instrument. Expressing fluxes per beam ensures that the measurements are referenced to the area over which the instrument is sensitive. The dust luminosity, $L_{IR}$ (W beam$^{-1}$), was calculated from the monochromatic flux, $S_v$ (Jy beam$^{-1}$), using

$$L_{IR} = 10^{-26} \times \frac{1}{0.13} \nu_{em} \frac{4\pi D_L^2}{1+z} S_v, \quad (1)$$

where $\nu_{em}$ is the rest-frame emission frequency at 158 microns (hertz), $D_L$ is the luminosity distance (metres), and $z$ is the redshift. The factor $10^{-26}$ converts jansky to watts per square metre per hertz, while the $1/(1+z)$ term accounts for cosmological redshift. We adopted the empirical conversion factor of 0.13 ± 0.02 from Khusanova et al. (2021) to convert the 158 µm rest-frame monochromatic luminosity into the total infrared luminosity, $L_{IR}$. This value corresponds to the calibration for the redshift

bin $4 < z < 5$ and is directly applicable to our sample, which predominantly lies within this range. The other sources are in the $5 < z < 6$ bin, in which Khusanova et al. (2021) measured a conversion factor of 0.12 ± 0.03. Since the two values are compatible with no redshift dependence, we use the more accurate value found at $4 < z < 5$ for the full sample. Both factors were calibrated using stacked SEDs. A compilation of SED templates from the literature was fitted to these SEDs and only the good fits were kept (reduced $\chi^2 < 1.5$). The conversion factor was then estimated using the mean and dispersion of these templates. The typical dust temperatures corresponding to these stacked SEDs are $T_d = 41 \pm 1$ K at $z \sim 4.5$ and $T_d = 43 \pm 5$ K at $z \sim 5.5$ for galaxies with SFR > 10 $M_{\odot}$ yr$^{-1}$ (see Appendix A of Khusanova et al. 2021). We note that adopting a constant conversion factor assumes a uniform dust temperature across the sample, and variations could introduce systematic effects on $L_{IR}$ and thus $\Sigma_{SFR}$ estimates (Faisst et al. 2017; Cochrane et al. 2022). We include this IR dust luminosity in the SED modelling process as a constraint on dust attenuation models (see Sect. 3.1).

The [CII] surface brightness ($\Sigma_{[CII]}$ in solar luminosities per square kiloparsec) was derived from the [CII] moment-0 map ($m_{[CII]}$ in jansky kilometre per second per beam) using

$$\Sigma_{[CII]} = \frac{m_{[CII]}}{D_A^2 \Omega_{\text{Beam}}} \left( 1.04 \times 10^{-3} \frac{L_{\odot}}{\text{GHz Mpc}^2 \text{Jy km s}^{-1}} \right) D_L^2 \nu_{\text{obs}}, \quad (2)$$

where $D_A^2 \Omega_{\text{Beam}}$ is the physical area associated with the synthesised beam, defined by $D_A$, the angular distance (in kiloparsecs), and $\Omega_{\text{Beam}}$, the solid angle of the beam. The remaining factors correspond to the conversion from line flux to luminosity (Solomon et al. 1992), where $D_L$ is the luminosity distance (in megaparsecs) and $\nu_{\text{obs}}$ the observed frequency (in gigahertz).

## 2.5. Uncertainty estimation

For HST and JWST images, we incorporated the convolution systematic error from Section 2.3. In addition, we applied a
relative uncertainty of 3% to JWST images to account for potential calibration uncertainties (Gordon et al. 2018). The conversion from continuum flux to $L_{IR}$ relies on an empirical factor of 0.13 ± 0.02, which introduces a relative systematic uncertainty of approximately 15% (0.02/0.13) to the derived $L_{IR}$ values in our continuum images. For images lacking error maps (HST data and $L_{IR}$), we assumed white noise at the native resolution and estimated the measurement uncertainty using the standard deviation of the noise in the convolved then reprojected image. For $\Sigma_{[CII]}$, we used the pixel-by-pixel error maps provided by the CRISTAL team.

JWST maps are provided with their associated error maps. We produced variance maps of the resolution-homogenised maps (in jansky per pixel) by convolving the noise variance map (error map squared) by the kernel:

$$ \sigma_{\text{conv}}^2 = \sigma_{\text{pixel}}^2 * K^2, \qquad (3) $$

where $\sigma_{\text{conv}}^2$ and $\sigma_{\text{pixel}}^2$ are the noise variance of the convolved and original maps, respectively, and K is the homogenisation kernel. We compared the results with the standard deviation of the noise in an empty area of the homogenised map, and we found a small excess, likely due to a small contribution from correlated noise already present in the original maps and not taken into account by our computation. This excess remains negligible for the short-wavelength filters (less than 6% for F115W and F150W), but becomes increasingly significant at longer wavelengths, reaching up to 14% for F277W and F444W. We thus applied a scaling factor (1–1.14, depending on the source and band) to the resolution-homogenised error maps to correct for this effect.

To account for this error budget in our analysis, we generated 100 perturbed realisations of each image. This process involves adding white noise with a standard deviation equal to the measured pixel-level noise at the native resolution of each image, followed by the application of systematic errors. The HST and JWST perturbed images were subsequently convolved. This method accounts for the non-independence of the pixels from spatially correlated noise and the smoothing of the maps.

# 3. Methods

## 3.1. Resolved SED fitting

Our primary objective is to derive the SFR along each line of sight for the galaxies in our sample by combining multi-wavelength observations. At a spatial resolution of 2 kpc, the use of an energy-conserving SED modelling code is justified, as the energy balance assumption-where UV/optical energy absorbed by dust is re-emitted in the IR-remains physically meaningful at this scale (e.g. Boquien et al. 2015).

To this end, we used CIGALE, a versatile Python tool capable of modelling galaxy emission from X-ray to radio wavelengths (Boquien et al. 2019). Its modular design allows for the inclusion of diverse SFHs, dust attenuation laws, and emission line models tailored to specific scientific objectives. CIGALE uses a Bayesian approach to perform parameter estimation, generating a large library of theoretical SEDs based on user-defined models. Observed data were then compared to these models, and posterior probability distributions were calculated for each parameter. This Bayesian framework ensures that uncertainties and degeneracies between parameters are properly propagated into the resulting parameter estimates.

In this section, we present the SED modelling parameters used for our analysis. We followed Boquien et al. (2022), who

analysed integrated measurements of the ALPINE galaxies (our parent sample), and Buat et al. (2019), who focused on dust-rich galaxies at $z \sim 2$. We refer the reader to those works for a comprehensive discussion of parameter choices. Our focus here is to detail the differences. A summary of our parameters is given in Table B.1.

We adopted delayed SFH model with an additional recent variation, as is described by Ciesla et al. (2017). In this framework, star formation begins after an initial delay corresponding to the time between the Big Bang and the onset of the first stars (here set to 250 Myr). The main stellar population forms according to a delayed SFH, which rises linearly from the onset, peaks at a characteristic timescale, and then declines exponentially. To capture recent deviations from this smooth evolution, we introduced a modification to the SFR at a specified time in the recent past. The SFH is thus expressed as

$$ SFR(t) = \begin{cases} t \exp(-t/\tau), & t < t_{\text{var}} \\ r_{\text{SFR}} t \exp(-t/\tau), & t \ge t_{\text{var}} \end{cases} \qquad (4) $$

where t is the time since the onset of star formation, $\tau$ is the SFH timescale, $t_{\text{var}}$ is the time since the onset of the recent SFH variation, and $r_{\text{SFR}}$ is the factor by which the SFR is increased or decreased. This approach provides flexibility to model both the overall stellar mass assembly and any recent changes in star formation activity that a simple delayed SFH model would not account for.

For all sources, we adopted the two-component dust attenuation model of Charlot & Fall (2000). This choice is motivated by its physical basis and agreement with radiative transfer predictions, allowing it to be applicable across the full range of dust attenuation, from weakly to highly obscured regions. While the modified Calzetti et al. (1994, 2000) relation (Noll et al. 2009) is often used for star-forming galaxies, it is less appropriate for low attenuation or highly obscured systems (e.g. Boquien et al. 2022), and allowing its slope parameter ($\delta$) to vary could lead to overfitting given our limited number of photometric bands. Although for most of our sample there are no significant differences between the two models (as is illustrated by VC875; see Appendix C), the Charlot & Fall (2000) law provides a substantially better agreement between SFR estimators for CRISTAL_24 (see Fig. C.2), further consolidating our choice of attenuation law.

We also evaluated the impact of including the ALMA-derived total infrared luminosity ($L_{IR}$, computed from the 158 µm continuum in Sect. 2.4) as a dust luminosity constraint in our SED modelling. We find that omitting this constraint causes the models to systematically overestimate the SFR compared to estimates from combined UV and far-IR data. This finding is consistent with Li et al. (2024), who reported a similar discrepancy in an overlapping sample. Their analysis, using region-averaged measurements with the magphys (da Cunha & Charlot 2011) SED fitting code, likewise highlights the necessity of the dust luminosity constraint in SED modelling. Our results, derived from a pixel-based CIGALE approach, reinforce this conclusion.

In Appendix C, we discuss our final modelling choices in detail, using the galaxy VC875 as a case study. These choices include the dust luminosity constraint, the adopted dust attenuation law, and the treatment of recent star formation episodes. A similar analysis was performed for all other sources in the sample.

To ensure robust measurements, we used only pixels with $\Sigma_{[CII]} \ge 5\sigma_{[CII]}$ and S/N$\ge 2$ in at least three bands. The latter
<figure><img src="image_4.png" /><figcaption>Fig. 3. Example of physical parameters maps obtained with the CIGALE SED modelling for VC875. From left to right: Stellar mass surface density (solar masses per square kiloparsec), SFR surface density averaged over 10 million years (solar masses per year per square kiloparsec), mass-weighted age (million years), and attenuation in the FUV rest frame (magnitudes). The contours are those of the [CII] luminosity map, starting at 3σ and increasing by steps of 2σ. The ellipse in the bottom left of each panel represent the synthesised ALMA beam. The maps are shown for the pixels having $L_{[CII]} \ge 5\sigma_{[CII]}$ and S/N ≥ 2 in at least three bands.</figcaption></figure>

cut, applied to the data used for SFR estimates, excludes only about 6% of pixels that passed the $L_{[CII]}$ threshold, confirming that our selection is primarily driven by [CII] detectability. This ensures the robustness of the [CII] measurements, while minimising biases in the inferred SFR distribution.

The SFR uncertainties were estimated from the 16th, 50th, and 84th percentiles of the log(SFR) distributions obtained via resampling noise realisations (see Sect. 2.5). This approach accounts for pixel correlations introduced by resolution-matching. It avoids biases from assuming Gaussian errors, as the SFR distributions are typically skewed.

Figure 3 presents the maps of stellar mass, SFR, mass-weighted age, and FUV attenuation for VC875, derived from SED modelling after applying all pixel selection criteria. These spatially resolved maps help to verify the reliability of our results by revealing whether parameters vary smoothly and align with the galaxy’s morphology. These maps also make it easier to search for modelling artefacts, such as abrupt discontinuities or parameter saturation. This enables a direct comparison of physical property distributions with features seen at other wavelengths or in simulations, providing a useful check on our assumptions. Across the sample, we find that stellar mass and SFR peaks are generally aligned. In contrast, the mass-weighted ages and FUV attenuation exhibit more complex and asymmetric patterns (see Appendix F).

## 3.2. Fitting the SFR-[CII] relation

To study the relationship between the SFR obtained by the SED modelling and the [CII] luminosity, we parametrised this relation as a power law with an intrinsic scatter that accounts for the variations in the physical conditions:

$$ \log_{10}(\Sigma_{\text{SFR}}) = \log_{10}(C) + \beta \log_{10}(\Sigma_{[\text{CII}]}) \pm \sigma_i, \quad (5) $$

where $\beta$, $C$, $\Sigma_{[CII]}$, $\Sigma_{SFR}$, and $\sigma_i$ are the slope, the normalisation constant, the [CII] surface brightness, the SFR surface density, and the intrinsic scatter, respectively.

To estimate the slope and scatter of the [CII]-SFR relation, we used a Bayesian fit relying on a maximum likelihood estimator. We constructed a likelihood that properly accounts for both normal distribution of [CII] maps error and log-normal nature of SFR distributions, but also estimates the intrinsic scatter. To mitigate potential degeneracies between the normalisation factor and the slope when fitting the data, we introduced a normalised [CII] surface brightness:

$$ g_{\text{k,norm}} = g_k/g_0, \quad (6) $$

where $g_k$ represents the [CII] surface brightness measurement for each data point, $k$, and the subscript norm indicate normalisation by the reference value, $g_0$. The parameter $g_0$ is a fixed normalisation pivot introduced to reduce the degeneracies between the slope, $\beta$, and the normalisation constant, $C$, and was set to the median [CII] surface brightness in our sample, ($\log_{10}(g_0) = 7.5$ L$_{\odot}$ kpc$^{-2}$). The slope, $\beta$, is independent of the choice of $g_0$, while the normalisation constant, $C$, is expressed relative to this pivot (denoted $C_0$). The sample $\log_{10} g_k$ values ranges from 6.89 to 8.51, with an interquartile range of 7.23 to 7.92, confirming the representativeness of this choice.

Next, we addressed the statistical treatment of our data. For each data point, $k$, in the [CII]-SFR plane, we considered the normalised [CII] surface brightness ($g_{k,n}$) and SFR surface density ($s_k$). To account for the uncertainties in the [CII] measurements, we modelled the true value $g_t$ of a given $g_{k,n}$ measurement as a normal distribution with standard deviation $\sigma_{g,k}$. This uncertainty is incorporated into the likelihood through marginalisation over the $g_t$ parameter.

Since the SFR estimates followed a log-normal distribution, we worked with the log-quantity, denoted as $y_k = \ln(s_k)$, which then followed a normal distribution. The mean of this distribution is given by $\mu = \ln(C_0) + \beta \ln(g_t)$. Finally, the standard deviation of this new distribution, $\sigma_T$, accounts for both the SFR estimate error and the intrinsic scatter of the relation

$$ \sigma_{T,k} = \sqrt{\sigma_i^2 + \sigma_{y,k}^2}, \quad (7) $$

where $\sigma_{y,k}$ represents the SFR estimate error for a pixel, $k$, computed from the $y_k = \ln(s_k)$ distribution and $\sigma_i$, the intrinsic scatter of the SFR-[CII] relation.

With this description, the likelihood for a given pixel is expressed as

$$ p(y_k, \sigma_{y,k}, g_{k,n}, \sigma_{g,k} | C_0, \beta, \sigma_i) = \frac{\int_{-\infty}^{+\infty} p(y_k, \sigma_{y,k}, g_t | C_0, \beta, \sigma_i) p(g_t | g_{k,n}, \sigma_{g,k}) dg_t}{2\pi \sigma_T \sigma_{g,k} \int_{-\infty}^{+\infty} \exp\left(-\frac{(y_k - \mu)^2}{2\sigma_T^2}\right) \exp\left(-\frac{(g_{k,n} - g_t)^2}{2\sigma_{g,k}^2}\right) dg_t} \quad (8) $$

The first term in the likelihood, $p(y_k, \sigma_{y,k} | g_t, C_0, \beta, \sigma_i)$, represents the KS relation and its associated uncertainties, while the second term, $p(g_t | g_{k,n}, \sigma_{g,k})$, accounts for the observational uncertainty on the [CII] surface brightness measurements.
<table><caption>**Table 2.** Best-fit parameters and uncertainties for the $\Sigma_{SFR}-\Sigma_{CII}$ relation derived from the different fitting methods described in Sect. 3.2.</caption><thead><tr><th>Fitting Method</th><th>log(C<sub>0</sub>)<sup>a</sup></th><th>β</th><th>σ<sub>i</sub></th></tr></thead><tbody><tr><td colspan="4">Likelihood Model + L-BFGS-B integration</td></tr><tr><td>Original Fit</td><td>0.30</td><td>0.87</td><td>0.19</td></tr><tr><td>Source Bootstrapping</td><td>0.30<sup>+0.05</sup><sub>-0.04</sub></td><td>0.86<sup>+0.15</sup><sub>-0.11</sub></td><td>0.19<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>Noise Sampling</td><td>0.33<sup>+0.02</sup><sub>-0.02</sub></td><td>0.80<sup>+0.05</sup><sub>-0.06</sub></td><td>0.20<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Hybrid Resampling</td><td>0.34<sup>+0.05</sup><sub>-0.04</sub></td><td>0.80<sup>+0.15</sup><sub>-0.13</sub></td><td>0.19<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td colspan="4">Likelihood Model + MCMC</td></tr><tr><td>Original Likelihood</td><td>0.30</td><td>0.87</td><td>0.19</td></tr><tr><td>Normalised Likelihood</td><td>0.30<sup>+0.03</sup><sub>-0.03</sub></td><td>0.87<sup>+0.09</sup><sub>-0.09</sub></td><td>0.20<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td colspan="4">Linmix</td></tr><tr><td>All Pixels</td><td>0.37</td><td>0.92</td><td>0.18</td></tr><tr><td>Region Averages (B23)</td><td>0.22<sup>+0.04</sup><sub>-0.04</sub></td><td>1.02<sup>+0.12</sup><sub>-0.12</sub></td><td>0.06<sup>+0.08</sup><sub>-0.08</sub></td></tr></tbody></table>

**Notes.** Parameter uncertainties, when given, were computed using the 16th, 50th, and 84th percentiles of the associated distributions. <sup>(a)</sup>We present the log(C<sub>0</sub>) using a normalisation factor ($g_0 = 10^{7.5}$), since the reported uncertainties have been derived using this $g_0$ cannot easily be extrapolated back to the un-normalised space. However, a subtraction of 7.5 from log(C<sub>0</sub>) allows one to obtain the actual relation intercept.

By marginalising over $g_i$, we propagated these uncertainties into our parameter estimation. The marginalisation was performed numerically by integrating on the $g_{k,n} \pm 5(\sigma_{g,k})$ range spanned by a grid of 1000 linearly spaced values. This covers 99.99994% of the probability mass for a normal distribution, and thus captures nearly the entire range of plausible true [CII] surface brightness values. To estimate the model parameters, we maximised the total log-likelihood across all data points. This is equivalent to minimising the following negative log-likelihood:

$$ \sum_k -\ln(p(y_k|g_{k,n}, \sigma_{g,k}, C_0, \beta, \sigma_T)). \quad (9) $$

Physical constraints enforcing positivity ($C > 0$, $\sigma_i > 0$) were imposed via flat, uniform priors restricted to the physically meaningful domain of positive values. Parameter optimisation was performed using L-BFGS-B (from `scipy.minimize` function), a boundary-aware algorithm specifically designed for constrained optimisation.

To assess potential biases in this fitting methodology, we performed an end-to-end numerical simulation that replicates the full observational chain, from mock source generation, noise injection, and resolution homogenisation to the inference of the physical relation parameters. This revealed no systematic parameter estimation biases on a simulated population similar to our sample (see Appendix D). However, when considering individual galaxies with a limited amount of independent data points, we obtained unreliable results, highlighting the necessity of combining multiple objects to recover the proper physical relation.

## 3.3. Uncertainty quantification

To evaluate the robustness of our results and quantify uncertainties in the model parameters, we implemented a systematic resampling strategy. This approach used the ensemble of the 100 perturbed [CII] maps per source in our sample and their corresponding SFR estimates (derived in Sects. 2.5 and 3.1), combined with the unperturbed dataset. We investigated three distinct methods to disentangle sources of uncertainty, applying each method 1000 times to generate robust parameter distributions. The first one, called source bootstrapping, consisted

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

**Fig. 4.** Resolved [CII]-SFR relation for galaxies at $z \sim 5$. The grey contours show the individual pixel distribution for the whole sample. Each pair of coloured points represents a high- and low-[CII] density average region for each galaxy of the sample following B23 methodology. The orange line and shaded region shows our best-fit power-law model from the tailored likelihood-based approach with the parameters in the top left corner (where $log_{10}(C_0)$ is the value at $g_0 = 10^{7.5} L_\odot \text{ kpc}^{-2}$). The dashed black line presents the extrapolation of the relation derived from a spatially resolved sample of 46 nearby galaxies by Herrera-Camus et al. (2015). We represented the typical errors bars for the individual pixels across the sample with three black squares.

of selecting 13 sources with replacement of the original noise realisation. This was used to test the sample representativeness. The second one, called noise sampling, selected the original 13 sources but with a random noise realisation, enabling an assessment of the impact of observational noise. Finally, we performed a hybrid resampling by selecting 13 sources randomly as well as a random noise realisation to combine both effects. For the last two methods, we accounted for the additional noise introduced by the resampling process by using $\sqrt{2}$ times larger uncertainties in our fitting procedure.

In addition to our primary resampling methods, we explored alternative fitting methods. Details of these supplementary analyses, including `linmix`<fnref n="2" /> Bayesian regression and Markov chain Monte Carlo (MCMC) techniques, are presented in Appendix E. The results of these secondary approaches are summarised in Table 2 for direct comparison with the primary resampling methods, and are discussed in Sect. 4.1.

# 4. Results

## 4.1. Resolved SFR-[CII] relation

We began by fitting the resolved [CII]-SFR surface densities relation for our $z \sim 5$ galaxy sample using our tailored likelihood-based fitting routine, applied to the original noise realisation of all 13 galaxies. This yielded a best-fit power-law model with $log(C_0) = 0.30, \beta = 0.87$, and $\sigma_i = 0.19 \text{ dex}$. At this

<footnote n="2">² Python port of Kelly (2007) LINMIX_ERR IDL code, available at https://github.com/jmeyers314/linmix.git.</footnote>
stage, we just derived the best fit, and formal uncertainties were not directly obtained. To assess the precision of these results, we compared them to those obtained with our primary resampling strategies. Using population bootstrapping, we found parameter estimates that are consistent with the baseline fit and that show uniform relative uncertainties of ~15% for all parameters: $log(C_0) = 0.30^{+0.05}_{-0.04}$, $\beta = 0.86^{+0.15}_{-0.12}$, and $\sigma_i = 0.19^{+0.03}_{-0.03}$ dex. These uncertainties remain stable across 1000 resamplings, despite the modest sample size ($N = 13$), indicating that the measured relation is not dominated by individual outliers. In contrast, the noise sampling approach yields a slightly flatter slope ($\beta = 0.80^{+0.05}_{-0.06}$) and a marginally higher intrinsic scatter ($\sigma_i = 0.20^{+0.02}_{-0.02}$). The hybrid resampling method shows a combination of both trends ($\beta = 0.80^{+0.15}_{-0.13}$, $\sigma_i = 0.19^{+0.04}_{-0.03}$). The largest error bars on the slope arise from source bootstrapping, highlighting the influence of sample selection, whereas noise resampling yields lower slope values. Importantly, the intrinsic scatter remains remarkably stable across all methods, providing a reliable measure of how tightly [CII] emission traces star formation in our sample. These results are summarised in Table 2.

Figure 4 shows the relation between SFR surface density and [CII] surface brightness for our $z \sim 5$ sample. We present the best-fit model derived from our tailored likelihood-based approach for both the original sample and the region-averaged pairs following the B23 methodology. For comparison, we also plot the relation from Herrera-Camus et al. (2015), which is based on spatially resolved measurements at ~0.2–1.5 kpc scales for 46 nearby galaxies ($d < 30$ Mpc). While the majority of our $z \sim 5$ sample lies above their local [CII]–SFR relation, our resolved analysis yields a shallower slope and similar intrinsic scatter compared to the lower-redshift study of Herrera-Camus et al. ($\beta = 1.13 \pm 0.01$, $\sigma_i = 0.21$ dex, reported without uncertainty). In addition to the redshift difference, the parameter space covered by the two samples is very different: the Herrera-Camus et al. (2015) sample spans $4.4 < \log \Sigma_{[CII]} (L_{\odot} \text{ kpc}^{-2}) < 6.5$, whereas our sample covers $6.8 < \log \Sigma_{[CII]} < 8.6$. This difference likely contributes to the observed slope discrepancy as it may reflect the evolution in ISM physical conditions across cosmic time, such as increased turbulence, lower metallicity, or changes in ISM structure at high redshift, yielding different dynamical range in $\Sigma_{[CII]}$ as well as $\Sigma_{SFR}$. We emphasise that methodological differences also play a critical role, especially on observation-based studies. For instance, Herrera-Camus et al. (2015) employed a least-squares bisector, which neglects pixel covariance. In contrast, our tailored approach explicitly accounts for pixel covariance and highlights the sensitivity of the slope estimate to this effect, as is demonstrated by the various tests conducted below.

We also explored alternative fitting methods to further validate our findings, which are presented in Appendix E. A MCMC analysis of the unperturbed dataset reproduces the best-fit parameters of the baseline fit, but initially yields unrealistically narrow error bars. For example, the MCMC approach yields an intrinsic scatter of $\sigma_i = 0.19 \pm 0.002$, showing error bars an order of magnitude smaller than uncertainties obtained via hybrid resampling (~0.03). This discrepancy likely arises from neglected correlations in the data. To correct for this, we used the likelihood normalised by half the number of pixels per beam, resulting in more realistic uncertainties that closely match those from the hybrid resampling approach (e.g. $\sigma_i = 0.20 \pm 0.02$). Using the linmix Bayesian regression method on the original dataset produces a slightly steeper slope ($\beta = 0.92 \pm 0.02$) and comparable intrinsic scatter ($\sigma_i = 0.18 \pm 0.04$ dex), again with very narrow uncertainties. To mitigate the impact of pixel

covariance, we also applied linmix to region-averaged points following the B23 methodology that consists of averaging the SFR and $\Sigma_{[CII]}$ measurements over large regions to obtain a sufficiently elevated signal-to-noise ratio and by construction gives a pair of points per galaxy corresponding to the low and high regime of [CII] surface brightness (pairs of points in Fig. 4 and Fig. 5). This approach yields a steeper slope than the pixel-based fit ($\beta = 1.02 \pm 0.12$) but a smaller and poorly constrained scatter ($\sigma_i = 0.06 \pm 0.08$ dex).

The MCMC and resampling methods yield parameter estimates that agree within ~1$\sigma$, demonstrating the robustness of our resolved analysis. In contrast, the linmix method yields a slightly higher slope and a poorly constrained intrinsic scatter, indicating a mild tension with the other primary approaches. These comparisons highlight the importance of properly accounting for pixel covariance, as methods that neglect this effect can underestimate uncertainties or fail to reliably constrain the intrinsic scatter.

Because of the limitations of our model when applied on a per-galaxy basis – that is, with restricted $\Sigma_{[CII]}$ dynamic ranges (see Sect. 3.2) – we were unable to recover resolved [CII]–SFR relations per galaxy. We therefore adopted the linmix method for simplicity, acknowledging that it does not properly account for pixel correlations. This method yields an average slope of 1.27, which is consistent with the mean value of 1.21 derived from Li et al. (2024) using spatially bin-averaged measurements on the same ALPINE-CRISTAL sample. As this value is manually measured from their Fig. 12, no formal uncertainty is provided. This agreement suggests that, despite the methodological differences and the limited dynamic range within single galaxies, the underlying scaling relation inferred from spatially resolved data remains broadly compatible with results based on coarser spatial binning.

Our results offer a valuable comparison to the theoretical predictions of Ferrara et al. (2019), who used spatially resolved models to investigate the [CII]-SFR relation within galaxies. Their work explored how local ISM conditions (e.g. density, metallicity) shape the resolved $\Sigma_{[CII]} - \Sigma_{SFR}$ relation, finding that it steepens at high surface densities due to $\Sigma_{[CII]}$ saturation. We observe a similar trend on local scales, with the mean slope within individual galaxies aligning with their predictions. However, our relation is significantly shallower for the population as a whole, confirming that this steepening is a local effect. This contrast highlights the critical role of spatial scale in interpreting the [CII]-SFR relation. Furthermore, it underscores the need for both local and population-wide perspectives to fully understand the diverse ISM conditions at high redshift.

## 4.2. [CII]-to-gas conversion factors

To estimate the gas mass, we first used the prescription of Zanella et al. (2018): the conversion factor is constant with a value of $\alpha_{[CII]} = 31 M_{\odot} L_{\odot}^{-1}$, and the gas mass surface density was derived as

$$ \Sigma_{\text{gas}} = \alpha_{[\text{CII}]} \Sigma_{[\text{CII}]} \qquad (10) $$

It is well established that [CII] emission can originate from multiple ISM phases, primarily PDRs and molecular gas, yet the contributions of these phases to the total [CII] luminosity remain debated and are sensitive to local ISM conditions. Consequently, the dominant gas phase traced by [CII] can vary across different environments, particularly at high redshift where ISM properties often diverge from those in local galaxies. This diversity
<figure><img src="image_6.png" /><figcaption>Fig. 5. Kennicutt-Schmidt relation using the $\alpha_{[CII]}$ conversion factor (left) and the $W_{[CII]}$ version (right). The grey contours show the distribution of individual pixel measurements across the entire sample; the pair of points are the B23 region averages for individual galaxies, colour-coded by redshift. The dashed grey lines show constant gas depletion time. The dashed black line represents the Wang et al. (2022) KS relation for local Universe and high-redshift main-sequence galaxies. We show CO measurements of the global KS relation from the low-z COLD GASS sample (Saintonge et al. 2012, crosses) and at redshifts up to z ~ 2.5 from the PHIBBS (Tacconi et al. 2013, plus signs) and PHIBBS2 (Freundlich et al. 2019, three-branch stars) programmes, together with a resolved z = 2.6 lensed starburst (Sharon et al. 2013, red points).</figcaption></figure>

has a direct impact on the appropriate [CII]-to-gas conversion factor (e.g. Madden et al. 2020; Vizgan et al. 2022). Recent theoretical work reinforces this perspective: Vallini et al. (2025), using galaxies from the SERRA cosmological zoom-in simulation (Pallottini et al. 2022), derive a relation varying with the [CII] surface brightness:

$$ \log W_{[CII]} = -0.506 \log \Sigma_{[CII]} + 4.933, \quad (11) $$

where the $W_{[CII]}$ notation is used to express the resolved, surface-brightness-dependent conversion factor, in contrast to the global $\alpha_{[CII]}$ value. This relation exhibits a typical $1\sigma$ dispersion of 0.4 dex. The dependence on $\Sigma_{[CII]}$ encapsulates the expectation that brighter [CII] regions are denser and more metal-enriched, two conditions that drive higher [CII] emission efficiency per unit gas mass. Across the $6.8 < \log \Sigma_{[CII]} (L_{\odot} \text{ kpc}^{-2}) < 8.6$ range probed in our sample, $W_{[CII]}$ varies from $31 M_{\odot} L_{\odot}^{-1}$ in the faintest regions down to $3.8 M_{\odot} L_{\odot}^{-1}$ in the brightest regions. This prescription can therefore lead to a steeper slope in the KS relation derived from [CII], particularly at high surface densities. It is important to note that the conversion from Vallini et al. (2025) yields a relation computed for the total cold gas content. This total content encompasses both atomic and molecular phases, thereby accounting for the outer layers of PDRs and the fully molecular clumps within clouds.

To explore the implications of this conversion factor on our interpretation of [CII]-based observations, we included it in our analysis. Using the $W_{[CII]}$ formula from Eq. (11) in place of $\alpha_{[CII]}$ in Eq. (10) yields

$$ \Sigma_{[CII]} \propto \Sigma_{\text{gas}}^{1/0.494}. \quad (12) $$

Substituting this relation into Eq. 5 gives

$$ \Sigma_{\text{SFR}} \propto \Sigma_{[CII]}^{\beta} \quad \text{and thus} \quad \Sigma_{\text{SFR}} \propto \Sigma_{\text{gas}}^{\beta/0.494}. \quad (13) $$

This implies that for a given slope, $\beta$, of the SFR-[CII] relation, the slope of the SFR-gas relation is expected to be approximately twice as steep. In contrast, the constant $\alpha_{[CII]}$ prescription leaves the slope, $\beta$, unaffected.

### 4.3. Impact on the Kennicutt-Schmidt relation

Figure 5 presents the resolved KS relation for our $z \sim 5$ sample, contrasting two [CII]-to-gas conversion prescriptions. The left panel adopts the constant $\alpha_{[CII]}$ from Zanella et al. (2018), and the right applies the surface-brightness-dependent $W_{[CII]}$ from Vallini et al. (2025). In both cases, we overlay a set of reference data from CO-based global KS relations from the low-redshift COLD GASS sample (Saintonge et al. 2012), as well as higher-redshift measurements from PHIBBS (Tacconi et al. 2013) and PHIBBS2 (Freundlich et al. 2019). The resolved $z = 2.6$ lensed-starburst of Sharon et al. (2013) is also shown, providing a benchmark for extreme star-forming environments. The dashed black line in both panels represents the KS relation for local and high-redshift main-sequence galaxies from Wang et al. (2022). Due to the substantial dispersion inherent to both conversion factor prescriptions, our fitting methodology cannot robustly constrain the intrinsic scatter of the relation. As such, the following discussion provides a qualitative assessment of the observed trends, and reported parameter values should be interpreted with caution.

When adopting a constant $\alpha_{[CII]}$, our sample is compatible with the Wang et al. (2022) results, despite a small slope difference ($\beta = 1.13 \pm 0.09$ for Wang et al. 2022, $\beta = 0.87 \pm 0.06$ in this work). The data align with constant gas depletion times of ~0.5–1 Gyr, extending the locus of the low- and intermediate-redshift CO-based studies (Saintonge et al. 2012; Tacconi et al. 2013; Freundlich et al. 2019) to higher gas and SFR surface densities. Interestingly, these are ~2–3 times shorter than those usually found for local star-forming spirals at the similar physical scales using CO-based $H_2$ tracers (e.g. Schruba et al. 2011; Leroy et al. 2023; Villanueva et al. 2021, 2024). Within individual galaxies, we observe moderate internal variations, with denser regions showing shorter depletion times, but no evidence of extreme starburst-like behaviour, with depletion times above 0.25 Gyr throughout.

Applying the $W_{[CII]}$ prescription fundamentally alters the picture. The derived slope steepens to $\beta = 1.75 \pm 0.14$, and
the sample bridges the gap between the low-redshift CO-based
locus and the regime occupied by the $z = 2.6$ starburst of
Sharon et al. (2013). Several regions within our galaxies, partic-
ularly in CRISTAL_09, CRISTAL_13, and CRISTAL_24, over-
lap with these resolved starburst measurements. This change is
accompanied by a pronounced decrease in depletion times across
the most [CII]-emitting regions, with values reaching as low as
0.1 Gyr, characteristic of starburst systems, in the brightest [CII]
regions. Importantly, this shift towards shorter depletion times
is not restricted to isolated, compact star-forming regions, but is
observed globally across the galaxies. While increasing spatial
resolution is generally expected to broaden the dynamic range of
depletion times and reveal localised starburst-like environments,
we do not observe such a trend here. Instead, the $W_{[CII]}$ conver-
sion produces a systematic reduction in depletion times through-
out all regions, compressing the inferred gas mass range. This
raises the possibility that, under certain conversion prescriptions,
main-sequence galaxies could exhibit global properties reminis-
cent of starbursts, challenging the classical separation between
main-sequence and starburst systems.

Recent results from the THESAN-ZOOM simulations (Shen et al. 2025) show a similar trend towards shorter depletion timescales for $\Sigma_{\text{gas}} > 100 M_{\odot} \text{pc}^{-2}$. In their sample of galaxies with stellar masses up to $10^{11} M_{\odot}$ at $z \sim 3$ and $10^{9.6} M_{\odot}$ at $z \sim 5$, resolved depletion times for neutral gas (HI+H$_2$) can reach as low as 0.1 Gyr in the densest regions (see their Fig. 8). This is consistent with the short depletion times inferred in our analysis.

Additionally, we colour-code our sample by redshift in Figure 5 and find no clear redshift dependence of the relative position of the galaxies within the KS plane. Under the constant $\alpha_{[CII]} = 31 M_{\odot} L_{\odot}^{-1}$, high-redshift main-sequence galaxies exhibit KS behaviour remarkably similar to their low-redshift counterparts, as traced by CO. In contrast, the $W_{[CII]}$ approach reveals a population of short-depletion-time regions that bridge the loci of main-sequence galaxies and starburst regimes, demonstrating the strong sensitivity of inferred gas properties, and thus depletion time estimates, to the adopted [CII]-to-gas calibration. As a result, our ability to confidently study the KS relation at high redshift is compromised.

5. Discussion

To further investigate the possible physical drivers of the derived KS relations, we constructed a diagnostic plot shown in Fig. 6 encoding AGN activity, stellar mass, and merger status for our sample by changing symbol type, size, and colour, respectively. The AGN identification follows Ren et al. (in prep.), while the merger status is classified according to Herrera-Camus et al. (2025). Additionally, environmental classification distinguishes between field galaxies and those residing in the $z \sim 4.57$ protocluster structure (CRISTAL_13, CRISTAL_15, and VC875; see Lemaux et al. 2018; Staab et al. 2024). Visual inspection of the KS diagram reveals no clear systematic offset or trend associated with AGN presence, merger stage, stellar mass, or environment across the KS plane. Within the dynamic range probed by our study, none of these properties seems to have a significant impact on the star formation. Nevertheless, we caution that the uncertainties and limited sample size may reduce sensitivity to subtle or systematic effects.

We also show KS relations for the most massive 25% of
galaxies at $z = 4$ in both the NewHorizon (Dubois et al. 2021;
Kraljic et al. 2024) and NewCluster (Han et al. 2025) simula-
tions, presented in their respective effective gas surface density
regimes. The main difference between the two simulations is

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

**Fig. 6.** Similar to Fig. 5. This time, symbols indicate integrated measurements for individual galaxies, colour-coded by interaction state: light blue for major mergers, light yellow for interacting systems, and orange for isolated galaxies, following the classification of Herrera-Camus et al. (2025). Symbol size encodes the total stellar mass and symbol shape denotes AGN identification: circles for non-AGNs, squares for AGN candidates, and triangles for robust AGNs as per Ren et al. (in prep.). The solid blue line shows the $z = 4$ relation from Kraljic et al. (2024), based on the 25% most massive galaxies in the NewHorizon simulation, while the solid red line indicates the New-Cluster analogue (Han et al. 2025). Both relations are plotted over their effective gas mass ranges for direct comparison with our sample.

environmental: NewCluster targets a high-density proto-cluster
region, while NewHorizon samples an average-density field.
Both use similar physical models for star formation and feed-
back. In the constant $\alpha_{[CII]}$ case, our observed sample aligns well
with predictions from the New Horizon simulation, though this
simulation primarily probes lower stellar masses ($10^{8.0-9.5} M_{\odot}$)
than our galaxies ($10^{9.5-10.6} M_{\odot}$). In comparison, the most mas-
sive 25% of galaxies in the NewCluster simulation span a similar
stellar mass range ($M_* = 10^{9.5-10.5} M_{\odot}$). We find that employ-
ing a density-dependent $W_{[CII]}$ conversion factor aligns with their
results. This compatibility could support the fact that the stellar
mass affects the slope of the KS relation for high-redshift galax-
ies, specifically a steepening of the slope as the total stellar mass
of the population increases, as is hinted at by the $z = 2-3$ rela-
tions presented in Kraljic et al. (2024).

The contrasting behaviours obtained assuming the two differ-
ent gas-conversion prescriptions, and their respective alignments
with different simulations, show the need for reliable [CII]-to-
gas mass conversion methods that account for the ISM condi-
tions prevailing at a given redshift and in a specific population of
galaxies. Even though we cannot fully exclude it at this stage, the
lack of correlation with AGN or merger status further points to
intrinsic ISM properties, such as density, turbulence, and metal-
licity, as the dominant factors shaping the [CII]-to-gas calibra-
tion and, by extension, the inferred star formation law.

These results highlight the limitations of adopting a uni-
versal [CII]-to-gas conversion factor and motivate the develop-
ment of physically motivated, simulation-calibrated or locally
calibrated prescriptions that account for galaxy mass, environ-
ment, and ISM state. Joint analysis of resolved observations and
simulations, as is demonstrated here, is essential for building an accurate picture of the star formation mechanisms in the high-redshift Universe.

## 6. Conclusions

We present a spatially resolved analysis of the [CII]–SFR and KS relations in a sample of 13 main-sequence galaxies at redshifts of $4 < z < 6$, combining ALMA, JWST, and HST data. Our pixel-by-pixel SED modelling approach applied to resolution-homogenised multi-wavelength datasets allows us to probe the physical conditions of star formation down to correlated kiloparsec scales in the early Universe.

We have developed a new comprehensive statistical framework that explicitly accounts for correlated noise and pixel covariance introduced by resolution homogenisation. By combining noise sampling, source bootstrapping, and hybrid resampling techniques, we reliably estimate the intrinsic scatter and slope of the spatially resolved [CII]-SFR relations. Comparison with classical fitting methods demonstrates that neglecting pixel correlations leads to underestimated uncertainties. This shows the relevance of our approach for robust parameter inference.

Our analysis reveals that the resolved [CII]–SFR relation exhibits a slope of $0.87 \pm 0.15$ and an intrinsic scatter of $0.19 \pm 0.03$ dex. These values are both shallower and tighter than those reported in many studies of the local Universe, suggesting a possible saturation of the SFR emission at high [CII] surface densities, population wise, while confirming a strong correlation between [CII] and SFR at $z \sim 5$. This tight relation supports the use of [CII] as a reliable tracer of star formation in early main-sequence galaxies.

In contrast, the resolved KS relation is highly sensitive to the assumption of the [CII]-to-gas conversion factor. Using a constant $\alpha_{[CII]}$ (Zanella et al. 2018) yields a slope of 0.87 and depletion timescales of roughly 0.5 to 1 Gyr, consistent with expectations for main-sequence galaxies. However, adopting the $\Sigma_{[CII]}$-dependent $W_{[CII]}$ prescription (Vallini et al. 2025) results in a much steeper slope of 1.75. The regions with our lowest probed gas surface densities still behave as the disc sequence with depletion timescales of the order of 0.5 Gyr, but high-density regions exhibit shorter depletion times ($<0.1$ Gyr), placing them in the starburst regime.

Notably, these two scenarios bracket the predictions from recent simulations: the constant $\alpha_{[CII]}$ case is in broad agreement with NewHorizon results at $z = 4$ and low mass ($10^{8.0-9.5} M_{\odot}$), while the $W_{[CII]}$ prescription produces KS slopes and depletion times compatible with the galaxies seen in NewCluster in a similar mass regime as the one of our sample ($10^{9.5-10.5} M_{\odot}$). The range of predicted depletion times produced by these two conversion factors points to a fundamental limitation in current gas mass estimates from [CII] and illustrates the need for physically motivated conversion factors that capture the evolving ISM conditions, metallicity, and dynamical processes at high redshift.

**Acknowledgements.** Cédric Accard acknowledges that this work of the Interdisciplinary Thematic Institute IRMIA++, as part of the ITI 2021-2028 programme of the University of Strasbourg, CNRS and Inserm, was supported by IdEx Unistra (ANR-10-IDEX-0002), and by SFRI-STRAT'US project (ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program. This work was supported by the « action thématique » Cosmology-Galaxies (ATCG) of the CNRS/INSU PN Astro. This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Université Côte d'Azur under reference number ANR-15-IDEX-01. M.Bo. acknowledges support from the ANID BASAL project FB210003. L.V. acknowledges support from the INAF Minigrant "RISE:

Resolving the ISM and Star formation in the Epoch of Reionization" (PI: Vallini, Ob. Fu. 1.05.24.07.01). M.A. is supported by FONDECYT grant number 1252054, and gratefully acknowledges support from ANID Basal Project FB210003 and ANID MILENIO NCN2024_112. E.d.C. acknowledges support from the Australian Research Council through project DP240100589. A.F. is partly supported by the ERC Advanced Grant INTERSTELLAR H2020/740120, and by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics. R.H.-C. thanks the Max Planck Society for support under the Partner Group project "The Baryon Cycle in Galaxies" between the Max Planck for Extraterrestrial Physics and the Universidad de Concepción. R.H.-C. also gratefully acknowledge financial support from ANID - MILENIO - NCN2024_112 and ANID BASAL FB210003. J.M. gratefully acknowledges support from ANID MILENIO NCN2024_112. A.N., P.S. acknowledge support from the Narodowe Centrum Nauki (NCN), Poland, through the SONATA BIS grant UMO-2020/38/E/ST9/00077. M.P. acknowledges financial support from the project "LEGO - Reconstructing the building blocks of the Galaxy by chemical tagging" granted by the Italian MUR through contract PRIN2022LLP8PKT_001. M.Re. acknowledges support from project PID2023-150178NB-I00 financed by MCIU/AEI/10.13039/501100011033. V.V. acknowledges support from the ANID BASAL project FB210003 and from ANID - MILENIO - NCN2024_112. W.W. acknowledges the grant support through JWST programs. Support for programs JWST-GO-03045 and JWST-GO-03950 were provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. S.K.Y. acknowledges support from the Korean National Research Foundation (RS-2025-00514475 and RS-2022-NR070872). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00428L, ADS/JAO.ALMA#2019.1.00226.S, ADS/JAO.ALMA#2021.1.00280.L, ADS/JAO.ALMA#2022.1.01118.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope and NASA/ESA Hubble Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST and NAS 1432 5-26555 for HST. NewCluster was granted access to the HPC resources of KISTI under the allocations KSC-2022-CRE-0344 and KSC-2023-CRE-0343, and of GENCI under the allocation A0150414625. We thank Diana Ismail for sharing her data compilation of NewCluster. We acknowledge the use of Astropy (Astropy Collaboration 2018), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020).

## References

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

* Béthermin, M., Accard, C., Guillaume, C., et al. 2023, A&A, 680, L8

* Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846

* Boquien, M., Calzetti, D., Aalto, S., et al. 2015, A&A, 578, A8

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

* Boquien, M., Buat, V., Burgarella, D., et al. 2022, A&A, 663, A50

* Buat, V., Ciesla, L., Boquien, M., Malek, K., & Burgarella, D. 2019, A&A, 632, A79

* Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582

* Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

* Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

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

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

* Ciesla, L., Elbaz, D., & Fensch, J. 2017, A&A, 608, A41

* Cochrane, R. K., Hayward, C. C., & Anglés-Alcázar, D. 2022, ApJ, 939, L27

* Condon, J. J. 1997, PASP, 109, 166

* Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

* da Cunha, E., & Charlot, S. 2011, MAGPHYS: Multi-wavelength Analysis of Galaxy Physical Properties, Astrophysics Source Code Library [record ascl:1106.0100]

* De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62

* de los Reyes, M. A. C., & Kennicutt, R. C. 2019, ApJ, 872, 16

* Dessauges-Zavadsky, M., Ginolfi, M., Pozzi, F., et al. 2020, A&A, 643, A5

* D'Eugenio, C., Daddi, E., Liu, D., & Gobat, R. 2023, A&A, 678, L9

* Devereaux, T., Cassata, P., Ibar, E., et al. 2024, A&A, 686, A156

* Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

* Faisst, A. L., Capak, P. L., Yan, L., et al. 2017, ApJ, 847, 21

* Faisst, A. L., Yan, L., Béthermin, M., et al. 2022, Universe, 8, 314

* Ferrara, A., Vallini, L., Pallottini, A., et al. 2019, MNRAS, 489, 1
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Freundlich, J. 2024, Fundam. Plasma Phys., 11, 100059

Freundlich, J., Combes, F., Tacconi, L. J., et al. 2013, A&A, 553, A130

Freundlich, J., Combes, F., Tacconi, L. J., et al. 2019, A&A, 622, A105

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Gordon, K., Boyer, M., Sloan, G., Muzerolle, J., & Volk, K. 2018, An End-to-End Calibration Error Budget for the JWST Science Instruments, Technical Report JWST-STScI-001007

Han, S., Yi, S. K., Dubois, Y., et al. 2025, A&A, submitted [arXiv:2507.06301]

Hashimoto, T., Álvarez Márquez, J., Fudamoto, Y., et al. 2023, ApJ, 955, L2

Hasinger, G., Capak, P., Salvato, M., et al. 2018, ApJ, 858, 77

Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, AJ, 800, 1

Herrera-Camus, R., González-López, J., Förster Schreiber, N., et al. 2025, A&A, 699, A80

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jiménez-Donaire, M. J., Brown, T., Wilson, C. D., et al. 2023, A&A, 671, A3

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

Kennicutt, R. C., Jr 1998, ApJ, 498, 541

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

Khusanova, Y., Bethermin, M., Le Fèvre, O., et al. 2021, A&A, 649, A152

Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Kraljic, K., Renaud, F., Dubois, Y., et al. 2024, A&A, 682, A50

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1

Lemaux, B. C., Le Fèvre, O., Cucciati, O., et al. 2018, A&A, 615, A77

Leroy, A. K., Walter, F., Sandstrom, K., et al. 2013, AJ, 146, 19

Leroy, A. K., Bolatto, A. D., Sandstrom, K., et al. 2023, ApJ, 944, L10

Li, J., Da Cunha, E., González-López, J., et al. 2024, ApJ, 976, 70

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

Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141

Mitsuhashi, I., Tadaki, K.-I., Ikeda, R., et al. 2024, A&A, 690, A197

Molina, J., Ibar, E., Smail, I., et al. 2019, MNRAS, 487, 4856

Nagy, D., Dessauges-Zavadsky, M., Messa, M., et al. 2023, A&A, 678, A183

Noll, S., Pierini, D., Cimatti, A., et al. 2009, A&A, 499, 69

Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621

Pessa, I., Schinnerer, E., Belfiore, F., et al. 2021, A&A, 650, A134

Romano, M., Morselli, L., Cassata, P., et al. 2022, A&A, 660, A14

Saintonge, A., Tacconi, L. J., Fabello, S., et al. 2012, ApJ, 758, 73

Schäerer, D., Ginolfi, M., Béthermin, M., et al. 2020, A&A, 643, A3

Schmidt, M. 1959, ApJ, 129, 243

Schruba, A., Leroy, A. K., Walter, F., et al. 2011, AJ, 142, 37

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Sharon, C. E., Baker, A. J., Harris, A. I., & Thomson, A. P. 2013, ApJ, 765, 6

Shen, X., Kannan, R., Puchwein, E., et al. 2025, MNRAS, submitted, [arXiv:2503.01949]

Solomon, P. M., Downes, D., & Radford, S. J. E. 1992, ApJ, 387, L55

Staab, P., Lemaux, B. C., Forrest, B., et al. 2024, MNRAS, 528, 6934

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Vallini, L., Gallerani, S., Ferrara, A., Pallottini, A., & Yue, B. 2015, ApJ, 813, 36

Vallini, L., Witstok, J., Sommovigo, L., et al. 2024, MNRAS, 527, 10

Vallini, L., Pallottini, A., Kohandel, M., et al. 2025, A&A, 700, A117

Vanderhoof, B. N., Faisst, A. L., Shen, L., et al. 2022, MNRAS, 511, 1303

Veraldi, E., Vallini, L., Pozzi, F., et al. 2025, A&A, 693, A34

Villanueva, V., Bolatto, A., Vogel, S., et al. 2021, ApJ, 923, 60

Villanueva, V., Bolatto, A. D., Vogel, S. N., et al. 2024, ApJ, 962, 88

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Vizgan, D., Greve, T. R., Olsen, K. P., et al. 2022, ApJ, 929, 92

Wang, T.-M., Magnelli, B., Schinnerer, E., et al. 2022, A&A, 660, A142

Wang, B., Leja, J., Labbé, I., et al. 2023, ApJS, 270, 12

Wolfire, M. G., Vallini, L., & Chevance, M. 2022, ARA&A, 60, 247

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976

Zavala, J. A., Bakx, T., Mitsuhashi, I., et al. 2024, ApJ, 977, L9

3 Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

4 INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, 40129 Bologna, Italy

5 Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

6 Millenium Nucleus for Galaxies (MINGAL), Concepción, Chile

7 Departamento de Astronomía, Universidad de Concepción, Barrio Universitario, Concepción, Chile

8 Dipartimento di Fisica e Astronomia, Università di Padova, Vicolo dell'Osservatorio 3, Padova, Italy

9 INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

10 International Centre for Radio Astronomy Research (ICRAR), The University of Western Australia, M468, 35 Stirling Highway, Crawley, WA 6009, Australia

11 ARC Center of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Canberra, Australia

12 Sterrenkundig Observatorium, Ghent University, Krijgslaan 281 - S9, B9000 Ghent, Belgium

13 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland

14 Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, 75014 Paris, France

15 IPAC, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA

16 Center for Frontier Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan

17 Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, 50019 Sesto Fiorentino, Firenze, Italy

18 INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

19 Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan

20 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

21 Gemini Observatory, NSF NOIRLab, 670 N. A'ohoku Place, Hilo, Hawaii 96720, USA

22 Department of Physics and Astronomy, University of California, Davis, One Shields Ave., Davis, CA 95616, USA

23 Department of Physics and Astronomy and George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, 4242 TAMU, College Station, TX 77843-4242, USA

24 Department of Physics and Astronomy, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA

25 Instituto de Física y Astronomía, Universidad de Valparaíso, Avda. Gran Bretaña 1111, Valparaíso, Chile

26 INAF – Osservatorio astronomico d’Abruzzo, Via Maggini SNC, 64100 Teramo, Italy

27 National Centre for Nuclear Research, ul. Pasteura 7, 02-093 Warsaw, Poland

28 Dipartimento di Fisica e Astronomia “Augusto Righi”, Alma Mater Studiorum, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

29 Dept. Fisica Teorica y del Cosmos, Universidad de Granada, Granada, Spain

30 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

31 Scuola Internazionale Superiore Studi Avanzati (SISSA), Physics Area, Via Bonomea 265, 34136 Trieste, Italy

32 Departamento de Astronomía, Universidad de Concepción, Barrio Universitario, Concepción, Chile

33 Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, Republic of Korea

34 Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan

<footnote n="1">¹ Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="2">² Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06000 Nice, France</footnote>
## Appendix A: Homogenised data

The resolution-homogenised maps of VC875 is shown in Fig. A.1 as an example.

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

**Fig. A.1.** Resolution homogenised version of Fig. 1, with the ALMA resolution as the target. The first row presents, from left to right, HST/ACS F814W (rest-frame 145 nm), JWST/NIRCam F115W (207 nm), F150W (270 nm), and F277W (499 nm). The second row shows JWST/F444W (800 nm), ALMA dust continuum emission around 156.9 µm, and ALMA [CII] 158 µm emission. The black contours are the (3+2k)σ levels ($k \ge 0$) of the [CII] luminosity map, while the ellipses in the bottom left corners all panels indicate the synthesised ALMA beam. These images illustrate the filter coverage available across all galaxies in our sample.

## Appendix B: SED modelling parameters

<table><caption>**Table B.1.** Parameter values for different modules used for the CIGALE runs.</caption><thead><tr><th>Module</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td rowspan="4">sfhdelayedbq</td><td>age_main (Myr)</td><td>(Universe age at source redshift) - 250 Myr</td></tr><tr><td>tau_main (Myr)</td><td>100, 300, 500, 700, 900, 1100, 1300</td></tr><tr><td>age_bq (Myr)</td><td>0, 5, 10, 25, 50, 100, 150, 200, 250, 300, 350, 400, 450</td></tr><tr><td>r_sfr</td><td>33 log-spaced values from 0 to 900.0</td></tr><tr><td rowspan="2">bc03</td><td>imf</td><td>1 (Chabrier 2003)</td></tr><tr><td>metallicity</td><td>0.008</td></tr><tr><td>nebular</td><td>logU</td><td>-2.0</td></tr><tr><td rowspan="4">dustatt_modified_CF00</td><td>Av_ISM</td><td>101 linearly spaced values from 0.3 to 5</td></tr><tr><td>mu</td><td>0.5</td></tr><tr><td>slope_ISM</td><td>-0.7</td></tr><tr><td>slpoe_BC</td><td>-0.7</td></tr><tr><td rowspan="2">dustatt_calzleit<br>(tested but not used)</td><td>E_BVs_young</td><td>100 linearly spaced values from 0.005 to 0.500</td></tr><tr><td>powerlaw_slope</td><td>0 or [-2.0, 0.5]</td></tr><tr><td>redshifting</td><td>redshift</td><td>[CII] spectroscopic redshift from ALPINE (Béthermin et al. 2020)</td></tr></tbody></table>

**Notes.** The values for `dustatt_modified_CF00` are from Buat et al. (2019) whereas the ones for `dustatt_calzleit` come from Boquien et al. (2022). Any parameter of a given module not shown here or not discussed in Sect. 3.1 has been left to the default value.
# Appendix C: Assessing the validity of the SED modelling

We assess the impact of different parametrisations in the CIGALE SED modelling process on SFR estimates. To evaluate the reliability of each model configuration, we compare the SFR values derived by CIGALE ($SFR_{SED}$) with those independently inferred from combined UV and far-infrared (FIR) measurements, following the methodology of Béthermin et al. (2023). Figure C.1 presents this comparison for VC875 as an example with each panel corresponding to a different set of CIGALE parameters. The one-to-one relation (dashed line) and the ±50% deviation lines (dotted) are included to facilitate a direct visual assessment of model performance.

It is important to emphasise that agreement between SED-based and UV+FIR SFR estimates does not necessarily guarantee the absence of systematic uncertainties or biases, since both methods may be subject to common assumptions regarding dust attenuation laws, SFHs, or IMFs.

The top panel shows results from our fiducial model, which employs the Charlot & Fall (2000) dust attenuation law combined with a delayed SFH including a recent episode of star formation variation. This configuration achieves the closest agreement between SFR estimators, with most data points clustering near the one-to-one line.

When the recent SFH variation is excluded (second panel), the SED-based SFR estimates are systematically lower, with the majority of pixels lying above the one-to-one relation, showing a deterioration in fit quality.

Excluding the dust luminosity constraint from the modelling (third panel) results in systematic overestimation of SED-based SFRs relative to the UV+FIR values, with many points falling well beyond the ~50% deviation line. In this situation, the error bars increase significantly, reflecting greater uncertainties in the fits. This highlights the critical role of including infrared constraints to properly model dust attenuation and derive reliable SFRs.

The fourth panel presents results obtained using the Calzetti et al. (1994) dust attenuation law. Across the sample, differences relative to the fiducial model are minimal. However, for CRISTAL_24 (see Fig. C.2), we observe the largest discrepancy between attenuation prescriptions: the Charlot & Fall (2000) law (top panel) achieves the best agreement between SFR estimators, with most data points lying close to the one-to-one line, while the fiducial model (centre panel) systematically underestimates SFRs and yields a lower fraction of well-fitted pixels. We verified that this difference is not attributable to the choice of the $\delta$ parameter in the Calzetti law, as letting $\delta$ vary produces identical results (bottom panel). Given the improved consistency for CRISTAL_24 and the more physically motivated nature of the Charlot & Fall (2000) law, we adopt it as our default attenuation model for the entire sample.

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

**Fig. C.1.** Comparison between SFR derived from SED modelling (SFRSED) and those measured from IR+UV observations (SFRIR + UV) for VC875, shown for different model parametrisations. From top to bottom: (first) the fiducial model using the Charlot & Fall (2000) dust attenuation law and a delayed SFH with a recent variation episode; (second) the same, but without the recent variation episode; (third) the fiducial model without ALMA dust luminosity constraints; (fourth) the Calzetti et al. (1994) attenuation law model. Blue points represent individual pixels, with error bars corresponding to the 16th and 84th percentiles of the SED-derived SFR distribution (see Sect. 3.1). Red arrows indicate upper limits. Black points show region-averaged SFR values computed using the B23 method. The dashed line indicates the one-to-one relation, and dotted lines mark the ±50% range.
<figure><img src="image_10.png" /><figcaption>**Fig. C.2.** Same as Fig. C.1, but for CRISTAL_24. From top to bottom: Results using the Charlot & Fall (2000) dust attenuation model, results using the modified Calzetti et al. (1994) model with $\delta = 0$, same but letting the $\delta$ parameter vary between -2.0 and 0.5</figcaption></figure>

## Appendix D: Validation of the likelihood model

To evaluate the accuracy and robustness of our fitting methodology, we perform an end-to-end simulation designed to replicate the observational conditions and intrinsic properties of our dataset. Specifically, we generate 1000 independent mock realisations. For each of them, we create first gas surface density maps composed of a random number of elliptical Gaussian clumps with varying amplitudes and positions. We then convolve this map with the ALMA beam. Using a fixed star formation relation characterised by a known slope, intercept, and intrinsic scatter, we derive the corresponding SFR surface density maps on a pixel-by-pixel basis. We then simulated the intrinsic scatter by multiplying this map by a map containing log-normal fluctuations at the scale of the ALMA beam with dispersion corresponding to the intrinsic scatter and a mean value of unity. This step is important, since we measure the intrinsic scatter at the scale of the ALMA beam, and applying a pixel-based scatter would lead to a smaller scatter after smoothing the maps by the beam.

We then generate instrumental noise maps scaled to the noise properties of our ALMA data. We first convolve these noise maps to the same resolution, rescale their standard deviation, and added them to the gas maps. For the SFR maps, we applied the typical log-normal errors found by CIGALE for the real

<figure><img src="image_11.png" /><figcaption>**Fig. D.1.** Distributions of the fitted parameters obtained from 1000 simulated datasets (see Appendix D). From top to bottom: intercept, slope, and intrinsic scatter. The blue histograms show the results from our likelihood-based fitting method, with the solid and dotted blue lines indicating the mean and the $\pm 1\sigma$ intervals. The true input value used to generate the simulations is shown as a dashed green line. Each panel is annotated with the mean and standard deviation of the fitted parameter and the true input value in the top left corner.</figcaption></figure>
data. This process was repeated 13 times to produce a mock
population comparable in size and complexity to our observed
sample. Each mock dataset is subsequently analysed using our
likelihood-based fitting procedure of Sect. 3.2, identical to what
used for the real observations. We recover maximum likelihood
solution of the slope, intercept, and intrinsic scatter parameters.
The distributions of the fitted parameters across 1000 mock pop-
ulation realisations were then examined to assess the fidelity of
the fitting method, see Fig. D.1.

We find that, for each parameter, the true input value lies
well within the 1σ confidence interval of the recovered distribu-
tion. Still, small systematic offsets between the mean fitted val-
ues and the true parameters are observed. The intercept is over-
estimated by approximately 0.15 dex relative to the true value,
while the slope and intrinsic scatter are underestimated by about
0.017 and 0.06 dex, respectively. These relative differences cor-
respond to roughly 0.45, 0.39, and 0.43 times the respective 1σ
uncertainties, indicating that the deviations are minor compared
to the statistical errors. We note that the 1σ intervals from the
mock ensemble are systematically smaller than those from the
real data. It could be due to a better S/N or the larger number of
beams with a [CII] signal in our simulation. We therefore con-
clude that our approach does not systematically bias the esti-
mates of the key parameters of the star formation relation.

Appendix E: Alternative Fitting Methods

In addition to this primary analysis, we use two other fitting
methods to assess the validity of the likelihood model and its
maximisation approach: the Bayesian linear regression package
`linmix` and the MCMC sampler `emcee` (Foreman-Mackey et al.
2013).

For the MCMC analysis, we configure 32 walkers with 1,000
steps per chain, discarding the first 100 steps as burn-in and thin-
ning chains by a factor of 20 to mitigate autocorrelation. To
take into account the correlation between pixels, we normalise
the likelihood from Eq. 9 by the number of independent pixels
per beam. This approach aims to approximate the number of
independent points per resolution element. However, since beam
smoothing may effectively yield twice as many statistically inde-
pendent noise samples as naively expected (see Condon 1997;
Condon et al. 1998), we thus divide the likelihood by twice
the number of pixels per beam. The impact of this normalisa-
tion is illustrated by the following MCMC fitting results. With-
out any likelihood normalisation, we obtain $\beta = 0.79_{-0.02}^{+0.01}$ and
$\sigma_i = 0.17_{-0.00}^{+0.00}$. When normalising the likelihood by twice the
number of pixels per beam, the results broaden to $\beta = 0.78_{-0.08}^{+0.09}$
and $\sigma_i = 0.17_{-0.02}^{+0.02}$. This demonstrates that accounting for the
effective number of independent data points leads to larger and
more realistic uncertainty estimates on the fitted parameters.

The *linmix* method uses Bayesian regression to account for
measurement errors in both variables and intrinsic scatter, pro-
viding robust estimates of the slope, intercept, and intrinsic scat-
ter of the studied [CII]-SFR relation. However, like in classical
MCMC approach, it assumes uncorrelated pixels, ignoring spa-
tial correlations within the beam and thus underestimate uncer-
tainties in the presence of significant beam smoothing.

While both the *linmix* and MCMC fitting procedures pro-
vide estimates of the relation parameters and their uncertainties,
they do not fully account for all sources of correlation, noise
propagation, or population sampling effects present in our data.
In particular, neither method can rigorously handle the spatial
correlations and noise structure introduced by beam smooth-

ing, nor the potential biases arising from the small sample size.
For these reasons, we consider our primary resampling-based
approaches described in Sect. 3.2 and Sect. 3.3 to provide more
reliable and conservative estimates of the uncertainties and best-
fit parameters for the [CII]–SFR relation.
## Appendix F: Other sources property maps

As in Fig. 3, we show here example of spatially resolved physical parameters maps obtained with the CIGALE SED modelling for the other sources in the sample in Figs. F.1, F.2 and F.3. In each case, we have: stellar mass surface density ($M_\odot$ kpc$^{-2}$), SFR surface density averaged over 10 Myrs ($M_\odot$ yr$^{-1}$ kpc$^{-2}$), mass-weighted age (Myr) and attenuation in FUV rest-frame (mag). Results are shown with the final selection used for the [CII]-SFR and KS relation analysis, that is pixels with S/N ≥ 2 in at least 3 photometric bands and where $\Sigma_{[CII]} > 5\sigma_{[CII]}$. Across the sample, we find that the peaks of stellar mass and SFR are generally aligned. In contrast, regions of highest SFR tend to coincide with diminished FUV attenuation, and the mass-weighted age maps reveal that the oldest stellar populations are preferentially located in regions of highest stellar mass surface density.

<figure><img src="image_12.png" /><figcaption>Fig. F.1. Similar to Fig. 3 for CRISTAL_01, CRISTAL_02, CRISTAL_03, and CRISTAL_04 (top to bottom)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. F.2. Similar to Fig. 3 for CRISTAL_06, CRISTAL_07, CRISTAL_09, CRISTAL_11 and CRISTAL_13 (top to bottom)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. F.3. Similar to Fig. 3 for CRISTAL_15, CRISTAL_19, and CRISTAL_24 (top to bottom)</figcaption></figure>