11institutetext: Departamento de Astronomía, Universidad de Concepción, Barrio
Universitario, Concepción, Chile
email<EMAIL_ADDRESS>22institutetext: Instituto de Astrofísica,
Facultad de Física, Pontificia Universidad Católica de Chile, Santiago
7820436, Chile 33institutetext: Instituto de Estudios Astrofísicos, Facultad
de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441,
Santiago 8370191, Chile 44institutetext: National Radio Astronomy Observatory,
520 Edgemont Road, Charlottesville, VA 22903, USA 55institutetext:
INAF–Istituto di Radioastronomia, Via Gobetti 101, 40129 Bologna, Italy
66institutetext: Chemistry Department, Sapienza University of Rome, P.le A.
Moro, 00185 Rome, Italy 77institutetext: Jodrell Bank Centre for Astrophysics,
Department of Physics and Astronomy, School of Natural Sciences, The
University of Manchester, Manchester, M13 9PL, UK 88institutetext:
International Centre for Radio Astronomy Research, ICRAR M468, 35 Stirling
Hwy, Crawley 6009, Western Australia 99institutetext: Sterrenkundig
Observatorium, Ghent University, Krijgslaan 281 – S9, B9000 Ghent, Belgium
1010institutetext: Department of Physics & Astronomy, University College
London, Gower Street, London WC1E 6BT, UK 1111institutetext: Institute of
Astrophysics, Foundation for Research and Technology–Hellas (FORTH),
Heraklion, GR-70013, Greece 1212institutetext: School of Sciences, European
University Cyprus, Diogenes street, Engomi, 1516 Nicosia, Cyprus
1313institutetext: Max-Planck-Institut fuer extraterrestrische Physik,
Giessenbachstrasse 1, 85748 Garching, Germany 1414institutetext: Hiroshima
Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-
Hiroshima, 739-8526 Hiroshima, Japan 1515institutetext: Department of
Astronomy, School of Science, SOKENDAI (The Graduate University for Advanced
Studies), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 1616institutetext:
Department of Astronomy, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan 1717institutetext: National Astronomical Observatory of Japan,
2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 1818institutetext: Max Planck
Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany
1919institutetext: Departamento Física Teórica y del Cosmos, Universidad de
Granada, E-18071 Granada, Spain 2020institutetext: Instituto Universitario
Carlos I de Física Teórica y Computacional, Universidad de Granada, E-18071
Granada, Spain 2121institutetext: 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 2222institutetext: Dipartimento di Fisica e Astronomia, Universita di
Bologna, via Gobetti 93/2, I-40129 Bologna, Italy 2323institutetext: INAF -
OAS, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via
Gobetti 93/3, I-40129 Bologna, Italy 2424institutetext: Department of Physics
and Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA 15260,
USA 2525institutetext: Faculty of Engineering, Hokkai-Gakuen University,
Toyohira-ku, Sapporo 062-8605, Japan 2626institutetext: Kavli Institute for
Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
2727institutetext: Cavendish Laboratory, University of Cambridge, 19 JJ
Thomson Avenue, Cambridge CB3 0HE, UK 2828institutetext: Scuola Normale
Superiore, Piazza dei Cavalieri 7, I- 50126 Pisa, Italy 2929institutetext: Las
Campanas Observatory, Carnegie Institution of Washington, Casilla 601, La
Serena, Chile;

# The ALMA-CRISTAL survey: Dust temperature and physical conditions of the
interstellar medium in a typical galaxy at $z=5.66$

V. Villanueva , ALMA-ANID Postdoctoral Fellow11 R. Herrera-Camus 11 J.
González-López 222828 M. Aravena 33 R. J. Assef 33 Mauricio Baeza-Garay 11 L.
Barcos-Muñoz 44 S. Bovino 5566 R. A. A. Bowler 77 E. da Cunha 88 I. De Looze
991010 T. Diaz-Santos 11111212 A. Ferrara 2828 N. Foerster-Schreiber 1313 H.
Algera 1414 R. Ikeda 15151717 M. Killi 33 I. Mitsuhashi 16161717 T. Naab 1818
M. Relano 19192020 J. Spilker 2121 M. Solimano 33 M. Palla 22222323 S.H. Price
2424 A. Posses 33 K. Tadaki 2525 K. Telikova 33 H. Übler 26262727

(Received July XX, 2024; accepted October XX, 2024)

We present new $\lambda_{\rm rest}=77$ $\mu$m dust continuum observations from
the Atacama Large Millimeter/submillimeter Array of HZ10 (CRISTAL-22), a dusty
main-sequence galaxy at $z$=5.66 as part of the [CII] Resolved Ism in STar-
forming Alma Large program, CRISTAL. The high angular resolution of the ALMA
Band 7 and new Band 9 data ($\sim{0}^{\prime\prime}.4$) reveals the complex
structure of HZ10, which comprises two main components (HZ10-C and HZ10-W) and
a bridge-like dusty emission between them (the Bridge). Using a modified
blackbody function to model the dust spectral energy distribution (SED), we
constrain the physical conditions of the interstellar medium (ISM) and its
variations among the different components identified in HZ10. We find that
HZ10-W (the more UV-obscured component) has an SED dust temperature of $T_{\rm
SED}$$\sim$51.2$\pm 13.1$ K; this is $\sim$5 K higher (although still
consistent) than that of the central component and previous global estimations
for HZ10. Our new ALMA data allow us to reduce by a factor of $\sim$2.3 the
uncertainties of global $T_{\rm SED}$ measurements compared to previous
studies. HZ10 components have [CII]-to-far-infrared (FIR) luminosity ratios
and FIR surface densities values consistent with those for local starbursts
galaxies. However, HZ10-W shows a lower [CII]/FIR ratio compared to the other
two components (although still within the uncertainties), which may suggest a
harder radiation field destroying polycyclic aromatic hydrocarbon associated
with [CII] emission (e.g., active galactic nuclei or young stellar
populations). While HZ10-C appears to follow the tight IRX-$\beta_{\rm UV}$
relation seen in local UV-selected starburst galaxies and high-$z$ star-
forming galaxies, we find that both HZ10-W and the Bridge depart from this
relation and are well described by dust-screen models with holes in front of a
hard UV radiation field. This suggests that the UV emission (likely from young
stellar populations) is strongly attenuated in the more dusty components of
the HZ10 system.

We present new $\lambda_{\rm rest}=77$ $\mu$m dust continuum observations from
the Atacama Large Millimeter/submillimeter Array of HZ10 (CRISTAL-22), a dusty
main-sequence galaxy at $z$=5.66 as part of the [CII] Resolved Ism in STar-
forming Alma Large program, CRISTAL. The high angular resolution of the ALMA
Band 7 and new Band 9 data ($\sim{0}^{\prime\prime}.4$) reveals the complex
structure of HZ10, which comprises two main components (HZ10-C and HZ10-W) and
a bridge-like dusty emission between them (the Bridge). Using a modified
blackbody function to model the dust spectral energy distribution (SED), we
constrain the physical conditions of the interstellar medium (ISM) and its
variations among the different components identified in HZ10. Our new ALMA
data allow us to reduce by a factor of $\sim$2.3 the uncertainties of global
$T_{\rm SED}$ measurements compared to previous studies. We find that HZ10-W
(the more UV-obscured component) has an SED dust temperature of $T_{\rm
SED}$$\sim$51.5$\pm 13.1$ K; this is $\sim$5 K higher (although still
consistent) than that of the other two components and previous global
estimations for HZ10. HZ10 components have [CII]-to-far-infrared (FIR)
luminosity ratios and FIR surface densities values consistent with those for
local starbursts galaxies. However, HZ10-W shows a lower [CII]/FIR ratio
compared to the other two components (although still within the
uncertainties), which may suggest a harder radiation field destroying
polycyclic aromatic hydrocarbon associated with [CII] emission (e.g., active
galactic nuclei or young stellar populations). While HZ10-C appears to follow
the tight IRX-$\beta_{\rm UV}$ relation seen in local UV-selected starburst
galaxies and high-$z$ star-forming galaxies, we find that both HZ10-W and the
Bridge depart from this relation and are well described by dust-screen models
with holes in front of a hard UV radiation field. This suggests that the UV
emission (likely from young stellar populations) is strongly attenuated in the
more dusty components of the HZ10 system.

###### Key Words.:

galaxies: evolution – galaxies: ISM – submillimeter: galaxies – ISM: lines and
bands

##  1 Introduction

Accurate measurements of dust temperatures ($T_{\rm dust}$) in star-forming
main sequence galaxies are critical to determine their infrared luminosities
($L_{\rm IR}$), star formation rates (SFR), and dust attenuation properties
(e.g., the excess of infrared emission compared to the ultra-violet, UV, or
Balmer decrement, among others), which are fundamental quantities in the
context of galaxy evolution. However, precise estimations of $T_{\rm dust}$
require the detection of the dust continuum emission near the peak of the
infrared spectral energy distribution (SED; e.g., Hodge & da Cunha, 2020; da
Cunha et al., 2021). Although this has been partially achieved at low- and
high-redshift galaxies (e.g., Bakx et al. 2021; Witstok et al. 2022; Akins et
al. 2022; Tsukui et al. 2023; Algera et al. 2024), we still need better
constraints around the peak of the dust SED in individual main-sequence or
typical star-forming main-sequence galaxies at $z\gtrsim 4$.

The advent of the new generation telescopes (e.g., the Karl G. Jansky Very
Large Array, VLA; the NOrthern Extended Millimetre Array, NOEMA; Atacama Large
Millimeter/submillimeter Array, ALMA, James Webb Space Telescope, JWST; among
others) have revolutionized the exploration of the physical properties of
dust, the extragalactic cold neutral/molecular gas, and their close relation
with the star formation activity. Nevertheless, we still lack a consensus
about the disagreement between $T_{\rm dust}$ estimates at low ($z\lesssim
1-2$) and high ($z\gtrsim 4$) redshift. Although dust peak temperature derived
from SED fitting (i.e., the temperature at the wavelength where the SED peaks,
$T_{\rm peak}$) in low-$z$ galaxies show a good agreement with predictions
from models, estimations in star-forming galaxies at high-$z$ find very low
$T_{\rm peak}$ (e.g., compared to predictions by Viero et al. 2022 at
low-$z$). The latter suggests that the physical conditions of dust at high
redshift may differ significantly compared to those in the local Universe. In
addition, while dust temperatures are well constrained by SEDs densely sampled
in wavelength in local galaxies (e.g., Villanueva et al. 2017; Herrera-Camus
et al. 2018b), high-$z$ galaxies usually lack a proper SED coverage (e.g.,
REBELS, Inami et al. 2022; SERENADE, Mitsuhashi et al. 2023a), which
translates a $T_{\rm dust}$ susceptible to severe biases (e.g., Bakx et al.,
2021; Algera et al., 2024). This can produce significant differences in the
derived far IR luminosities ($L_{\rm FIR}$, up to 3 times and more; Bouwens et
al. 2020), which not only affects our understanding of the dust physical
properties in high-$z$ galaxies, but also for the interpretation of their
interstellar medium (ISM) and star formation properties (e.g., Faisst et al.
2020; Herrera-Camus et al. 2021).

Galaxy surveys on the cold neutral gas deepened our knowledge of the star-
formation main-sequence (MS) in the local Universe (e.g., Brinchmann et al.
2004; Whitaker et al. 2012; Cano-Díaz et al. 2016; Saintonge et al. 2016;
Colombo et al. 2020; Villanueva et al. 2024), and at high-$z$ ($z\gtrsim 4$;
e.g., Capak et al. 2015, and the Alma Large Program to INvestigate C+ at Early
times, ALPINE; Le Fèvre et al. 2020). For instance, spectral studies of [CII]
data in high-$z$ galaxies have revealed signatures of high speed outflows
($\sim$400-500 km s-1) with mass-outflows rates comparable to their SFRs
(e.g., Gallerani et al., 2018; Ginolfi et al., 2020; Herrera-Camus et al.,
2021), probably associated to diffuse and extended [CII] components around
galaxies (or [CII] halos; e.g., Fujimoto et al. 2019, 2020; Fudamoto et al.
2022; Pizzati et al. 2020; Solimano et al. 2024). Moreover, while low-$z$,
normal star-forming galaxies show a tight correlation between the [CII]
luminosity and the SFR surface densities (the $\Sigma_{\rm[CII]}$-$\Sigma_{\rm
SFR}$ relation; e.g., De Looze et al. 2014; Herrera-Camus et al. 2015; Lupi et
al. 2019), more extreme star-forming systems depart from the
$L_{\rm[CII]}/L_{\rm IR}-L_{\rm IR}$ relation due to a deficit in their [CII]
content. At low-$z$ (e.g., $z\lesssim 0.2$), these [CII]-deficient galaxies
are typically dusty, dense starbursts characterized by their hard radiation
fields (see Malhotra et al. 2001; Graciá-Carpio et al. 2011; Díaz-Santos et
al. 2017; Herrera-Camus et al. 2018b and reference therein). However, the
modest spatial resolution achieved by most high-$z$ galaxy surveys (typically
$\sim$5-10 kpc at best) does not allow us to verify this effect at the
relevant physical scales ($\sim$1 kpc; e.g., Shibuya et al. 2015). A detailed
characterization of the ISM in high-$z$ galaxies at physical scales comparable
to those at low redshift is therefore crucial for a more comprehensive
understanding of the nature of this effect.

Given that the shape of the dust SED is very sensitive to the temperature
(i.e., $T_{\rm SED}$, which refers to the estimation of the $T_{\rm dust}$
from SED modeling), adopting an average value derived from local galaxies can
potentially underestimate the IR luminosity by a factor of 5 (e.g., Faisst et
al. 2017) or even higher (e.g., Hodge & da Cunha, 2020). The observational
evidence suggests that galaxies dust temperatures are on average warmer at
high-$z$ (e.g., Magdis et al. 2012; Magnelli et al. 2014; Béthermin et al.
2015; Ferrara et al. 2017; Schreiber et al. 2018; Liang et al. 2019; Sommovigo
et al. 2020), which could be the result of either higher star formation
activity and/or lower metal content compared to those at low $z$. More
accurate constraints on $T_{\rm SED}$ estimations are thus necessary for a
better characterization of the IR SEDs and the derivation of the physical
conditions of the ISM, and particularly for the dust, in galaxies at $z>4$.

In order to check how ISM properties (including dust heating) vary on
spatially resolved scales in high-$z$ galaxies, we present a systematic
analysis of the dust and [CII] in HZ10, a main-sequence galaxy at $z\approx
5.66$ (see Table 1 for more details), as part of the [CII] Resolved Ism in
STar-forming galaxies with ALma survey, CRISTAL (Mitsuhashi et al. 2023b;
Solimano et al. 2024; Posses et al. 2024; Herrera-Camus et al. in prep.).
Based on new ALMA Band 9+7 data (77 and 158 $\mu$m dust continuum) and earlier
ALMA Band 8+6 (110 and 198 $\mu$m dust continuum) observations, we cover the
region close to the peak of the dust SED. This allows us to obtain $T_{\rm
SED}$ measurements for HZ10 at kpc-scales, significantly improving the
accuracy from the previous studies (derivations). The paper is organized as
follows: Section 2 presents the main features of HZ10, data processing and the
ancillary data. In Section 3 we explain the methods applied to analyze the
data and the equations used to derive the key physical quantities. Finally, in
Section 4 we present our results and discussion, and in Section 5 we summarize
the main conclusions. Throughout this work, we assume a $\Lambda$CDM
cosmology, adopting the values $\Omega_{\Lambda}=0.7$, $\Omega_{\rm M}=0.3$
and H${}_{\rm o}=70$ km s-1 Mpc-1, thus resolving physical scales $\approx$
6.02 kpc per arcsec.

Property | Value
---|---
(1) $z$ | 5.659
(2) $D_{\rm L}$ [Mpc] | 54988
(3) Angular Scale Conversion [kpc/′′] | 6.023
(4) log[SFRIR/(M⊙ yr-1)] | 2.48${}^{+0.15}_{-0.25}$
(5) log[SFRUV/(M⊙ yr-1)] | 1.56${}^{+0.06}_{-0.06}$
(6) log[$M_{\star}$/M⊙] | 10.39${}^{+0.17}_{-0.17}$

Table 1: Compilation of the previous-best global physical quantities of HZ10.
Row (1): Ly$\alpha$ redshift. Rows (2) and (3): Luminosity distance and the
angular scale conversion factor, respectively.Row (4): IR star-formation rate
(SFRIR). Row (5): rest-UV SFR derived from rest-UV luminosity. Row (6):
stellar mass derived from SED fitting to the photometry available in the
COSMOS or GOODS-S fields using LE_PHARE. Rows (1), (5) and (6) are taken from
Capak et al. (2015). Row (4) are drawn from Faisst et al. (2020). Figure 1:
HZ10 HST, ALMA Band 9, 8, and 7, [CII], and Band 6 morphologies in cutouts of
${1.3}^{\prime\prime}\times{2.2}^{\prime\prime}$. Panels A, B, and C contain
the WFC3 F105W, F125W, and F160W images, respectively, in grey scale. The
black contours in the three panels are the [0.9, 1.2, 2.1]$\times 10^{-21}$
erg s-1 cm-2 $\rm\AA^{-1}$ levels of the F105W, F125W, and F160W filters. The
blue-dashed and red-solid contours correspond to [3$\sigma$, 4$\sigma$,
5$\sigma$] and [6$\sigma$, 12$\sigma$, 18$\sigma$] levels for the 77 $\mu$m
continuum and [CII] integrated intensity, respectively. From D to H: Panels
include the 77 $\mu$m ALMA continuum, 110 $\mu$m ALMA continuum, 158 $\mu$m
ALMA continuum, and 158 $\mu$m [CII], and 198 $\mu$m ALMA continuum images,
respectively. Contours levels in panels D and G are the same as in panel A, B,
and C, for 77 $\mu$m continuum, and [CII] integrated intensity, respectively.
While blue-solid contours in panels E and H [3$\sigma$, 4$\sigma$, 5$\sigma$]
levels for 110 $\mu$m and 198 $\mu$m ALMA continuum maps, respectively, the
white-solid contours in panel G are the [6$\sigma$, 4$\sigma$, 5$\sigma$]
levels for the 158 $\mu$m ALMA continuum maps. Finally, while green crosses in
panel B are the centers of the Seérsic profiles computed for the three sources
(e.g., HZ10-C, HZ10-W, and the Bridge) analysed in this work (see §3.2), the
green circles correspond to the apertures used to perform the angular resolved
analysis (see §4.2 for more details).

##  2 Observations

###  2.1 The ALMA-CRISTAL sample

CRISTAL is an ALMA Cycle-8 large program (2021.1.00280.L; PI: R. Herrera-
Camus), aiming to get spatially resolved [CII] line and 870 $\mu$m dust-
continuum emission data for star-forming main-sequence galaxies at $z\sim 4-6$
(Herrera-Camus et al. in prep.).

The CRISTAL sample is drawn from [CII] detected galaxies from the ALPINE
survey (Le Fèvre et al., 2020), all located in the COSMOS or GOODS-S fields.
Galaxies were selected based on spectral energy distribution (SED) modelling
to have: i) specific star-formation rates (sSFR) within a factor of three from
the MS; ii) ancillary HST/WFC3 rest-frame UV data; and iii) stellar masses
log$[M_{\star}/$M${}_{\odot}]\geq 9.5$. Additionally, six extra galaxies from
the COSMOS field (which meet the selection criteria mentioned above) were
added: HZ04, HZ07, HZ10, DC818760, DC873756, VC8326, all with comparable
spatial resolutions and sensitivities to those of the main sample (from ALMA
programs 2018.1.01359.S and 2019.1.01075.S; PI: Manuel Aravena,
2018.1.01605.S; PI: Rodrigo Herrera-Camus, 2019.1.00226.S; PI: Edo Ibar). The
latter also meet the selection criteria mentioned above.

###  2.2 The HZ10 data

HZ10 stands as one of the best CRISTAL sources in order to investigate the
physical conditions of the dust at $z=4-6$ due to the rich multi-wavelength
data available. The latter come mostly from previous studies primarily
oriented to derive the main features of HZ10 using integrated quantities (see
Table 1 and references therein). To analyze HZ10 (CRISTAL-22 in the CRISTAL
survey), we use data from different ALMA cycles/projects listed and described
as follows:

  * •

New Band 9 data: Observations of the rest-frame 77 $\mu$m continuum were taken
on August 19th, 2022, during Cycle 8 ($\nu_{\rm band9}=682$ GHz; panel D in
Fig. 1) as part of project 2022.1.00678.S (P.I.: R. Herrera-Camus). The
integration time was $\sim$77 minutes on-source, achieving a sensitivity of
$\sim 0.24$ mJy beam-1 over $\Delta\nu=13$ GHz and a beamsize of
$\sim{0}^{\prime\prime}.414\times{0}^{\prime\prime}.293$.

  * •

Band 8: Observations of the rest-frame 110 $\mu$m continuum were taken on
January 9th, 2019, during Cycle 6 ($\nu_{\rm band8}=411.4$ GHz; panel E in
Fig. 1) as part of the project 2018.1.00348.S (P.I.: A. Faisst). The
integration time was $\sim 48$ minutes on-source, achieving a sensitivity of
$\sim 0.07$ mJy beam-1 over $\Delta\nu=8$ GHz, and a beamsize of
$\sim{0}^{\prime\prime}.663\times{0}^{\prime\prime}.601$.

  * •

Band 7: Observations of the rest-frame [CII] 158 $\mu$m line emission and the
corresponding continuum were taken between March 26th and March 30th, 2021
during Cycle 7 ($\nu_{\rm band7}=285.3$ GHz; panels F and G, respectively, in
Fig. 1) as part of the project 2019.1.01075.S (P.I.: M. Aravena). The
integration time on source-was $\sim 93$ minutes. While the sensitivity and
the beamsize achieved for the 158 $\mu$m continuum are $\sim 1.6\times
10^{-2}$ mJy beam-1 over $\Delta\nu\approx 6$ GHz, and
$\sim{0.252}^{\prime\prime}\times{0.212}^{\prime\prime}$, respectively, for
the [CII] emission line are $\sim 0.14$ mJy beam-1 at $\Delta\nu\approx 47$
MHz, and $\sim{0}^{\prime\prime}.264\times{0}^{\prime\prime}.225$,
respectively.

  * •

Band 6: Observations of the rest-frame 198 $\mu$m continuum were taken on
January 5th 2015 during Cycle 3 ($\nu_{\rm band6}=220.5$ GHz; panel H in Fig.
1) as part of the project 2015.1.00388.S (P.I.: N. Lu). The integration time
was $\sim 54$ minutes on-source, achieving a sensitivity of $\sim 0.03$ mJy
beam-1 over $\Delta\nu=8$ GHz, and a beamsize of
$\sim{1}^{\prime\prime}.31\times{1}^{\prime\prime}.03$.

HZ10 was included in the CRISTAL sample as part of one of the CRISTAL pilot
programs that resulted as a combination of the ALMA programs 2019.1.01075.S
(P.I.: M. Aravena), and 2012.1.00523.S (P.I.: P. Capak; see Capak et al. 2015
for more details); however, the latter was not considered in this work due to
its coarser angular resolution compared to the former. The details of the data
reduction of the Band 7 observations are presented in Solimano et al. (2024);
among the main features their data processing, the self-calibrated and
combined measurement sets were processed with CRISTAL’s reduction pipeline as
described in Herrera-Camus et al. (in prep.). Briefly, it starts by
subtracting the continuum on the visibility space using using the Common
Astronomy Software Application (CASA; CASA Team et al. 2022) uvcontsub task.
After that, it runs tclean with automasking multiple times, producing cubes
with different weightings and channel widths. In all cases the data are
cleaned down to 1$\sigma$. In this paper, we use datacubes with 20 km s-1
channel width and Briggs (robust=0.5) weighting. The ancillary ALMA Band 8 and
Band 6 observations are presented in Faisst et al. (2020), using Briggs
weighting (Briggs 1995) for the image reconstruction with a robust =0.5.

In addition, we include HST data for HZ10 as part of the project 13641 (P.I.:
Peter Capak) and retrieved from the Barbara A. Mikulski Archive for Space
Telescopes
(MAST111https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html). The
data comprise Wide Field Camera 3 images (WFC3), including F105W (panel A in
Fig. 1), F125W, and F160W bands. These bands cover a rest-frame wavelength
range between $\sim$1200 to 2200$\rm\AA$.

##  3 Methods and products

###  3.1 Basic equations and assumptions

To compute the integrated [CII] line emission and the rest-frame 77 $\mu$m,
110 $\mu$m, 158 $\mu$m, and 198 $\mu$m continuum fluxes, we use the following
equation:

$S_{\rm i}=\int_{A}I_{\rm i}(r)dA,$ (1)

where $A$ is the area of a circular aperture with a diameter equivalent to the
major axis of the 198 $\mu$m continuum beamsize (i.e., the coarsest angular
resolution among the full dataset; $\sim{1}^{\prime\prime}.2$) and centered at
the position of the Bridge (see §4.2 for more details). Here, $I_{\rm i}$ is
the velocity integrated flux density (in Jy km/s beam-1) for [CII], and flux
density (in Jy beam-1) for dust-continuum bands. Finally, $i=$ [CII],
77$\mu$m, $110\mu$m, $158\mu$m, and $198\mu$m.

We also use Equation 1 to compute resolved values for the three components of
HZ10 identified in this work. In this case, $A$ is the area of a circular
aperture with a diameter equal to the major axis of the 77 $\mu$m continuum
beamsize (i.e., the coarser angular resolution between ALMA Band 7 and 9 data;
$\sim{0}^{\prime\prime}.4$) and centered at the locations of the three
components identified in HZ10 (green circles in panels B and D from Fig.1).

We use the SED parametrization from Casey (2012) to compute the best spectral
energy distribution (SED) fitting, which corresponds to:

$S_{\lambda}=N_{\rm bb}\times f(\lambda;\beta_{\rm d},T_{\rm SED}),$ (2)

with

$f(\lambda;\beta_{\rm d},T_{\rm
SED})\equiv\frac{(1-e^{{-(\lambda_{0}/\lambda_{\rm c})}^{\beta_{\rm
d}}})(\frac{c}{\lambda})^{3}}{e^{{(hc)/(\lambda kT_{\rm SED})}}-1},$ (3)

where $N_{\rm bb}$ is the normalization, $\beta_{\rm d}$ is the emissivity
index, $T_{\rm SED}$ is the SED dust temperature, $\lambda_{0}$ is the
wavelength at optical depth $\tau=1$, $\lambda_{\rm c}$ is the power-law
turnover wavelength, and $c$ is the speed of light. We use $\lambda_{0}=100$
$\mu$m for the global and resolved HZ10 SED fittings as adopted in Faisst et
al. (2020). From Casey (2012), we also adopt the parametrization of
$\lambda_{\rm c}=\frac{3}{4}L(\alpha,T)$, where $L(\alpha,T_{\rm
SED})=[(b_{1}+b_{2}\alpha)^{-2}+(b_{3}+b_{4}\alpha)\times T_{\rm SED}]^{-1}$
(with $b_{1}=26.68$, $b_{2}=6.246$, $b_{3}=1.905\times 10^{-4}$, and
$b_{4}=7.243\times 10^{-5}$, and where $\alpha=2.0$; see Casey 2012 for more
details).

We compute the peak dust temperature ($T_{\rm peak}$; e.g. Béthermin et al.
2015; Schreiber et al. 2018), which is derived from the IR emission by the
Wien’s displacement law,

$T_{\rm peak}=\frac{2.898\times 10^{3}({\rm\mu m\,K})}{\lambda_{\rm
peak}({\rm\mu m})}.$ (4)

Since contribution from background CMB heating could potentially affect dust
temperatures at $z>5$ (e.g., da Cunha et al. 2013; Faisst et al. 2020), we
apply CMB corrections to our $T_{\rm SED}$ estimations using Equation (12)
from da Cunha et al. (2013):

$T_{\rm dust}(z)=[(T^{z=0}_{\rm dust})^{4+\beta_{\rm d}}+(T^{z=0}_{\rm
CMB})^{4+\beta_{\rm d}}([1+z]^{4+\beta_{\rm d}}-1)]^{\frac{1}{4+\beta_{\rm
d}}},$ (5)

where $T^{z=0}_{\rm dust}$ and $T^{z=0}_{\rm CMB}$ are the dust temperature
and the CMB temperature measured at $z=0$, respectively. Along this paper, we
adopt $T^{z=0}_{\rm CMB}=2.73$ K .

To compute the total far-IR luminosity ($L_{\rm FIR}$), we perform a numerical
integration of Equation 2 in the wavelength range between 42.5 and 125.5
$\mu$m (e.g., Helou et al., 1988). To do so, we adopt the best-SED fitting
parameters of the source (see §4.2 for more details), and we compute the total
far-IR luminosity, $L_{\rm FIR}$, by integrating the flux between 42.5 and
122.5$\mu$m (as described in Faisst et al. 2020). Similarly, we also compute
the total IR luminosity ($L_{\rm IR}$) by integrating numerically the best-
fitting SED but in the wavelength range between 8 and 1000 $\mu$m.

We obtain the [CII] luminosity (in $L_{\odot}$) using the following equation
(Solomon & Vanden Bout, 2005):

$L_{\rm[CII]}=3.25\times 10^{7}S_{\rm{CII}}\Delta v\,\nu^{-2}_{\rm
obs}\,D^{2}_{\rm L}\,(1+z)^{-3},$ (6)

where $S_{\rm{CII}}\Delta v$ is the velocity integrated flux (in Jy km/s
beam-1), $D^{2}_{\rm L}$ is the luminosity distance (in Mpc), $\nu_{\rm obs}$
is the observed frequency (in GHz), and $z$ is the redshift.

The UV spectral slope ($\beta_{\rm UV}$) is derived from HST WFC3 F125W and
F160W images and using Equation 1 in Liang et al. (2021):

$\beta_{\rm
UV}=\frac{\log(f_{\lambda,0.19})-\log(f_{\lambda,0.23})}{\log(\lambda_{0.19})-\log(\lambda_{0.23})},$
(7)

where $f_{\lambda,0.19}$ and $f_{\lambda,0.23}$ are the specific fluxes at
${\rm 1877\AA}$ and ${\rm 2311\AA}$ rest frame taken from the F125W and F160W
images, respectively. We compute $\beta_{\rm UV}$ by taking advantage of its
almost constant value along the wavelength range ${\rm 1260<\lambda<3200\AA}$,
avoiding contamination by the ${\rm 2175\AA}$ “bump” feature.

Finally, we compute the IR excess (IRX) using the equation (e.g., Meurer et
al. 1999; Popping et al. 2017)

${\rm IRX}=L_{\rm IR}/L_{\rm UV},$ (8)

where $L_{\rm UV}$ is the monochromatic rest-frame UV luminosity at ${\rm
1600\AA}$. In order to derive $L_{\rm 1600\AA}$, we use the HST WFC3/F105W
image (see panel A in Fig. 1) and the PHOTFLAM constant to convert the flux
from e-/s to erg s-1/cm2/$\rm\AA$.

###  3.2 Parametric 2D fitting

As Fig. 1 shows, the ALMA Band 9 and 7 high angular resolution observations
reveal that HZ10 is a complex system with (at least) two main components: a
central one (HZ10-C) and one in the west direction (HZ10-W). To derive the
morphological parameters of these two components, we use the 2D light profile
modeling code PYAUTOGALAXY222https://github.com/Jammy2211/PyAutoGalaxy
(Nightingale et al., 2023), which is based on PYAUTOFIT (Nightingale et al.,
2021). Implementing the image-based filter PYAUTOGALAXY mode for a faster
workflow, we generate the noise-map by feeding the full covariance matrix into
the calculation of the likelihood using the code
ESSENCE333https://github.com/takafumi291/ESSENCE (Tsukui et al., 2023). We
model the HZ10-C and HZ10-W morphologies by fitting a single 2D Sérsic profile
(Sersic, 1968), with a total of seven free parameters: the coordinates of the
source’s center (R.A. and Dec.), the coordinate of the vector of the Sérsic
profile’s minor and major axes $(x,y)$, the effective radius ($R_{\rm e}$),
the Sérsic’s index $(n)$, and the intensity at the center ($I_{0}[r=0]$).

We perform three independent fits to derive the morphologies of the rest-frame
77 $\mu$m continuum (Band 9), and the 158 $\mu$m continuum and [CII] line
emission maps (Band 7). While the continuum maps are directly generated after
running the tclean task in CASA, we produce the [CII] intensity map (or moment
0 maps) by fitting a Gaussian function to the line profile and integrating the
emission in the spectral range [$\mu$-FWTM,$\mu$+FWTM] (where $\mu$ is the
central frequency and FWTM is the full width at one tenth of maximum of the
Gaussian profile). We use the following method to compute the best
morphological parameters of the emission:

  1. 1.

We estimate the centroid using 2D Gaussian functions (R.A. and Dec.) of the 77
$\mu$m continuum emission of HZ10-C and HZ10-W. We use these coordinates as an
initial guess to look for the centers of the two 2D Sérsic profiles.

  2. 2.

Then, we use the task minimize from the PYTHON package lmfit (Newville et al.,
2015). We perform a two step best-parameters searching: we adopt the
least_squares method, followed by the emcee method (the latter looks for the
maximum likelihood via a Monte-Carlo Markov Chain). To do so, we constrain the
values of $n$ and $I_{0}$ within the ranges $[0.1,3]$ and $[0,2I_{\rm max}]$,
respectively ($I_{\rm max}$ is the maximum value of the intensity map).

  3. 3.

Finally, we repeat step 2 for the 158 $\mu$m continuum and the [CII] line
emission maps, using this time the 77 $\mu$m continuum best parameters as a
prior (i.e., centroids, Sérsic indexes, $I_{\rm max}$).

Table 2 lists the best parameters and uncertainties of the fitting procedure.
We note that the differences in DEC in Table 2 are larger than the spatial
resolutions of the observations in each band (see §2.2).

##  4 Results and Discussion

Parameter |  Value
---|---
|  R.A.:10h00m, Dec.: +01${}^{\circ}{33}^{\prime}$

|  |  | HZ10-C |  |  |  | HZ10-W
---|---|---|---|---|---|---|---

77 $\mu$m continuum

Center (R.A.) | 59.304${}^{\rm s}\pm$0.037s | 59.255${}^{\rm s}\pm$0.091s
---|---|---
Center (Dec.) |  ${19.507}^{\prime\prime}$$\pm$${0.038}^{\prime\prime}$ |  ${19.408}^{\prime\prime}$$\pm$${0.037}^{\prime\prime}$
Effective radius $R_{\rm e}$ [kpc] | 1.26$\pm$0.26 | 1.43$\pm$0.34
Sérsic index ($n_{\rm 77\mu m}$) | 1.70$\pm$0.42 | 1.84$\pm$0.38
Axis ratio (min/maj) | 0.66$\pm$0.10 | 0.60$\pm$0.08
$I_{0}$(r=0) [mJy/bm] | 0.61$\pm$0.25 | 0.22$\pm$0.24

158 $\mu$m continuum

Center (R.A.) | 59.297s$\pm$$0.039^{\rm s}$ | 59.254s$\pm$$0.040^{\rm s}$
---|---|---
Center (Dec.) |  ${19.486}^{\prime\prime}$$\pm$${0.037}^{\prime\prime}$ |  ${19.396}^{\prime\prime}$$\pm$${0.038}^{\prime\prime}$
Effective radius $R_{\rm e}$ [kpc] | 0.78$\pm$0.25 | 1.00$\pm$0.23
Sérsic index ($n_{\rm 158\mu m}$) | 1.63$\pm$0.49 | 1.54$\pm$0.19
Axis ratio (min/maj) | 0.58$\pm$0.07 | 0.59$\pm$0.04
$I_{0}$(r=0) [mJy/bm] | 0.24$\pm$0.04 | 0.12$\pm$0.03

[CII] line emission

Center (R.A.) | 59.298${}^{\rm s}\pm$0.039s | 59.253s$\pm$$0.048^{\rm s}$
---|---|---
Center (Dec.) |  ${19.492}^{\prime\prime}$$\pm$${0.040}^{\prime\prime}$ |  ${19.405}^{\prime\prime}$$\pm$${0.039}^{\prime\prime}$
Effective radius $R_{\rm e}$ [kpc] | 1.49$\pm$0.23 | 1.03$\pm$0.10
Sérsic index ($n_{\rm[CII]}$) | 2.76$\pm$0.12 | 2.21$\pm$0.31
Axis ratio (min/maj) | 0.64$\pm$0.02 | 0.48$\pm$0.07
$I_{0}$(r=0) [Jy/bm km/s] | 1.18$\pm$0.09 | 1.30$\pm$0.27

Table 2: Results of the parametric 2D fitting of Sérsic profiles, using
PYAUTOGALAXY, of the 77 $\mu$m continuum (top), 158 $\mu$m continuum (middle),
and [CII] 158 $\mu$m line emission (bottom) for HZ10-C and HZ10-W. Figure 2:
Results from the parametric 2D modeling of HZ10 using PYAUTOGALAXY, in panels
of ${3}^{\prime\prime}\times{3}^{\prime\prime}$. While left panels contain the
observed emission (colormap and black-solid contours), middle and right panels
show the best model using 2D Sérsic profiles (see Table 2 for more details)
and the residuals after subtracting the maximum likelihood model,
respectively. From top to bottom, panels contain the rest-frame 77 $\mu$m ALMA
continuum (top row), 158 $\mu$m ALMA continuum (middle row), and [CII] 158
$\mu$m line emission (top row). The beamsize of the data are represented by
ellipses at the bottom-left of middle panels. For all the subplots, the
contour are the [0 (dashed), 4$\sigma$, 5$\sigma$, and 6$\sigma$] levels. The
figure confirms the binary nature of HZ10, which can be decomposed in HZ10-C
(at the center), HZ10-W (to the left), and the Bridge in between the two main
components (with at least a 5$\sigma$ significance on the residual 158 $\mu$m
continuum map). The latter seems to reflect the extended dusty component
connecting HZ10-C and HZ10-W.

###  4.1 The structure of HZ10

Figure 2 shows the results of the 2D parametric fitting, confirming the
complex structure of HZ10. The figure shows this system has (at least) two
main components, HZ10-C and HZ10-W, each of them slightly showing different
morphologies depending on the datasets. On the one hand, HZ10-W shows a
slightly more extended distribution of the [CII] emission than HZ10-C; on the
other hand, the two components have surprisingly similar Sérsic profiles when
analyzing the 77 and 158 $\mu$m continuum maps. These results may reflect
significant dust content (relative to the gas) in HZ10-W that causes severe
dust attenuation, as evidenced by its faint UV emission (shown in panel A of
Fig. 1).

Figure 3: Best-fitting IR SEDs for the global ALMA Bands 9, 8, 7, and 6
continuum emission of HZ10 (77 $\mu$m, 110 $\mu$m, 158 $\mu$m, and 198 $\mu$m,
respectively). The wavelength is given in rest frame. The blue-solid line
corresponds to the best-modified SED given by equations 2-4 and using the
maximum likelihood parameters (see top-left corner). The combination between
new ALMA Band 9 and 7 continuum measurements (77 $\mu$m and 158 $\mu$m rest
frame, respectively) allow us to probe the peak of the dust SED, and hence
constrain these parameters accurately. Errorbars are smaller than symbols.

Figure 4: Best-fittings for far-IR SEDs for 77 and 158 $\mu$m continuum
emission, using equations 2-3 and adopting the $\beta_{\rm d}$ obtained in
Fig. 3, for the three components identified in the analysis from §4.1: HZ10-C
(left panel), HZ10-W (middle panel), and the Bridge (right panel). Conventions
are as in Fig. 3. As shown in Figs. 1 and 2, HZ10-C and the Bridge have lower
77 $\mu$m continuum emission compared to that for HZ10-W (although still
within the uncertainties); this may be reflecting a strong UV dust absorption
in HZ10-W, which translates into a higher $T_{\rm SED}$ temperature when
compared to those of the other two components (although all of them with the
same temperature within 1$\sigma$).

Interestingly, the middle row of Figure 2 show that the HZ10 system cannot be
described just considering two Sérsic profiles. A visual inspection of the
residual map of the 158 $\mu$m continuum emission reveals a third component
connecting HZ10-C with HZ10-W. Applying the procedure described in §3.2 to the
residual 158 $\mu$m map, we note that the extra component (or the Bridge)
cannot be well described by a Sérsic profile (after trying to fit a Sérsic
profile to the residuals). To analyze the physical properties of the dust at
the Bridge, we derive the coordinates of the centroid of its 158 $\mu$m
continuum emission (R.A.Int, Dec.Int); we obtain (10h00m59s.279$\pm$$0^{\rm
s}$.037,
+01${}^{\circ}{33}^{\prime}{19}^{\prime\prime}$.463$\pm$${0}^{\prime\prime}$.038).

The Bridge may suggests either the interaction or a tidal tale connecting the
two main components of the HZ10 system. We are conducting a more detailed
analysis of the [CII] morphology to further probe the interacting nature of
the system to be published in a forthcoming paper (Telikova et al. in prep).
This third component is also presented in new and deep JWST/NIRSpec IFU
observations of the main nebular lines in HZ10 (Jones et al., 2024).

###  4.2 SED fitting, $T_{\rm SED}$ and $T_{\rm peak}$

Source | $S_{\rm 77\mu m}$ | $S_{\rm 158\mu m}$ | $T_{\rm SED}$ | $T_{\rm peak}$ | log[$L_{\rm FIR}$] | log[$L_{\rm IR}$] | log[$L_{\rm[CII]}$] | log[$\Sigma_{\rm FIR}$] | $\beta_{\rm UV}$ | IRX
---|---|---|---|---|---|---|---|---|---|---
| (mJy) | (mJy) | (K) | (K) | ($L_{\odot}$) | ($L_{\odot}$) | ($L_{\odot}$) | ($L_{\odot}$ kpc-2) |  |
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
Global | 3.78$\pm$0.63 | 1.76$\pm$0.51 | 46.7$\pm$6.8 | 35.5$\pm$7.0 | 12.42$\pm$0.12 | 12.60$\pm$0.27 | 10.48$\pm$0.14 | 9.59$\pm$0.14 | -1.59$\pm$0.24 | 2.27$\pm$0.98
HZ10-C | 1.04$\pm$0.38 | 0.42$\pm$0.20 | 46.4$\pm$11.7 | 35.3$\pm$11.9 | 11.83$\pm$0.38 | 12.01$\pm$0.32 | 9.89$\pm$0.15 | 9.94$\pm$0.15 | -1.63$\pm$0.03 | 1.71$\pm$1.01
HZ10-W | 1.06$\pm$0.39 | 0.36$\pm$0.13 | 51.2$\pm$13.1 | 37.8$\pm$13.3 | 11.89$\pm$0.45 | 12.06$\pm$0.33 | 9.65$\pm$0.20 | 9.98$\pm$0.20 | -1.80$\pm$0.14 | 2.71$\pm$1.32
Bridge | 0.79$\pm 0.37$ | 0.44$\pm 0.17$ | 39.2$\pm$9.8 | 35.9$\pm$10.0 | 11.71$\pm$0.33 | 11.88$\pm$0.37 | 9.76$\pm$0.42 | 9.80$\pm$0.42 | -1.74$\pm$0.15 | 2.38$\pm$1.59

Table 3: Main properties of the HZ10 system derived in this work. Column (1):
component name. Column (2) and (3): 77 and 158 $\mu$m continuum emission flux,
respectively. Column (4) and (5): logarithmic of the SED and peak
temperatures, respectively. Column (6): logarithmic of the total FIR
luminosity computed in the range of 42.5 and 125.5 $\mu$m. Column (7):
logarithmic of the total IR luminosity computed in the range of 8 and 1000
$\mu$m. Column (8): total luminosity of the [CII] emission line computed by
using Equation 6. Column (9): FIR luminosity surface density computed after
dividing column (7) by the area of a circular region with a diameter
equivalent to the major axis of the ALMA 77 $\mu$m continuum beamsize (i.e.,
$\sim 2.5$ kpc). Column (10): UV spectral slope. Column (11): IR excess.

We perform far-IR SED fitting in order to determine the physical conditions of
the dust in HZ10. We use Equation 1 to compute the global values of the 77
$\mu$m, 110 $\mu$m, 158 $\mu$m, and 198 $\mu$m dust-continuum fluxes, which
correspond to the sum of all the flux within a circular area with a diameter
equivalent to the major axis of the 198 $\mu$m continuum beamsize. We recall
that the latter is the coarsest angular resolution among all the dataset
(i.e., $={1.}^{\prime\prime}2\approx 7$ kpc at HZ10’s distance). The aperture
is placed at the 158 $\mu$m continuum emission centroid of the Bridge (see
§4.1). To do this, we convolve all the maps to the coarsest angular resolution
among the continuum maps (i.e., 198 $\mu$m dust-continuum beamsize). We then
use these fluxes to look for the best-fit parameters of equations 2 and 3 by
performing the task minimize from the PYTHON package lmfit (Newville et al.,
2016). We adopt the emcee method, which determines for the maximum likelihood
of the parameters via a Monte-Carlo Markov Chain. We finally choose fixed
values for $\alpha=2.0$ and $\lambda_{\rm 0}=100$ $\mu$m to perform the
fitting (see §3.1). The results, including the best parameters for the far-IR
SED fitting and the $1\sigma$ uncertainties curves, are shown in Figure 3. We
highlight that ALMA Band 9 data, which is closer to the SED peak of HZ10 (see
Fig. 4), allow us to reduce the uncertainties by almost three times of $T_{\rm
SED}$ estimations when compared to global estimations included in previous
studies (e.g., Faisst et al. 2020; Mitsuhashi et al. 2023b).

We also perform the far-IR SED fitting for the three sources described in
§4.1. Similarly to way than for the global SED fitting described above, we use
a circular aperture with a diameter equivalent to the major axis of the 77
$\mu$m continuum (i.e., the coarsest angular resolution among the [CII] line,
77, and 158 $\mu$m continuum emission maps; $\sim{0}^{\prime\prime}.4$, or
$\sim 2.5$ kpc) to compute the integrated fluxes within such aperture
(included in Table 3). The aperture is located at three different positions
given by the centers of the 77 $\mu$m continuum Sérsic profiles for HZ10-C and
HZ10-W (see Table 2), and the 158 $\mu$m continuum emission centroid of the
Bridge. The integrated fluxes are shown in Table 3.

When comparing our results from the integrated SED fitting with previous
studies, we find a dust emissivity index $\beta_{\rm d}=2.00\pm 0.14$ slightly
lower (although still consistent) than that derived from Faisst et al. (2020)
($\beta_{\rm d}=2.15^{+0.41}_{-0.54}$; based on the SED modelling of Band 6, 7
and 8 data). Conversely, we obtain a consistent global SED dust temperature
when compared that derived by Faisst et al. (2020) to ours
($46.2^{+16.2}_{-8.5}$ and $46.7\pm 6.8$ K, respectively). We remark that
$\beta_{\rm d}$, $T_{\rm SED}$, and $T_{\rm peak}$ (derived from Equation 4)
are very sensitive to both the assumptions for the fixed values and the
completeness of the dust-continuum dataset.

For the resolved SED far-IR SED fitting, and considering the caveat above, we
adopt a fixed value for $\beta_{\rm d}=2.0$ (i.e., the value derived from the
global SED fitting). We find slightly dissimilar $T_{\rm dust}$ estimations
when compared to the integrated quantities included in Faisst et al. (2020).
Our results of the far-IR SED fitting for the three components, including the
best parameters and the $1\sigma$ uncertainty curves, are shown in Figure 4.
Although still within the uncertainties, we note that the resolved structure
of the dust revealed by the data presented in this work shows that the $T_{\rm
SED}$ depends critically on the component of the HZ10 system. While the SED
dust temperature in HZ10-C is consistent with integrated values presented by
previous studies (46.4 K), we note that HZ10-W’s $T_{\rm SED}$ is about 5 K
higher than that found by Faisst et al. (2020) (51.2 and 46.2 K,
respectively). Interestingly, we also find that the Bridge has a lower $T_{\rm
SED}$ values (although still consistent; 39.2$\pm$9.8) compared to HZ10-C and
HZ10-W. The latter seems to be reflecting the detached nature of the Bridge
respect to the two main components (perhaps due to outflows) and without close
sources of hard radiation fields required to increase its $T_{\rm dust}$ up to
temperatures similar to those of HZ10-C and HZ10-W.

Figure 5: Peak dust temperature ($T_{\rm peak}$) evolution with redshift for
HZ10-C, HZ10-W, and the Bridge as blue, orange, and green squares,
respectively. The figure also includes galaxy samples at different redshift
ranges: $z=0.2-4$ (gray unfilled circles from ALESS; da Cunha et al. 2015;
yellow triangles, Béthermin et al. 2015), and $z>6$ (purple circles; Knudsen
et al. 2016; Hashimoto et al. 2019; turquoise pentagons, Algera et al. 2024).
The red triangle, black circle, and brown circle are galaxies at $z\sim 5$
included in Faisst et al. (2020). Magenta stars correspond to $z\sim 7$
galaxies selected from the Reionization Era Bright Emission Line Survey,
REBELS (Sommovigo et al., 2022). Blue large dark-blue circles correspond to
galaxies from at $0.5<z<4$ selected from the deep CANDELS fields, and the
dashed-black line correspond to their linear best-fit, both extracted from
Schreiber et al. (2018). Green inverted triangles correspond to results based
on the stacking analysis and the expected redshift evolution from Viero et al.
(2022). The blue dot-dashed line is the expected peak temperature evolution
from hydrodynamic simulation (Liang et al. 2019; Ma et al. 2019). Red unfilled
circles correspond to galaxies at $z\gtrsim 6$ selected from the Systematic
Exploration in the Reionization Epoch using Nebular And Dust Emission
(SERENADE; Harikane et al. in prep.), as included in Mitsuhashi et al.
(2023a). The sample also contains three galaxies selected from Harikane et al.
(2020). Finally, red and green dashed areas are correspond to the parameter
space covered by galaxies selected at different redshfits from the Herschel
Multi-tiered Extragalactic Survey (HerMES; Magnelli et al. 2014) and
ultraluminous infrared galaxies selected from the Great Observatories All-sky
LIRG Survey (GOALS; U et al. 2012), respectively.

To put the results from HZ10 in a more general context, Figure 5 shows the
evolution of peak dust temperature in star-forming galaxies (expected
evolution of the peak temperature derived from Eq. 4, $T_{\rm peak}$) as a
function of redshift. The combination of $T_{\rm peak}$ measurements for
low-$z$ galaxies and the interpolation of these estimations at high redshift
suggests that temperatures increase up to $z\sim 4$ and then flatten off and
then flatten off following hydrodynamical simulations from Liang et al. (2019)
(blue dashed-dotted line in Fig. 5). The scatter from the linear relation
expected for galaxies at $0.5<z<4$ (see black-dashed line in Fig. 5) has been
shown to depend on many factors. Some of them have been identified in galaxies
from the local Universe (e.g., KINGFISH; Skibba et al. 2011, GOALS; U et al.
2012), including changes in the dust mass density, effects on the dust
opacity, and/or variability in the UV luminosity of a central source (e.g.,
young stars or AGN activity). For instance, Faisst et al. (2017) propose that
metallicity also has a significant effect on the $T_{\rm peak}$ measured. In
particular, their results indicate that galaxies with faint IR emission and
low metallicity could have $T_{\rm peak}$ values similar to those for high IR
luminous galaxies. In addition, metallicity can potentially alter the dust
properties (e.g., Pak et al. 1998; Misselt et al. 1999; Sommovigo et al.
2022), mainly produce environments with lower opacities (e.g., Issa et al.
1990; Lisenfeld & Ferrara 1998) due to harder stellar radiation fields are
expected for low metallicity environments altering the dust temperature.

Based on the HST WFC3/F160W data of HZ10-W shown in panel A of Figure 1, the
higher $T_{\rm SED}$ may be responding to UV emission severely attenuated by
dust, producing a significant increase of the dust temperature. Spectroscopic
observations of the main nebular lines with JWST/NIRSpec (e.g., Jones et al.
2024) and future ALMA Band 10 continuum and JWST observations could help us to
get better constraints on the dust emission and metallicity, allowing us to
break down the potential degeneracy of the $T_{\rm SED}$ estimations.

###  4.3 The Physical Properties of the ISM in HZ10

Figure 6 presents the relation between far-IR luminosity, $L_{\rm FIR}$, and
$T_{\rm SED}$ for the three components of HZ10 identified in this work. The
values of $L_{\rm FIR}$ for HZ10-C, HZ10-W, and the Bridge are consistent with
the global far-IR luminosity derived by Faisst et al. (2020). However,
compared to the expected relation from the optically thick case (i.e., where
$L\propto T^{4}$; grey-dashed line in Fig. 6), the HZ10 components seem to be
either too warm (i.e., to the right of the relation) or underluminous in the
far-IR (thus, below the relation). In addition, when compared to dusty star-
forming galaxies at similar redshifts (DSFG, included as gray solid circles in
Fig. 6; Riechers et al. 2020), our $L_{\rm FIR}$ values may be reflecting
intrinsic differences of the dust properties between some CRISTAL and DSFG
systems, but most likely due different spatial configuration and/or optically
thickness. Some of these differences could respond as well to the factors
described at the end of §4.2, including the effects of metallicity, dust
abundance on the dust opacity (Faisst et al., 2017), or changes in the
photoelectric efficiency of the dust (e.g., Nath et al. 1999; Malhotra et al.
2017; McKinney et al. 2021; Glatzle et al. 2022).

The [CII]-to-FIR luminosity ratio, $\mbox{[C{\small
II}]}/$FIR=$L_{\rm[CII]}/L_{\rm FIR}$, has been found to be closely related to
the physical properties of the ISM. For instance, while in low-metallicity
environments the typical values are log($\mbox{[C{\small II}]}/$FIR)$\sim-2$
(e.g., dwarf galaxies or star-forming disks; Madden et al. 2013; Smith et al.
2017; Herrera-Camus et al. 2018a), in nuclear regions, starburst systems, and
AGNs are around log($\mbox{[C{\small II}]}/$FIR)$\sim$[-4,-3] (e.g., Malhotra
et al. 2001; Díaz-Santos et al. 2013; Herrera-Camus et al. 2018a).

Figure 6: Comparison of far-IR luminosity ($L_{\rm FIR}$) and SED temperatures
($T_{\rm SED}$) for the three components of HZ10 covered in this work. The
figure includes a sample of dust star-forming galaxies at $z>5$ (gray solid
circles; Riechers et al. 2020), strongly gravitational lensed dusty star-
forming galaxies at $1.9<z<6.9$ included in Reuter et al. (2020), and previous
comparisons for HZ09 (black unfilled circle), HZ04 (red empty triangle), and
HZ10 (purple filled circle), as included in Faisst et al. (2020). The gray
dashed line is the $L\propto T^{4}$ relation, also extracted from Faisst et
al. (2020). Figure 7: [CII]/FIR ratio as a function of the FIR surface density
($\Sigma_{\rm FIR}$) for HZ10-C (blue square), HZ10-W (orange square), and the
Bridge between them (green empty square). For comparison, the figure also
includes nearby star-forming and starbursts galaxies from the SHINING sample
(grey filled stars; Herrera-Camus et al. 2018a), $\sim 100$ pc scale regions
in the central disk of M82 (green unfilled dots; Contursi et al. 2013;
Herrera-Camus et al. 2018a), $\sim 400$ pc scale regions from central regions
of M83 (purple filled circles), lensed dusty star-forming galaxies at $z\sim
1.9-5.7$ (DSFGs; red-solid diamonds; Spilker et al. 2016) , and four kpc-size
regions extracted across the disk of HZ04 (red unfilled triangles; Herrera-
Camus et al. 2021). The The solid grey line corresponds to the best quadratic
fit to the SHINING data (as included in Herrera-Camus et al. 2018a).

Figure 7 shows $\mbox{[C{\small II}]}/$FIR as a function of the far-IR
luminosity surface density, $\Sigma_{\rm FIR}$, for the three components of
HZ10 analyzed in this work. We note that HZ10-C, HZ10-W, and the Bridge have a
very smooth distribution of [CII]-to-FIR values, with a variation at most as
$25$% between the faintest (i.e., the Bridge) and the strongest continuum
emission (i.e., HZ10-W). Although differences in the [CII] emission are
slightly more significant, these are at 35-45% between the faintest (i.e., the
Bridge and HZ10-W) and the strongest [CII] emitter (HZ10-C). Our results are
consistent with those of $\sim 100$ pc scale regions in the central disk of
the nearby starburst galaxy M82 (green unfilled dots; Contursi et al. 2013;
Herrera-Camus et al. 2018a). Following a similar methodology to that described
in §4.2 to the physical parameters of the dust in a typical star-forming
galaxy at $z\sim 5.5$ ($\beta_{\rm d}=1.5$, $T_{\rm SED}=45$ K; e.g., Pavesi
et al. 2016; Faisst et al. 2017), Herrera-Camus et al. (2021) compute the
$\mbox{[C{\small II}]}/$FIR-$\Sigma_{\rm FIR}$ relation for HZ04 (CRISTAL-20
in the CRISTAL survey). Splitting the analysis in four independent beams (red
unfilled triangles in Fig. 7), they find log($\mbox{[C{\small
II}]}/$FIR)$\sim$$[-2.7,-2.3]$ and log($\Sigma_{\rm FIR}$/[L⊙
kpc-2])$\sim$10\. Although $\mbox{[C{\small II}]}/$FIR values for our sources
are comparable to those for HZ04, we note that our sources have $\Sigma_{\rm
FIR}$ values (directly the star-formation rate, SFR) around 4 times higher
than those derived in Herrera-Camus et al. (2021). These results may suggest a
deficit in the [CII] content in the HZ10 system, specially in HZ10-W, although
we note that the dust SED setup used in Herrera-Camus et al. (2021) is
different to ours and therefore this may affect the comparison between the two
sources given its impact on the derived IR luminosities.

Several studies have shown the strong impact that AGNs could potentially have
on the [CII] line emission, producing low [CII]-to-FIR luminosity ratios in
nearby galaxies (e.g., Stacey et al. 2010; Sargsyan et al. 2012; Herrera-Camus
et al. 2018b). Díaz-Santos et al. (2013) suggest that in galaxies with a
strong [CII] deficit (log($\mbox{[C{\small II}]}/$FIR)$<-3$) the AGNs can play
an important role in destroying significant fractions of polycyclic aromatic
hydrocarbon (PAH) molecules (i.e., associated with [CII] emission) due to
their hard radiation field (e.g., Lai et al. 2023). This is consistent with
the evidence (although marginal) from the new JWST/NIRSpec observations that
suggest that HZ10-W may contain nuclear activity (Jones et al., 2024). On the
other hand, by modeling the emission of the molecular, neutral, and ionized
gas of galaxies selected from the Survey with Herschel of the Interstellar
medium in INfrared Galaxies (SHINING), Graciá-Carpio et al. (2011) found that
a decrease in $\mbox{[C{\small II}]}/$FIR ratios can be explained by
increasing the value of the ionization parameter ($U$) on the surface of
molecular clouds. As $U$ increases, a larger fraction of UV photons are
absorbed by dust in the ionized region and reemitted in the form of infrared
emission. The net effect is that the fraction of UV photons available to
ionize and excite the gas is reduced at high $U$, decreasing the relative
intensity of the fine structure lines compared to the FIR continuum (e.g,
Voit, 1992; Abel et al., 2009; Graciá-Carpio et al., 2011; Herrera-Camus et
al., 2018a). However, the UV photons could be also interacting with gas
itself, having a potential impact on its physical properties (i.e., the [CII]
line emission).

###  4.4 The IRX-$\beta_{\rm UV}$ relation

Unbiased galaxy star-formation rates at low and high redshift are essential to
get a complete picture of the mechanisms changing the physical properties of
the ISM as a function of cosmic times. It is critical thus to account for both
the dust thermal emission and the UV light to mitigate the natural biases on
the derivation of accurate SFR estimations. In this sense, the IRX-$\beta_{\rm
UV}$ relation (e.g., Meurer et al. 1995, 1999) has shown to be a strong
observational tool to link the IR emission to UV measurements, particularly to
the tight correlations of the IRX-$\beta_{\rm UV}$ plane exhibited at
different redshift ranges (e.g., Heinis et al. 2013; McLure et al. 2018;
Fudamoto et al. 2020; Bowler et al. 2024). In the last decades, however,
several studies have revealed that some local LIRGs/ULIRGs (e.g., Howell et
al. 2010) and high-$z$ galaxies (e.g., Álvarez-Márquez et al. 2016; Bouwens et
al. 2016; Reddy et al. 2018) can depart from this relation due to changes in
intrinsic dust properties. Such variations can be produced by the composition
of the dust, the spatial distribution of the dust and UV emission, or due to
ISM turbulence, among others (see Liang et al. 2021 and references therein).

To test which is the more adequate scenario in the HZ10 system, we use
equations 7 and 8 to drawn the IRX-$\beta_{\rm UV}$ relation for HZ10-C,
HZ10-W, and the Bridge, as shown in Figure 8. In addition to several galaxy
samples at different redshifts, we have also overplotted some of the
schematics from Figure 11 in Popping et al. (2017), which reflect the distinct
physical mechanisms affecting the properties of the emission sources (see
light-blue and light-orange shaded areas, and the green arrow at the bottom of
Fig. 8). While HZ10-C is close to the best linear-fit for UV-selected near
starburst galaxies (black solid line; Meurer et al. 1999), we note that both
HZ10-W and the Bridge lie significantly above it. According to Popping et al.
(2017) models, the location of HZ10-W and the Bridge in the IRX-$\beta_{\rm
UV}$ relation seems to favor scenarios where a screen of dust is placed with
holes in between a relative young stellar population and the observer.
Although a small level of turbulence within the dust could also increase its
optical depth (e.g. Fischera et al. 2003; Popping et al. 2017), the dust
screen appears to be a simple (yet reasonable) explanation for the dust
emission in the more dusty components of the HZ10 system. For example,
although with higher $\beta_{\rm UV}$ than those for the three components
analyzed in this work, the star-forming main-sequence galaxies analyzed by
Fudamoto et al. (2020) have similar IRX values than HZ10’s. They propose that
the high dust-attenuation properties shown by those galaxies may correspond to
supernovae (SNe) driven dust production at $z\geq 2-3$. Such SNe dust could be
consistent with the steeper dust curve observed at $z\sim 4-6$ galaxy sample
compared to the attenuation inferred at $z<3$ for sources similar to those
included in Meurer et al. (1999)(e.g., Maiolino et al. 2004; Hirashita et al.
2005; Gallerani et al. 2010).

Numerical simulations have also found that [CII]/FIR ratios can decrease with
increasing $\Sigma_{\rm FIR}$, leading to an apparent [CII] deficit. In
particular, Bisbas et al. (2022) show that this could reflect the thermal
saturation of [CII] as a consequence of the strong far-UV heating related to
the high SFR. They propose that while the [CII] emissivity increase
asymptotically in this regime, the FIR emission increases linearly, leading
the deficit of the [CII].

The IRX-$\beta_{\rm UV}$ relation has been shown to be a powerful tool to
characterize the complex structure of the HZ10 system. However, some studies
have suggested the that IRX-$\beta_{\rm UV}$ relation may fail in probing the
physical conditions of the ISM in galaxies with large IR-to-UV flux ratios
compared to the $\beta_{\rm UV}$ (e.g., Ferrara et al. 2022). To address this
problem, Ferrara et al. (2022) introduce the non-dimensional “molecular
index”, $I_{\rm m}=(F_{158\mu m}/F_{\rm 1500\AA})/(\beta_{\rm UV}-\beta_{\rm
int})$; here, $F_{158\mu m}$ and $F_{\rm 1500\AA}$ are the fluxes at $1500\AA$
and the observed far-infrared continuum flux at 158 $\mu$m, respectively, and
$\beta_{\rm int}$ is the intrinsic UV spectral slope (typically, $\beta_{\rm
int}\approx-2.406$; see Ferrara et al. 2022 for more details). When applying
to sources selected from REBELS, they note that $I_{\rm m}$ is a good
predictor of the simultaneous presence of optically thin and thick regimes,
showing that galaxies with a two-phase medium have $I_{\rm m}>1120$. When
computing the molecular index for the HZ10 system, we obtain $I_{\rm
m}=1156\pm 104$. The latter is consistent with the IRX-$\beta_{\rm UV}$
relation for HZ10 discussed above, supporting the scenario of a two-phase
structure with young stars embedded in optically thick giant molecular clouds
and an older stellar populations immersed in a more transparent medium.

Future ALMA-CRISTAL studies (Killi et al. in prep.) will analyze the
variations of the IRX-$\beta_{\rm UV}$ relation among CRISTAL galaxies on
$\sim$kiloparsec scales, therefore allowing us to derive more statistically
significant results and to perform a meticulous comparison with other similar
studies from the literature.

Figure 8: Comparison between the IR excess (IRX) and the UV spectral slope
($\beta_{\rm UV}$); i.e., the IRX-$\beta_{\rm UV}$ relation for HZ10-C,
HZ10-W, and the Bridge between them. Conventions are as in Fig. 7. Green
diamonds correspond to the MASSIVEFIRE sample at $z=6$, as included in Liang
et al. 2021. The figure also encompasses a series of galaxy samples from the
literature: yellow triangles are taken from Heinis et al. (2013), red stars
from Álvarez-Márquez et al. (2016), orange stars from Bouwens et al. (2016),
blue-edge diamonds from Reddy et al. (2018), magenta-edged circles from McLure
et al. (2018), and cyan-edged squares from Fudamoto et al. (2020). The black
solid line is the best linear-fit for the IRX-$\beta_{\rm UV}$ relation for
UV-selected starburst galaxies as shown in Meurer et al. (1999). The light-
blue/light-orange shaded areas, and the green arrow at the bottom are
extracted from the schematic figure included in Popping et al. (2017), which
summarizes the different physical mechanisms affecting the properties of the
emission sources.

##  5 Summary and conclusions

We present a study of the dusty, main-sequence galaxy HZ10 at $z\approx 5.66$,
as part of the CRISTAL survey. We present new ALMA [CII] line emission and
Band 9 dust continuum data, which is key to constrain the peak of the dust
SED. In combination with ALMA Band 6, 7, 8 and HST WFC3 archival data, we
conduct a systematic analysis of the morphology and the physical conditions of
dust and gas in such multi-component system. We characterize the main
properties of the structures identified in HZ10, such as the SED and peak dust
temperatures, FIR luminosities, [CII]/FIR ratios, UV spectral slopes, and
their IR excesses, and we compare our results with the current literature. Our
main conclusions are enumerated as follows:

  1. 1.

We perform a parametric 2D fitting to derive the morphological parameters of
the two main components of HZ10 (HZ10-C and HZ10-W), which are well described
by Sérsic profiles in the 77 and 158 $\mu$m dust continuum and the [CII] line
emission data (see Table 2). Interestingly, the residual map of the 158 $\mu$m
dust continuum data reveals a third component (Bridge), which seems to
correspond to a bridge-like dusty structure connecting the central and west
components.

  2. 2.

We carry out a global modified blackbody SED fitting of HZ10 using the 77,
110, 158, and 198 $\mu$m dust continuum fluxes. We derive a global dust
emissivity index $\beta_{\rm d}\approx 2.0\pm 0.14$ and a global dust
temperature of $T_{\rm SED}=46.7\pm 6.8$ K. We adopt the global value of
$\beta_{\rm d}$ to perform modified SED fittings of the three main structures
identified in the HZ10 system using the 77 and 158 $\mu$m dust continuum maps
(since these two datasets allow us to resolve spatially the three components
of HZ10), obtaining spatially-resolved estimations of their SED and peak dust
temperatures. We find that HZ10-W, which is the component that shows the
higher obscuration of the rest-frame UV emission, has a dust temperature of
$T_{\rm SED}=51.2\pm 13.1$ K, which is about $\sim 5$ K higher than the two
other components (although all of them with the almost the same SED
temperature within 1$\sigma$). More importantly, the inclusion of the new ALMA
Band 9 continuum data allow us to reduce the uncertainties in the global dust
temperature measurements by a factor of $\sim 2.3$.

  3. 3.

We compute the [CII]-to-FIR luminosity ratio, [CII]/FIR, for HZ10-C, HZ10-W,
and the Bridge. When we compare [CII]/FIR with the FIR surface density,
$\Sigma_{\rm FIR}$, we find that our sources cover a similar parameter space
to that of local starburst galaxies. We note that HZ10-W shows signs of a
[CII] deficit, suggesting possibles scenarios such as hard radiation field
destroying PAHs associated with [CII] emission (e.g., young stellar
populations or AGN activity), or variations in the dust photoelectric
efficiency.

  4. 4.

We calculate the IR excesses and the UV spectral slopes (the IRX-$\beta_{\rm
UV}$ relation), for the three components in HZ10. While HZ10-C has IRX and
$\beta_{\rm UV}$ values consistent with those of UV-selected starburst local
galaxies and other high-$z$ galaxies, both HZ10-W and the Bridge clearly
depart from the observed sequence of global galaxies the IRX-$\beta_{\rm UV}$
relation. According to theoretical models from previous studies, our results
suggest that the UV-emission in HZ10-W and the Bridge may be strongly
attenuated by a dust screen in between young stellar populations and the
observer.

As mentioned previously, $T_{\rm SED}$ estimations (and related quantities)
are very sensitive to the completeness of the integrated FIR continuum fluxes
set available to model the dust SED. Complementary continuum fluxes
measurements at higher frequencies than those covered by ALMA Band 9 data
could thus allow us to get better constraints on the dust temperature of the
different components in HZ10. In particular, it is necessary to check
possibles scenarios where the dust peak is located at higher frequencies than
those covered by Band 9. ALMA Band 10, centered at $\lambda_{\rm Band10}=$350
$\mu$m, is a good candidate to address this issue for $z\sim 4-6$ galaxies.
Since it gives spectral coverage in the wavelength range $\sim[45,60]$ $\mu$m
at the redshift of HZ10, ALMA Band 10 provides a valuable option for a better
characterization of dust continuum emission.

Upcoming ALMA-CRISTAL studies will analyze the [CII], kinematics and
morphologies of HZ10, and the variations of the IRX-$\beta_{\rm UV}$ relation
among CRISTAL galaxies in more detail (e.g., Telikova et al. in prep.; Ikeda
et al. in prep.; Killi et al. in prep.). In addition, future ALMA Band 10
continuum and JWST (e.g., Jones et al. 2024) data will help us to get better
constraints on the dust emission and metallicity of HZ10, therefore allowing
us to break down the potential degeneracy of our dust SED temperature
estimations and related physical quantities. In addition, ALMA band 10 data
would allow us to trace the small $T_{\rm SED}$ differences and reducing their
uncertainties in the HZ10 system.

###### Acknowledgements.

V. V. acknowledges support from the ALMA-ANID Postdoctoral Fellowship under
the award ASTRO21-0062. 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 BASAL projects
FB210003. R. I. is supported by Grants-in-Aid for Japan Society for the
Promotion of Science (JSPS) Fellows (KAKENHI Grant Number 23KJ1006). K. T.
acknowledges support from JSPS KAKENHI grant No. 23K03466. M. K. was supported
by the ANID BASAL project FB210003. K. T. was supported by ALMA ANID grant
number 31220026 and by the ANID BASAL project FB210003. M. R. acknowledges
support from project PID2020-114414GB-100, financed by
MCIN/AEI/10.13039/501100011033. M. S. was financially supported by Becas-ANID
scolarship #21221511, and also acknowledges ANID BASAL project FB210003. M. A.
acknowledges support from FONDECYT grant 1211951, and ANID BASAL project
FB210003. R. J. A. was supported by FONDECYT grant number 1231718 and by the
ANID BASAL project FB210003. R.B. acknowledges support from an STFC Ernest
Rutherford Fellowship [grant number ST/T003596/1]. This paper makes use of the
following ALMA data: ADS/JAO.ALMA #2022.1.00678.S, ADS/JAO.ALMA
#2019.1.01075.S, ADS/JAO.ALMA #2018.1.00348.S, ADS/JAO.ALMA #2015.1.00388.S.
This research is based on observations made with the NASA/ESA Hubble Space
Telescope obtained from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5–26555. These observations are associated with
program 13641 (P.I.: Peter Capak). The National Radio Astronomy Observatory is
a facility of the National Science Foundation operated under cooperative
agreement by Associated Universities, Inc.

Software: Astropy (Astropy Collaboration et al., 2018), MatPlotLib (Hunter,
2007), NumPy (Harris et al., 2020), PYAUTOGALAXY (Nightingale et al., 2023),
SciPy (Virtanen et al., 2020), seaborn (Waskom, 2021), Scikit-learn (Pedregosa
et al., 2011).

## References

  * Abel et al. (2009) Abel, N. P., Dudley, C., Fischer, J., Satyapal, S., & van Hoof, P. A. M. 2009, ApJ, 701, 1147
  * Akins et al. (2022) Akins, H. B., Fujimoto, S., Finlator, K., et al. 2022, ApJ, 934, 64
  * Algera et al. (2024) Algera, H. S. B., Inami, H., Sommovigo, L., et al. 2024, MNRAS, 527, 6867
  * Álvarez-Márquez et al. (2016) Álvarez-Márquez, J., Burgarella, D., Heinis, S., et al. 2016, A&A, 587, A122
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bakx et al. (2021) Bakx, T. J. L. C., Sommovigo, L., Carniani, S., et al. 2021, MNRAS, 508, L58
  * Béthermin et al. (2015) Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113
  * Bisbas et al. (2022) Bisbas, T. G., Walch, S., Naab, T., et al. 2022, ApJ, 934, 115
  * Bouwens et al. (2020) Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112
  * Bouwens et al. (2016) Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72
  * Bowler et al. (2024) Bowler, R. A. A., Inami, H., Sommovigo, L., et al. 2024, MNRAS, 527, 5808
  * Briggs (1995) Briggs, D. S. 1995, PhD thesis, New Mexico Institute of Mining and Technology
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Cano-Díaz et al. (2016) Cano-Díaz, M., Sánchez, S. F., Zibetti, S., et al. 2016, ApJ, 821, L26
  * Capak et al. (2015) Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455
  * CASA Team et al. (2022) CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501
  * Casey (2012) Casey, C. M. 2012, MNRAS, 425, 3094
  * Colombo et al. (2020) Colombo, D., Sanchez, S. F., Bolatto, A. D., et al. 2020, A&A, 644, A97
  * Contursi et al. (2013) Contursi, A., Poglitsch, A., Graciá Carpio, J., et al. 2013, A&A, 549, A118
  * da Cunha et al. (2013) da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13
  * da Cunha et al. (2021) da Cunha, E., Hodge, J. A., Casey, C. M., et al. 2021, ApJ, 919, 30
  * da Cunha et al. (2015) da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110
  * De Looze et al. (2014) De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62
  * Díaz-Santos et al. (2017) Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32
  * Díaz-Santos et al. (2013) Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68
  * Faisst et al. (2017) Faisst, A. L., Capak, P. L., Yan, L., et al. 2017, ApJ, 847, 21
  * Faisst et al. (2020) Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020, MNRAS, 498, 4192
  * Ferrara et al. (2017) Ferrara, A., Hirashita, H., Ouchi, M., & Fujimoto, S. 2017, MNRAS, 471, 5018
  * Ferrara et al. (2022) Ferrara, A., Sommovigo, L., Dayal, P., et al. 2022, MNRAS, 512, 58
  * Fischera et al. (2003) Fischera, J., Dopita, M. A., & Sutherland, R. S. 2003, ApJ, 599, L21
  * Fudamoto et al. (2020) Fudamoto, Y., Oesch, P. A., Magnelli, B., et al. 2020, MNRAS, 491, 4724
  * Fudamoto et al. (2022) Fudamoto, Y., Smit, R., Bowler, R. A. A., et al. 2022, ApJ, 934, 144
  * Fujimoto et al. (2019) Fujimoto, S., Ouchi, M., Ferrara, A., et al. 2019, ApJ, 887, 107
  * Fujimoto et al. (2020) Fujimoto, S., Silverman, J. D., Bethermin, M., et al. 2020, ApJ, 900, 1
  * Gallerani et al. (2010) Gallerani, S., Maiolino, R., Juarez, Y., et al. 2010, A&A, 523, A85
  * Gallerani et al. (2018) Gallerani, S., Pallottini, A., Feruglio, C., et al. 2018, MNRAS, 473, 1909
  * Ginolfi et al. (2020) Ginolfi, M., Jones, G. C., Béthermin, M., et al. 2020, A&A, 633, A90
  * Glatzle et al. (2022) Glatzle, M., Graziani, L., & Ciardi, B. 2022, MNRAS, 510, 1068
  * Graciá-Carpio et al. (2011) Graciá-Carpio, J., Sturm, E., Hailey-Dunsheath, S., et al. 2011, ApJ, 728, L7
  * Harikane et al. (2020) Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hashimoto et al. (2019) Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71
  * Heinis et al. (2013) Heinis, S., Buat, V., Béthermin, M., et al. 2013, MNRAS, 429, 1113
  * Helou et al. (1988) Helou, G., Khan, I. R., Malek, L., & Boehmer, L. 1988, ApJS, 68, 151
  * Herrera-Camus et al. (2015) Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, ApJ, 800, 1
  * Herrera-Camus et al. (2021) Herrera-Camus, R., Förster Schreiber, N., Genzel, R., et al. 2021, A&A, 649, A31
  * Herrera-Camus et al. (2018a) Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2018a, ApJ, 861, 94
  * Herrera-Camus et al. (2018b) Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2018b, ApJ, 861, 95
  * Hirashita et al. (2005) Hirashita, H., Nozawa, T., Kozasa, T., Ishii, T. T., & Takeuchi, T. T. 2005, MNRAS, 357, 1077
  * Hodge & da Cunha (2020) Hodge, J. A. & da Cunha, E. 2020, Royal Society Open Science, 7, 200556
  * Howell et al. (2010) Howell, J. H., Armus, L., Mazzarella, J. M., et al. 2010, ApJ, 715, 572
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Inami et al. (2022) Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126
  * Issa et al. (1990) Issa, M. R., MacLaren, I., & Wolfendale, A. W. 1990, A&A, 236, 237
  * Jones et al. (2024) Jones, G. C., Bunker, A. J., Telikova, K., et al. 2024, arXiv e-prints, arXiv:2405.12955
  * Knudsen et al. (2016) Knudsen, K. K., Richard, J., Kneib, J.-P., et al. 2016, MNRAS, 462, L6
  * Lai et al. (2023) Lai, T. S. Y., Armus, L., Bianchin, M., et al. 2023, ApJ, 957, L26
  * Le Fèvre et al. (2020) Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1
  * Liang et al. (2021) Liang, L., Feldmann, R., Hayward, C. C., et al. 2021, MNRAS, 502, 3210
  * Liang et al. (2019) Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397
  * Lisenfeld & Ferrara (1998) Lisenfeld, U. & Ferrara, A. 1998, ApJ, 496, 145
  * Lupi et al. (2019) Lupi, A., Volonteri, M., Decarli, R., et al. 2019, MNRAS, 488, 4004
  * Ma et al. (2019) Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844
  * Madden et al. (2013) Madden, S. C., Rémy-Ruyer, A., Galametz, M., et al. 2013, PASP, 125, 600
  * Magdis et al. (2012) Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6
  * Magnelli et al. (2014) Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86
  * Maiolino et al. (2004) Maiolino, R., Schneider, R., Oliva, E., et al. 2004, Nature, 431, 533
  * Malhotra et al. (2001) Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766
  * Malhotra et al. (2017) Malhotra, S., Rhoads, J. E., Finkelstein, K., et al. 2017, ApJ, 835, 110
  * McKinney et al. (2021) McKinney, J., Armus, L., Pope, A., et al. 2021, ApJ, 908, 238
  * McLure et al. (2018) McLure, R. J., Dunlop, J. S., Cullen, F., et al. 2018, MNRAS, 476, 3991
  * Meurer et al. (1999) Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64
  * Meurer et al. (1995) Meurer, G. R., Heckman, T. M., Leitherer, C., et al. 1995, AJ, 110, 2665
  * Misselt et al. (1999) Misselt, K. A., Clayton, G. C., & Gordon, K. D. 1999, ApJ, 515, 128
  * Mitsuhashi et al. (2023a) Mitsuhashi, I., Harikane, Y., Bauer, F. E., et al. 2023a, arXiv e-prints, arXiv:2311.16857
  * Mitsuhashi et al. (2023b) Mitsuhashi, I., Tadaki, K.-i., Ikeda, R., et al. 2023b, arXiv e-prints, arXiv:2311.17671
  * Nath et al. (1999) Nath, B. B., Sethi, S. K., & Shchekinov, Y. 1999, MNRAS, 303, 1
  * Newville et al. (2015) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2015, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python
  * Newville et al. (2016) Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python, Astrophysics Source Code Library, record ascl:1606.014
  * Nightingale et al. (2021) Nightingale, J., Hayes, R., & Griffiths, M. 2021, The Journal of Open Source Software, 6, 2550
  * Nightingale et al. (2023) Nightingale, J. W., Smith, R. J., He, Q., et al. 2023, MNRAS, 521, 3298
  * Pak et al. (1998) Pak, S., Jaffe, D. T., van Dishoeck, E. F., Johansson, L. E. B., & Booth, R. S. 1998, ApJ, 498, 735
  * Pavesi et al. (2016) Pavesi, R., Riechers, D. A., Capak, P. L., et al. 2016, ApJ, 832, 151
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pizzati et al. (2020) Pizzati, E., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 495, 160
  * Popping et al. (2017) Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315
  * Posses et al. (2024) Posses, A., Aravena, M., González-López, J., et al. 2024, arXiv e-prints, arXiv:2403.03379
  * Reddy et al. (2018) Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56
  * Reuter et al. (2020) Reuter, C., Vieira, J. D., Spilker, J. S., et al. 2020, ApJ, 902, 78
  * Riechers et al. (2020) Riechers, D. A., Hodge, J. A., Pavesi, R., et al. 2020, ApJ, 895, 81
  * Saintonge et al. (2016) Saintonge, A., Catinella, B., Cortese, L., et al. 2016, MNRAS, 462, 1749
  * Sargsyan et al. (2012) Sargsyan, L., Lebouteiller, V., Weedman, D., et al. 2012, ApJ, 755, 171
  * Schreiber et al. (2018) Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30
  * Sersic (1968) Sersic, J. L. 1968, Atlas de Galaxias Australes
  * Shibuya et al. (2015) Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15
  * Skibba et al. (2011) Skibba, R. A., Engelbracht, C. W., Dale, D., et al. 2011, ApJ, 738, 89
  * Smith et al. (2017) Smith, J. D. T., Croxall, K., Draine, B., et al. 2017, ApJ, 834, 5
  * Solimano et al. (2024) Solimano, M., González-López, J., Aravena, M., et al. 2024, arXiv e-prints, arXiv:2401.04919
  * Solomon & Vanden Bout (2005) Solomon, P. M. & Vanden Bout, P. A. 2005, ARA&A, 43, 677
  * Sommovigo et al. (2020) Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 497, 956
  * Sommovigo et al. (2022) Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122
  * Spilker et al. (2016) Spilker, J. S., Marrone, D. P., Aravena, M., et al. 2016, ApJ, 826, 112
  * Stacey et al. (2010) Stacey, G. J., Hailey-Dunsheath, S., Ferkinhoff, C., et al. 2010, ApJ, 724, 957
  * Tsukui et al. (2023) Tsukui, T., Iguchi, S., Mitsuhashi, I., & Tadaki, K. 2023, Journal of Astronomical Telescopes, Instruments, and Systems, 9, 018001
  * U et al. (2012) U, V., Sanders, D. B., Mazzarella, J. M., et al. 2012, ApJS, 203, 9
  * Viero et al. (2022) Viero, M. P., Sun, G., Chung, D. T., Moncelsi, L., & Condon, S. S. 2022, MNRAS, 516, L30
  * Villanueva et al. (2024) Villanueva, V., Bolatto, A. D., Vogel, S. N., et al. 2024, ApJ, 962, 88
  * Villanueva et al. (2017) Villanueva, V., Ibar, E., Hughes, T. M., et al. 2017, MNRAS, 470, 3775
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Voit (1992) Voit, G. M. 1992, ApJ, 399, 495
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021
  * Whitaker et al. (2012) Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29
  * Witstok et al. (2022) Witstok, J., Smit, R., Maiolino, R., et al. 2022, MNRAS, 515, 1751