Spatially resolved dust properties and quasar-galaxy decomposition of a hyper-luminous infrared galaxy at z = 4.4

ABSTRACT

We report spatially resolved dust properties of the quasar host galaxy BRI1335−0417 at redshift z = 4.4 constrained by the Atacama Large Millimetre/submillimetre Array observations. The dust temperature map, derived from a greybody fit to rest frame 90 and 161 μm continuum images, shows a steep increase towards the centre, reaching 57.1 ± 0.3 K and a flat median profile at the outer regions of ∼38 K. Image decomposition analysis reveals the presence of a point source in both dust continuum images spatially coincident with the highest temperature peak and the optical quasar position, which we attribute to warm dust heated by an active galactic nucleus (AGN). We show that a model including this warm component along with cooler dust heated by star formation describes the global spectral energy distribution better than a single-component model, with dust temperatures of 87.1|$^{+34.1}_{-18.3}$| K (warm component) and 52.6|$^{+10.3}_{-11.0}$| K (cold component). The star-formation rate (SFR) estimated from the cold dust component is |$1700_{-400}^{+500}\ \mathrm{M}_\odot$| yr⁻¹, a factor of three smaller than previous estimates due to a large AGN contribution (⁠|$53^{+14}_{-15}$| per cent). The unresolved warm dust component also explains the steep temperature gradient, as the temperature profile derived after the point source subtraction is flat. The point source subtraction also reduces the estimated central SFR surface density Σ_SFR by over a factor of three. With this correction, spatially resolved measurements of Σ_SFR and the surface gas mass density Σ_gas form a roughly linear sequence in the Kennicutt–Schmidt diagram with a constant gas depletion time of 50–200 Myr. The demonstrated AGN-host galaxy decomposition reveals the importance of spatially resolved data for accurate measurements of quasar host galaxy properties, including dust temperature, SFRs, and size.

quasars: individual, galaxies: starburst, galaxies: ISM, galaxies: disc, galaxies: spiral, galaxies: bulges

1 INTRODUCTION

Star-formation activity in the Universe peaked at redshift 2 < z < 4 marking a critical stage of galaxy formation (Shapley 2011; Madau & Dickinson 2014), during which more massive galaxies are thought to form their stellar mass earlier (Renzini 2006). Hyper-luminous infrared galaxies (HyLIRGs, defined as galaxies with L_IR > 10¹³ L_⊙) are among the most extreme galaxies found during this epoch, and are thought to host star-formation rates (SFRs) of over 1000 M_⊙ yr⁻¹. Such objects are suggested to form a majority of stellar mass in the most massive elliptical galaxies in the local Universe over 100 Myr (Narayanan et al. 2015). This extraordinary phase of star formation and subsequent quenching largely decides the fate of early forming massive galaxies and their central black holes (BHs). In this context, it represents the most important life event of massive galaxy formation. However, we do not yet understand the detailed driving mechanism of this extreme phase of star formation.

One possibility is the gas-rich major merger paradigm, which suggests dynamic evolution of the system with rapid gas infall towards the centre driven by collisions or external gravitational torques of merging galaxies (Sanders & Mirabel 1996). The rapid radial inflow then triggers a nuclear starburst and feeds gas on to the central BH activating an active galactic nucleus (AGN). The initially dust-obscured AGN removes gas and dust within ∼100 Myr after the peak SFR (Davies et al. 2007; Hopkins 2012), at which point the AGN becomes visible as an optical quasar (Hopkins et al. 2008a, b). The powerful winds driven by the AGN and the burst of star formation eventually eject the remaining galactic gas and prevent subsequent accretion, leading to the rapid cessation of star formation (Silk & Rees 1998; Somerville & Davé 2015). An alternative to gas-rich major mergers suggested by recent studies is that continuous rapid accretion of gas and satellite galaxies from the large gas reservoir of the cosmic web makes an important contribution to fuelling starbursts and the growing BHs (McAlpine et al. 2019; Umehata et al. 2019; Mitsuhashi et al. 2021). Large gas fractions produced by rapid accretion cause the galactic disc to become gravitationally unstable and form non-axisymmetric substructures (such as clumps, a bar, or spiral arms; Hodge et al. 2019), which provide gravitational torques to transfer angular momentum and induce a further inflow of gas (violent disc instability; Dekel & Burkert 2014; Inoue et al. 2016).

In either scenario, the formation of the central BH and the build-up of stellar mass in the surrounding galaxy are co-regulated, with AGN and star-formation feedback acting together to produce the observed tight correlation between central BH masses and host galaxy properties such as central velocity dispersion and bulge stellar luminosity (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000). This remarkably tight correlation may encapsulate a fossil record of their co-evolution (Kormendy & Ho 2013). Early, rapidly evolving galaxies at redshift z > 2 are presumably on their way to establishing this BH-host relationship (Izumi et al. 2019; Pensabene et al. 2020). Studying them, therefore, provides an important probe of the interconnection of AGN and the surrounding star formation.

Our understanding of such starburst galaxies¹ is primarily based on their SFR. This quantity has been derived by assuming the measured far-infrared (FIR; |$\lambda_{\mathrm{rest}}\ \gt\ 100\ {\mu\rm{m}}$|⁠) emission is mainly due to the cold dust component heated by young massive stars with negligible contribution from the warm dust component heated by AGN [which typically dominates at rest frame mid-infrared (MIR), |$5\ \lt\ \lambda _{\mathrm{rest}}\ \lt\ 30\ {\mu\rm{m}}$|⁠, Mullaney et al. 2011; Stalevski et al. 2016; Hönig & Kishimoto 2017]. FIR photometric data are usually fitted with greybody functions to constrain the dust mass M_dust, temperature T_dust, and the power-law slope at a longer wavelength (opacity index β_dust), which together determine the functional shape. The total infrared (TIR) luminosity L_TIR is then derived by integrating the function over the wavelength range of 8–1000 |${\mu\rm{m}}$|⁠, and the SFR is derived using an L_TIR to SFR calibration (Kennicutt 1998). If there are photometric data points available over rest-frame MIR to FIR, the standard procedure is to decompose the spectral energy distribution (SED) into warm and cold dust components by multi-component SED fitting and to estimate the SFR after removing the contribution from the warm dust presumably heated by AGN (e.g. Farrah et al. 2003; Kirkpatrick et al. 2012; Leipski et al. 2014; Kokorev et al. 2021).

However, uncertainty regarding the potential contribution of AGN heating to L_TIR has hindered the understanding the exact nature of high-redshift starburst galaxies. There is growing observational evidence that the AGN contribution to L_TIR is not negligible in all galaxies, but instead increases with total galactic luminosity (Nardini et al. 2010; Yuan, Kewley & Sanders 2010; Alonso-Herrero et al. 2012; Stanley et al. 2017). Therefore, particularly for the brightest population of HyLIRGs, the nature of the heating source for the thermal dust radiating at FIR wavelengths remains an open question (McKinney et al. 2021; Di Mascia et al. 2023)—Does the AGN dominate only for the central region? Or does it heat the entire galaxy? (Symeonidis & Page 2021). To date, most attempts to separate these components have relied on spatially unresolved data, for which estimates of the AGN contribution depend on assumptions (e.g. spectral templates) made during SED modelling. Solving the problem requires a thorough investigation of spatially resolved dust properties (e.g. temperature, mass, and optical depth) and spatial separation of the star formation-dominated host galaxy from the AGN-dominated central region. The need for spatially resolved studies is underscored by the case of the heavily obscured nearby galaxy Arp 220, for which a significant fraction of L_TIR (33 per cent) originates from a compact region with a dust temperature of T_dust = 200 K and a radius of 15 pc, suggesting a significant AGN contribution to L_TIR (Scoville et al. 2017).

Spatially resolved measurements such as those for Arp 220 are however challenging at high redshift. Deriving the temperature distribution as a function of position requires multiple continuum band images near the greybody peak. The most commonly used instruments for studying the thermal dust properties in the previous decade were the Herschel Space Observatory (Herschel) photometers PACS (from 70 to 160 microns) and SPIRE (from 250 to 500 microns), which provided a spatial resolution [full width at half-maximum (FWHM) = 36 arcsec] sufficient only to resolve nearby galaxies (Galametz et al. 2012). Recently, however, the Atacama Large Millimeter/submillimeter Array (ALMA) has started to provide frequency coverage near the peak of the greybody function for z > 4 galaxies with enough angular resolution to permit spatially resolved analysis. While spatially resolved pixel-by-pixel temperature maps for such high-redshift galaxies are now attainable in principle, they have not been achieved so far in practice. This is mainly because high-angular resolution observations with a sufficient signal-to-noise ratio (S/N) in at least two frequency bands² are demanding even with ALMA for relatively faint high-redshift galaxies. A few exceptions exist where authors have constrained the temperature gradient within individual galaxies: Shao et al. (2022) derive one-dimensional (1D) radial profiles of the dust temperature, mass, and optical depth for a quasar host galaxy at redshift z = 6 using azimuthally averaged profiles at two FIR continuum bands, while Akins et al. (2022) derive a pixel-by-pixel temperature map, which shows evidence for a temperature gradient, in a strongly lensed star-forming galaxy at redshiftz = 7.

In this paper, we derive spatially resolved temperature maps for BRI1335−0417, a quasar host galaxy at a redshift of z = 4.4704 (Guilloteau et al. 1997), 1.4 Gyr after the big bang. This galaxy is one of the brightest unlensed submillimetre sources known at z > 4 (Jones et al. 2016), and is classified as a HyLIRG with an extraordinary infrared luminosity of 3.1 × 10¹³ L_⊙(Carilli et al. 2002). It was originally identified as an optical quasi stellar object (QSO) by optical imaging and spectroscopy from the Automatic Plate Measuring survey (Irwin, McMahon & Hazard 1991; Storrie-Lombardi et al. 1996). The optical QSO position is RA = 204.514 232 850 ± 0.000 000 092 deg, Dec = −4.543 050 970 ± 0.000 000 055 deg (International Celestial Reference System, ICRS; Gaia Collaboration 2016, 2021). The galaxy hosts a BH with a mass suggested to be 10^9.77 M_⊙, shining at |$\sim\!36\ {{\ \rm per\ cent}}$| of its Eddington luminosity (Shields et al. 2006). The SFR of the galaxy was estimated to be 5040 ± 1300M_⊙ yr⁻¹ from modelling of the spatially unresolved SED with thermal dust and synchrotron emission components (Wagg et al. 2014), a rate high enough to deplete its total molecular gas reservoir of ∼10¹¹ M_⊙ yr⁻¹ estimated by CO(2-1) line observation (Jones et al. 2016) in only ∼20 Myr.

Although such short depletion times are commonly attributed to gas-rich major mergers, the morphology and kinematics of the galaxy appear to be inconsistent with this scenario. Spatially resolved [C ii] and dust continuum observations (∼1.3 kpc resolution) show clear evidence for a rotating disc and spiral morphology (Tsukui & Iguchi 2021) with a further analysis of the gas motion indicating the presence of a compact mass structure in the centre of the galaxy. Indeed, BRI1335–0417 is the highest redshift galaxy thus far to show a spiral morphology. The spiral structure is visible in [C ii] line and dust emission, indicating ongoing star formation. BRI1335–0417 is not the only extremely luminous high-redshift galaxy to show a surprisingly quiet and well-ordered morphology. Contrary to earlier views that cold discs and spirals only begin to appear at z < 2 (Elmegreen & Elmegreen 2014), ALMA and JWST observations now suggest a surprisingly earlier epoch of galaxy settling at z ∼ 2–4, with cold gas discs found at z ∼ 4 (ALMA, Neeleman et al. 2020; Rizzo et al. 2020, 2021; Lelli et al. 2021; Tsukui & Iguchi 2021), stellar spiral structure detected in passive galaxies at 1 < z < 3 (JWST; Fudamoto, Inoue & Sugahara 2022), and grand design barred spirals already in place at z ∼2 (JWST; Guo et al.2023).

This contradiction between the ultrashort depletion time and the morphology of BRI1335–0417 suggests that it is worth revisiting the high SFR previously estimated from an unresolved SED by using the spatially resolved data to which we now have access. Doing so may provide insight not just on this particular source, but more broadly on the driving mechanism of high-z starbursts, early build-up of BHs and stellar bulges, and evolutionary links from the cold gas discs at z = 4 to stellar spirals at z ∼ 2–3. With this motivation in mind, this paper aims to (i) derive a resolved dust temperature map and clarify the heating source of the dust, (ii) separate the warm dust heated by AGN and cold dust associated with the host galaxy, and (iii) study the SFR distribution under the effects of the central quasar in BRI1335–0417. We present new Band 9 (⁠|$\sim\!484\ {\mu\rm{m}}$|⁠) and Band 4 (⁠|$\sim\!2080\ {\mu\rm{m}}$|⁠) ALMA observations in addition to the earlier Band 7 (⁠|$\sim\!869\ {\mu\rm{m}}$|⁠) data presented in Tsukui & Iguchi (2021; PI: J. González López). Multiple continuum images provide a number of resolution elements over the galaxy, making it possible to investigate the spatially resolved physical properties of dust (such as temperature and optical depth). The paper is organized as follows. In Section 2, we describe the observations and the data reduction. In Section 3, we present the results of dust SED modelling applied to each pixel of the spatially resolved continuum images. In Section 3.2, we describe our image decomposition method (point source and host galaxy) followed by the panchromatic SED modelling of BRI1335–0417 with the help of the decomposition results. In Section 5, we present the discussion of the result. In Section 6, we summarise and conclude the paper.

Throughout the paper, we adopt a flat lambda cold dark matter cosmology with a present-day Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹, and a density parameter of pressureless matter Ω_M = 0.3, providing an angular size distance D_A = 1375 Mpc and a luminosity distance D_L = 40 202 Mpc at the redshift of 4.4074.

2 OBSERVATIONS AND DATA REDUCTION

2.1 ALMA imaging

ALMA Band 4, Band 7, and Band 9 observations of BRI1335–0417 were carried out as part of the programs #2017.1.00394.S and #2018.1.01103.S (PI: J. González López). These observations targeted emission lines, including CO(7-6) (Band 4) and [C ii] (Band 7), along with underlying continuum emission at observing wavelength 2080 |${\mu\rm{m}}$| in Band 4, 869 |${\mu\rm{m}}$| in Band 7, and 484 |${\mu\rm{m}}$| in Band 9, corresponding to rest frame 385, 161, and 90 |${\mu\rm{m}}$|⁠, respectively. We performed standard calibration and data reduction using the Common Astronomy Software Application (casa; CASA Team et al. 2022) pipeline.³ The flux and bandpass were calibrated using the quasars J1337−1257 for Band 4 and Band 7, and J1256−0547 for Band 9 data. The phase was calibrated using the quasars J1332−0509 for Band 4 and J1336−0829 for Band 7 and Band 9. We identified line-free channels using the hif_findcont task in casa and additionally removed channels affected by atmospheric absorption. We estimated the flux density of the underlying continuum emission by fitting a linear function to the identified line-free channels and subtracted it from the data cube in the visibility plane. We then imaged the line-free channels to produce continuum images and the continuum-subtracted data to produce emission line cubes. The visibility data were weighted using a Briggs weighting scheme with a robust parameter of 1.0 for the CO(7-6) line data cube to improve the sensitivity and 0.5 for others, which provides a good compromise of spatial resolution and sensitivity (Briggs 1995). The resulting images and cube information such as angular resolution and point source sensitivity are summarised in Table 1. The emission line flux maps are made by summing the velocity channels from −400 to 400 km s⁻¹. The mean velocity and the velocity dispersion of the emission line are extracted by fitting a single Gaussian to the spectrum at each pixel. We only use pixels where the signal is detected more than 4σ over at least four channels.

Table 1.

ALMA observation data summary. The root mean square (RMS) noise level is measured using the emission-free region of the images. Synthesized beams are reported as [major axis, minor axis, and position angle (PA)]. RMS noise in the unit of SFR is calculated by assuming optically thin greybody emission with a dust temperature of 39 K, a dust emissivity index β = 2.14 (see Section 3.2), and a Chabrier (2003) initial mass function.

Data	RMS noise	Synthesized beam	Subtended area per beam	RMS noise in SFR
	(μJy beam−1)	maj (arcsec)/min (arcsec)/PA (deg)	kpc2 beam−1	M⊙ yr−1 kpc−2
|$I_{\lambda _{\mathrm{rest}}\ =\ 385\ \mu \mathrm{m}}$| (Band 4)	6.58	0.21/0.14/85	1.52	11.8
|$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| (Band 7)	25.9	0.19/0.16/75	1.56	7.5
|$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| (Band 9)	186	0.21/0.16/−53	1.69	5.4
|$\left[\mathrm{C}~{\small II}\right]$|	389	0.19/0.16/82	1.59	–
CO(7-6)	101	0.23/0.16/84	1.93	–

In this paper, the same phase centre and pixel size are used to image the visibility data so that all data products have the same pixel coordinates, with absolute positional accuracy of 10–20 mas.⁴ Before deriving the physical parameters using two or more images (e.g. dust SED modelling), we convolved images to have the same spatial resolution with the smallest common beam size, using common_beam in the spectral-cube package. We similarly convolve all other data products to the same resolution, so that physical parameters derived (e.g. velocity dispersions derived from emission line cubes) are resolution-matched. The matched resolution of the [C ii] cube, rest frame 161 |${\mu\rm{m}}$|⁠, and 90 |${\mu\rm{m}}$| continuum images is 0.21 arcsec × 0.18 arcsec at PA = 114 deg (beam area of 1.92 kpc² and effective radius of 0.78 kpc). Analysis and visualization of data in this paper are done at this resolution except for the radial profile analysis carried out in Section 3.5, which includes the poorer resolution CO(7-6) data. For the comparison with CO(7-6) all data are convolved to the common resolution of 0.24 arcsec × 0.18 arcsec at PA = 94 deg (beam area of 2.3 kpc² and effective radius of 0.85 kpc). The effects of the beam on our spatially resolved SED analysis are also discussed in Section 3.7.

2.2 Hubble Space Telescope imaging

We retrieved the STIS/50CCD image of BRI1335–0417 from the Hubble Space Telescope (HST) archive. The data were taken as a part of program GO-8572 (PI: L. Storrie-Lombardi) in January 2001, consisting of the four dither exposures with each integration time of 645 s. The detector has broad sensitivity from 2000 to 10 300 Å. Therefore, we only use the data for visualization purposes to assist the interpretation of the other optical photometric data of the galaxy shown in Table A1. We processed the fully calibrated subexposure images from the HST archive applying the geometrical distortion and world coordinate system corrections, and then drizzling on to the final pixel grid with the pixel size of 0.05 arcsec. The target acquisition is based on the Guide Star Catalog GSC 1.0, which is expected to have a pointing accuracy of 1–2 arcsec in worst cases. We calibrate the pointing offset of the image by translation to match the image peak of the BRI1335–0417 to the corresponding Gaia coordinate. We identify the BRI1335–0417 in the image using three Gaia sources available in the field of view, including BRI1335–0417. The calibrated offset length of 0.8 arcsec is within the expected pointing accuracy. Fig. 1 shows the 50CCD image overlaid with the contours of the ALMA Band 7 continuum image. The optical image is consistent with the point source emission from the AGN [see Fig. A1 for a comparison of the radial distribution of the emission and the point spread function (PSF)]. In contrast, the ALMA Band 7 image clearly resolves the extended host galaxy.

Figure 1.

HST STIS/50CCD image overlain with the ALMA Band 7 image contour. The optical band emission is consistent with the point source. In contrast, the ALMA Band 7 image provides the extended morphology of the host galaxy.

3 SPATIALLY RESOLVED DUST PROPERTIES

3.1 The resolved FIR continuum and [C ii] emission

In Fig. 2, we show continuum images of rest frame 385, 161, 90 |${\mu\rm{m}}$|⁠, which we denote |$I_{\lambda _{\mathrm{rest}}\ =\ 385\ \mu \mathrm{m}}$|⁠, |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|⁠, and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠, respectively, in the paper, together with the total flux and mean velocity and velocity dispersion derived from the [C ii] and CO(7-6) line data. In all images of this paper, the RA and Dec offsets are given relative to the position RA = 204.514 232 deg and Dec = −4.543 055 deg, which coincides with the peak of the continuum images and the quasar position within the positional accuracy of ALMA. The continuum emission is more centrally concentrated than the [C ii] line emission, as the 385, 161, and 90 |${\mu\rm{m}}$| continuum images and their contours are shown in log scale while the [C ii] image is shown in linear scale (we also show linear-scaled continuum images below, in Section 3.2). This difference is commonly seen in other spatially resolved observations of quasar host galaxies (e.g. Shao et al. 2022; Walter et al. 2022). As presented in Tsukui & Iguchi (2021), both |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and [C ii] images show Z-shaped spiral structures. The |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| image also has a disc-like morphology with a major axis similar to that in |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|⁠, but a fainter spiral-arm feature. The |$I_{\lambda _{\mathrm{rest}}\ =\ 385\ \mu \mathrm{m}}$| is less sensitive to SFR which leads to lower S/N, so the emission is not detected in the outer part of the galaxy. The mean [C ii] velocity map is qualitatively consistent with the moment 1 map (intensity-weighted velocity of the spectral line) presented in Tsukui & Iguchi (2021), but provides more accurate measurements than the moment 1 map, which can easily be affected by noise in the emission-free channels. The [C ii] velocity dispersion map shows a triangle-shaped velocity-enhanced region extending from the centre to the north-west direction. Such an asymmetric structure cannot be explained by the velocity gradient of a regularly rotating motion, suggesting that the origin may not be gravitational.

Figure 2.

ALMA observation for QSO source BRI1335–0417 at redshift z ∼ 4.4. The top panels show the continuum images at rest frame 385 |${\mu\rm{m}}$| (Band 4; a), 161 |${\mu\rm{m}}$| (Band 7; b), and 90 |${\mu\rm{m}}$| (Band 9; c) from left to right. The middle panels show the [C ii] integrated flux (d), mean velocity (e), and velocity dispersion maps (f). The bottom panels show the CO(7-6) flux (g), mean velocity (h), and velocity dispersion (i) maps. The mean velocity and velocity dispersion at each pixel are extracted by fitting a single Gaussian to the spectrum. We only use pixels where the signal is detected more than 4σ over at least four channels. The solid line contours in (a), (b), and (c) are plotted for 3|$\sigma \times \sqrt{3}^{n};n=\lbrace 0,1,2,3,...\rbrace$| to 27σ, 140.3σ, and 81σ, respectively. The solid line contours in (d) and (g) are plotted for 3σ + 4n; n = {0, 1, 2, 3,...} to 23σ and 31σ, respectively. Dotted line contours in all flux maps show −3σ. Contours in velocity and velocity dispersion map are shown every 20 km s⁻¹, where plus and minus values are shown in solid and dotted lines, respectively. The size of the synthesized beam (FWHM) is shown in the lower-left corner of each map.

3.2 Dust SED modelling

To derive dust properties from our observed spatially resolved FIR fluxes, we use a greybody function to describe the flux. For a uniform region of dust with an equilibrium temperature T_dust at redshift z, solving the equation of radiative transfer yields a predicted cosmic microwave background (CMB)-subtracted flux within the ALMA beam (Walter et al. 2022):

$$\begin{eqnarray} F_\nu =\Omega _\mathrm{a}\times [B_\nu (T_{\mathrm{dust,z}})-B_\nu (T_{\mathrm{CMB,z}})]\times [1-\exp (-\tau _\nu)](1+z)^{-3}, \\ \end{eqnarray}$$

(1)

where Ω_a is the solid angle of the ALMA synthesized beam in steradians, B_ν is the Planck function, T_dust,z is the dust temperature, and T_CMB,z = 2.73(1 + z) K is the CMB temperature at a given redshift, z = 4.4704 in our case. τ_ν is optical depth, which is related to the dust mass M_dust by τ_ν = κ₀(ν/ν₀)^βM_dustA⁻¹, where |$A=\Omega _\mathrm{a} D_\mathrm{A}^2$| is the physical area subtended by the beam, β is the dust opacity index, and κ₀ = 5.1 cm² g⁻¹ is the dust opacity at a reference frequency ν₀ = 1199.1 GHz (Draine & Li 2007). Note that while the observed flux depends on the actual dust temperature at the redshift of the source, T_dust,z, for convenience we compute and report the intrinsic dust temperature corrected for the CMB heating, T_dust, which is the temperature the dust would have at redshift zero, where CMB heating is much smaller. As described in da Cunha et al. (2013), these two temperatures are related as

$$\begin{eqnarray} T_{\mathrm{dust}, z}=\left\lbrace T_{\mathrm{dust}}^{\beta +4}+T_{\mathrm{CMB}, z=0}^{\beta +4}\left[(1+z)^{\beta +4}-1\right]\right\rbrace ^{\frac{1}{4+\beta }}. \end{eqnarray}$$

(2)

3.3 Global SED of the galaxy

In Fig. 3, we show the total (spatially integrated over the galaxy) SED of BRI1335–0417 with our new flux measurements at rest frame 385, 161, and 90 |${\mu\rm{m}}$|⁠. We measured the integrated flux densities by summing over the 2 arcsec × 2 arcsec region shown in Fig. 2. The available data points are listed in Table A1. We excluded available photometric measurements from the Herschel SPIRE 350 and 500 |${\mu\rm{m}}$| bands because BRI1335–0417 is not sufficiently separated from a nearby bright source at Herschel’s angular resolution (see each photometric image in Fig. A2). To constrain the global dust properties of the galaxy, we fit the SED with a single greybody in the optically thin limit (τ_ν → 0, large Ω_a). In this limit, equation (1) becomes

$$\begin{eqnarray} F_{\nu }=(1+z)D_\mathrm{L}^{-2} \kappa _{\nu }M_{\mathrm{dust}}[B_{\nu }\left(T_{\mathrm{dust,z}}\right)-B_{\nu }\left(T_{\mathrm{CMB,z}}\right)]. \end{eqnarray}$$

(3)

Figure 3.

The SED of BRI1335–0417 fitted with a single greybody in the optically thin limit, with the best-fitting parameters of |$T_{\mathrm{dust}} = 39.0_{-2.9}^{+3.4}$| K and β = 2.14 ± 0.17. The available data points are listed in Table A1, where the greybody fit is performed on bands longer than 160 |${\mu\rm{m}}$| (rest frame 36 |${\mu\rm{m}}$|⁠); we exclude the 160 |${\mu\rm{m}}$| data (the grey point) from the fit because a single greybody fit cannot reproduce it. The TIR luminosity derived from the fit is |$L_{\mathrm{TIR}}=2.7_{-0.3}^{+0.4}\times 10^{13}\ \mathrm{L}_\odot$|⁠, consistent with the previously reported value for this galaxy of 3.1 × 10¹³ L_⊙ (Carilli et al. 2002).

In the fitting, we treat T_dust, M_dust, and β_dust as free parameters. We find that a single greybody fit cannot reproduce the flux observed in the Herschel 160 |${\mu\rm{m}}$| band (rest frame 36 |${\mu\rm{m}}$|⁠), where warm dust is expected to dominate. Therefore, we repeat the fit with the 160 |${\mu\rm{m}}$| band flux excluded; doing so yields a total dust mass of |$1.9^{+0.3}_{-0.3}\times 10^9\ {\rm M}_{\odot }$|⁠, a dust temperature |$T_{\mathrm{dust}} = 39.0_{-2.9}^{+3.4}$| K, and an emissivity index β = 2.14 ± 0.17, which describes the power-law slope of the greybody function in the Rayleigh–Jeans tail (∝λ^{−2 − β} at |$\lambda \gg 100\ {\mu\rm{m}}$|⁠). The derived emissivity index and the computed total FIR luminosity |$L_{\mathrm{TIR}}=2.7_{-0.3}^{+0.4}\times 10^{13}\ \mathrm{L}_\odot$| are consistent with the previously derived values for this galaxy, β = 1.89 ± 0.23 (Wagg et al. 2014) and L_TIR ∼ 3.1 × 10¹³ L_⊙ (Carilli et al. 2002). The gas-to-dust mass ratio is found to be 54.2 ± 9.3 if we adopt a total gas mass derived by Jones et al. (2016) from CO (2-1) using the standard conversion factor for SMGs α_CO = 0.8 M_⊙ pc⁻²/(K km s⁻¹) (Bolatto, Wolfire & Leroy 2013) and a ratio of |$\mathrm{r}_{21}=L_{\mathrm{CO}(2 \rightarrow 1)}^{\prime }/L_{\mathrm{CO}(1 \rightarrow 0)}^{\prime }=0.85$| for SMGs (Carilli & Walter 2013).

The flux density excess, which cannot be captured by a single greybody function in rest frame from ∼8 to |$50\ {\mu\rm{m}}$|⁠, has previously been attributed to the presence of higher temperature components particularly in nuclear regions (Casey, Narayanan & Cooray 2014). Different prescriptions have been employed in the literature to reproduce the excess flux, including the sum of two greybody functions with different dust temperatures (Dunne & Eales 2001; Farrah et al. 2003) or a greybody function with the Wien part replaced by the power law function (Casey et al. 2012). In the former method, the warmer component has been commonly attributed to the AGN-heated dust (Farrah et al. 2003; Kirkpatrick et al. 2012; Leipski et al. 2014; Kokorev et al. 2021), an interpretation we will explore in detail with the spatially resolved data in this study. The latter approach of adding a power-law component is somewhat phenomenological, but has the advantage that it introduces only a single additional free parameter—an advantage that is not negligible when the set of measurements available to fit the data is very small, as it is here.

3.4 Intrinsic dust properties derived from the spatially resolved dust SED

Spatially resolved images at different frequency bands allow us to derive the intrinsic dust properties at each individual pixel. The galaxy is spatially resolved with 14, 32, and 20 resolution elements detected at ≥3σ significance in the rest frame 385 |${\mu\rm{m}}$| (Band 4), 161 |${\mu\rm{m}}$| (Band 7), and 90 |${\mu\rm{m}}$| (Band 9) continuum images, respectively. The rest frame 161 and 90 |${\mu\rm{m}}$| bands are close to the peak of the greybody, and their ratio is therefore sensitive to the dust temperature. The spatially resolved rest frame 385 |${\mu\rm{m}}$| may help to constrain the emissivity index β or the SED slope at a longer wavelength. However, the rest frame 385 |${\mu\rm{m}}$| data has the worst sensitivity in terms of the minimal detectable SFR (see Table 1), and the 3σ detected emission area is smaller than in the other bands. Therefore, we choose to use the best estimate of β = 2.14 ± 0.17 obtained with the global SED modelling as a fiducial value in this paper, and not to include the rest frame 385 |${\mu\rm{m}}$| band for spatially resolved SED modelling. With flux densities measured in two bands and the fiducial dust emissivity β = 2.14 ± 0.17, we can fully constrain the two remaining free parameters Σ_dust and T_dust in equations (1) and (2). We propagate the uncertainty on the spectral index, β = 2.14 ± 0.17, to the derived dust properties via Monte Carlo.

Figs 4a–c show the best-fitting dust temperature, the optical depth of the dust at rest frame |$161\ {\mu\rm{m}}$| (or the dust surface density Σ_dust), and the integrated dust luminosity density maps of the galaxy, respectively. We derive our 1σ confidence intervals on the fit parameters (T_dust and Σ_dust) by Monte Carlo resampling: We repeat the fit 300 times using fluxes randomly drawn from a Gaussian distribution centred on the best-fitting measured value with a dispersion equal to the 1σ observational uncertainty, and a spectral index β drawn from a Gaussian centred at β = 2.14 with a dispersion of 0.17. We report the width of the central 68 per cent of these 300 trials as our 1σ uncertainties on T_dust and Σ_dust. We only show and use pixels where the 68 per cent confidence interval for the dust temperature is smaller than 10 K. Figs 4d–f show how the measured flux densities at rest frame |$161\ {\mu\rm{m}}$| (black dot–dashed line) and |$90\ {\mu\rm{m}}$| (black dashed line) band constrain dust physical parameters Σ_dust and T_dust at fixed β = 2.14, for three example regions: (d; left-hand panel) the region with the lowest temperature, (e; middle panel) the central region that has the highest temperature, and (f; right-hand panel) the region with the second highest temperature. The intersection of the two black lines corresponds to the most likely solution of Σ_dust and T_dust and the resampled distributions are shown in grey points. The dust temperature T_dust and surface density Σ_dust are correlated but still constrained around the intersection of the isoflux lines of the two bands.

Figure 4.

The upper panels show the spatial distribution of the best-fitting dust temperature (a), optical depth at λ_rest = 161 μm or dust mass surface density (b), and surface density of the integrated TIR luminosity over 8–1000 μm, Σ_TIR (L_⊙ kpc⁻²) (c). The lower panels show the isoflux lines of the measured flux density (λ_rest = 161 μm; dot–dashed line, λ_rest = 90 μm; dashed line) in the parameter space Σ_dust and T_dust, for three example regions: (d) the region with the lowest temperature, (e) the central region with the highest temperature, and (f) the region with the second highest temperature. The regions (d–f) are indicated with blue ellipses in panels (a–c) and labels in panel (b). The intersection of the two black lines corresponds to the most likely solution (blue solid lines), while the blue shades enclose the 68th percentile of the resampled distribution (grey points) reflecting the noise in the data and the confidence interval for the fiducial dust emissivity β = 2.14 ± 0.17 derived in the global SED fitting (Fig. 3). Contours in (a) are shown every 5 K from 30 to 55 K. Contours in (b) are shown every 0.1 from 0.1 to 0.5. Contours in (c) are shown for 10¹¹, 10^11.5, 10¹², and 10^12.5 L_⊙ kpc⁻². The size of the synthesized beam (FWHM) is shown in the lower-left corner of (a–c).

The temperature map reveals a central high-temperature region surrounded by a low-temperature region. The highest temperature peak coincides with the peak position of both dust continuum emission images |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|⁠. The dust temperature of the central pixel is well constrained to be 57.1 ± 0.3 K, which is higher than typical dust temperatures of quasar host galaxies (∼47 K; Beelen et al. 2006), indicating that the central high-temperature peak is likely due to the dust heated by the AGN. The temperature map also shows several high-temperature regions in the outer part of the galaxy. The surface density or optical depth maps show that the central region is moderately optically thick at rest frame 161 μm and optically thick at rest frame 90 μm (τ > 1). The surface density map shows the southern spiral structure, which corresponds to the relatively low temperature <40 K region in the temperature map.

Fig. 5 shows the radial distribution of the temperature, computed from data points extracted from the image (Fig. 4a). We calculate the de-projected radius using the thin disc geometry proposed in Tsukui & Iguchi (2021) (PA 7.6 deg, inclination 37.3 deg), derived from dynamical modelling of [C ii] emission kinematics with prior constraints from the axial ratio of the continuum image and the kinematic PA of [C ii] emission. The derived dust temperature steeply decreases as a function of the galactic radius from the centre to 2.5 kpc. At the outer disc (∼3 kpc), both high- and low-temperature regions can be seen, and the temperature difference is statistically significant, corresponding to the clumpy high-temperature regions seen in the map. At a larger radius, the temperature profile becomes flat with a median temperature of 38 K and a large spread of values from 30 to 45 K, roughly consistent with the typical dust temperature of the high redshift starburst galaxies (Magnelli et al. 2012) at L_TIR ∼ 10¹³ L_⊙. The presence of the high-temperature regions at a larger radius may be common in other quasar host galaxies and is consistent with the observational signature of J0305−3150 at z = 6.6, which shows an increase in the surface brightness ratio |$I_{\lambda _{\mathrm{rest}}\ =\ 175\ \mu \mathrm{m}}/I_{\lambda _{\mathrm{rest}}\ =\ 459\ \mu \mathrm{m}}$| with radius (Li et al. 2022), for which the authors suggested two possibilities: radial increase of the temperature or decrease of the opacity. We discuss the possible origin of the outer-disc high-temperature regions in Section 3.6.

Figure 5.

Radial profile of the temperature of BRI1335–0417. The black points with error bars show temperature values, three data points per beam extracted from the temperature map (Fig. 4a) with a regular grid aligned with the central pixel. The grey-shaded region encloses 68 per cent of all available pixels in radial bins with a width of 1/3 of the beam FWHM. The grey line indicates the median in each bin. The horizontal axis shows the distance from the galactic centre (dust peak position) assuming the disc geometry presented in Tsukui & Iguchi (2021).

3.5 Comparison of the derived dust mass profile and other gas mass tracer profiles

Fig. 6 shows surface gas mass densities derived by three possible gas mass tracers: emission lines of [C ii] and CO(7-6), and the SED-derived dust mass. For the two former methods, we assume a constant emissivity per unit mass, while for the latter we assume a constant dust-to-gas ratio, and we normalize the profiles derived from each method to be equal at ∼2 kpc. Both the dust- and CO-derived gas mass profiles are steeper than that estimated from the [C ii] line. Based on global measurements of galaxies, Zanella et al. (2018) proposed the [C ii] luminosity as a molecular gas mass tracer with a mass-to-light ratio |$M_{\mathrm{gas}}/L_{\mathrm{[C~{\small II}]}}=30\ {\rm M}_{\odot }/\mathrm{L}_{\odot}$|⁠. However, our result suggests that the gas mass profiles estimated from [C ii] do not agree with either CO- or dust-based profiles under the assumption of a spatially constant mass-to-light ratio for [C ii]. Possible explanations for [C ii] showing a flatter profile than CO and dust are (i) [C ii] can be emitted from various gas phases including warm ionized gas and diffuse CO-dark gas, which are likely to extend to regions beyond those that can be traced by CO emission (Pineda et al. 2013), (ii) in the central region the C⁺ abundance decreases because the high surface density shields gas from UV photons, leading most carbon to transition to CO (Narayanan & Krumholz 2017), and (iii) the effect of infrared background radiation as proposed by Walter et al. (2022).

Figure 6.

Radial profile of normalized gas mass surface densities of BRI1335–0417. The black points with error bars and black shade show gas surface density assuming the constant gas-to-dust ratio, extracted in the same way as Fig. 5 from the dust surface density map (Fig. 4b). We also show two additional gas mass estimates from [C ii] (blue shade) and CO(J = 7→6) (red shade) using a constant mass-to-light ratio. All the profiles are normalized to 1 at a radius of ∼2 kpc.

The CO-derived gas mass profile shows a steeper profile than the dust-derived one, especially in the central ∼1 kpc. The excess in the centre may be due to the radially different excitation conditions, where the quasar radiation dominates in the central region and intense star formation dominates in the outer part of the galaxy (Carilli & Walter 2013). In addition to the excitation conditions, some of the same mechanisms that flatten the [C ii] profile may steepen the CO one: (i) the dust may trace not only molecular gas but also generally more extended atomic gas (Orellana et al. 2017) and (ii) in the central region the CO abundance may be enhanced by increased shielding against dissociating UV photons (Narayanan & Krumholz 2017). However, in contrast to [C ii], the infrared background radiation effect proposed by Walter et al. (2022) is not significant for CO(7-6) because the dust continuum optical depth at the CO(7-6) frequency is <0.1.

In addition to the normalized profiles shown in Fig. 6, we can also examine the absolute surface densities. If we derive the molecular gas mass profile using α_CO = 0.8 M_⊙ pc⁻²/(K km s⁻¹) and |$\mathrm{r}_{21}=L_{\mathrm{CO}(2 \rightarrow 1)}^{\prime }/ L_{\mathrm{CO}(1 \rightarrow 0)}^{\prime }=0.85$| as we did for our total molecular mass estimate, we automatically obtain the same estimate of the total molecular gas mass as Jones et al. (2016), 1.03 ± 0.07 × 10¹¹ M_⊙. To match the CO-estimated and dust-estimated gas masses at 2 kpc would require a gas-to-dust ratio 44 ± 24, with a relatively large error bar due to the large dispersion of the derived dust mass. This gas-to-dust mass ratio is consistent within the uncertainties with that derived from the whole-galaxy SED fit (see Fig. 3). Further exploration of gas emissivities per unit mass in various tracers is beyond the scope of this work, but this spatially resolved analysis may highlight that assumptions of a constant emissivity break down when applied to spatially resolved data. (See appendix C of Herrera-Camus et al. 2021 for further discussion of [C ii] emissivity per unit mass in the optically thin limit and negligible background emission.)

3.6 The correlation of the temperature map and velocity dispersion map

In Fig. 7, we compare the dust temperature and velocity dispersion maps overlaid with low- and high-temperature peaks. We see that high-temperature peaks coincide with regions of enhanced velocity dispersion extending from the centre to the north-west direction, while low-temperature peaks coincide with low-velocity dispersion regions along the disc major axis. Such anisotropic distributions preferentially aligned with the disc minor axis may indicate that an outflow from the central region heats up the dust and enhances the velocity dispersion. In a nearby luminous infrared galaxy, NGC6240, Saito et al. (2018) find a similar bipolar distribution for the high-velocity dispersion CO emission, which coincides with the spatial distribution of H α, NIR, and X-ray emission. To examine the potential bulk gas motion due to the outflow which deviates from the mean line-of-sight motion in BRI1335–0417, we fit the [C ii] spectrum at each position of the high-temperature peaks. The spectra can be fitted with only a single Gaussian component and do not require any additional components given the current noise level of the data.

Figure 7.

The temperature map (left-hand panel) and velocity dispersion map (right-hand panel) as shown in Fig. 4a and Fig. 2f, respectively. The identified high- and low-temperature peaks (magenta and cyan crosses) are overlaid on the maps. The size of the synthesized beam (FWHM) is shown in the lower-left corner of each panel.

Another possible origin of the correlated enhancement of the temperature and gas velocity dispersion is energy injection by intense star formation. Using adaptive optics assisted Keck observations, Oliva-Altamirano et al. (2018) suggested a spatial correlation between the peaks of the Pa α velocity dispersion and SFR for a relatively low-redshift galaxy sample (z = 0.07 to 0.2). If this scenario is correct, we expect high-velocity dispersion, high-temperature regions to show high gas densities in order to host and sustain the intense star formation; these regions might include giant star-forming clumps or interacting satellites. Neither the derived dust mass surface density (see Fig. 4b) nor the observed [C ii] and CO(7-6) maps (see Fig. 2) show such a distinct dense region, though they may show a broader triangle shape. However, Tadaki et al. (2018) identified 200 pc-scale clumpy structures in a highly star-forming submillimetre galaxy at redshift z = 4, both in the CO line and dust continuum. Whether similar small-scale clumpy structures occupy the temperature-enhanced regions in our target cannot be addressed with the current data due to sensitivity and resolution limits. To further constrain this scenario, higher resolution data capable of tracing denser gas are required.

3.7 Systematic error of the spatially resolved measurements

We confirm that the absolute flux uncertainty does not change the overall temperature distribution. Considering the worst possible case, when the true Band 7 flux and Band 9 flux are 10 per cent higher and lower than the observed values, respectively, the temperature would be underestimated by 2.7 K in the median. In the inverse scenario, where the true Band 7 flux and Band 9 flux are 10 per cent lower and higher than the observed values, the temperature would be overestimated by 2.0 K in the median. Therefore, the systematic error due to the absolute flux uncertainty is estimated to be  + 2.7 K/−2.0 K. The flux normalization errors only induce a pixel-to-pixel variation of 0.67 and 0.49 K, respectively.

We used the fiducial β value β = 2.14 ± 0.17 derived from the integrated SED modelling, which is close to the recently measured values for the high-redshift galaxies: median β = 1.9 ± 0.4 for ALESS well-constrained subsamples (da Cunha et al. 2021) and β = 2.1 ± 0.1 for 47 starburst galaxies (Birkin et al. 2021, reported as a private communication in da Cunha et al. 2021) both derived in the spatially integrated manner.

However, the β value can spatially vary over the galaxy, as shown in the nearby galaxies (Galametz et al. 2012; Smith et al. 2012). We investigate the effect on the derived dust temperature and mass profiles if the β value changes across the galaxy from 1 to 2.5. Adopting a small β of 1 increases the overall temperature by ∼20 K, while a large β of 2.5 decreases the overall temperature by ∼2.5 K. This means a radially decreasing β from 2.5 to 1.0 can erase the temperature gradient of 20 K seen in Fig. 5. However, the scenario is unlikely because it contradicts the generally observed anticorrelation of the temperature and β in nearby galaxies (warmer/colder dust has lower/higher β: e.g. Galametz et al. 2012; Smith et al. 2012). The expected radial increase in β would increase rather than decrease the temperature gradient. The variation of β also would not significantly change the mass profile. The induced changes in the dust mass profile are at most 40 per cent, which is smaller than the uncertainty of dust mass due to the dust opacity coefficient κ₀ (at least a factor of 2, see e.g. Clark et al. 2019).

As a final remark on the analysis in this section, note that the measured flux of an individual pixel represents the integrated value weighted by the beam, so the derived temperature and optical depth need to be considered as the luminosity-weighted spatially averaged solutions over the region subtended by the beam (∼1.92 kpc²). It is expected that such spatial smoothing by the beam acts to flatten any intrinsic distribution.

As we will see in the next section, there is a central unresolved component and an extended component in the dust continuum images, and the central unresolved component coincides with the high-temperature peak position. If the central unresolved component is due to the warm dust heated by the AGN, the derived high temperature is just due to the result of forcibly fitting two dust components with different temperatures with a single greybody function. This contamination is potentially significant out to the beam FWHM in radius. Similarly, we may overestimate the SFR in this area if we naively convert from our derived spatially resolved TIR luminosity to SFR, ignoring the possibility that the warm dust can contaminate our measurement. Therefore, it is important to decompose the image into the unresolved and the extended component and remove the effect of the unresolved component to correctly estimate the SFR distribution in the galaxy. We do exactly this in the next section.

4 AGN-HOST GALAXY DECOMPOSITION

4.1 Image decomposition of Band 7 and Band 9 continuum images

In the previous section, we demonstrated that the dust temperature shows a steep increase up to ∼57 K, which is much higher than the typical dust temperature assumed for quasar hosts (47 K; Beelen et al. 2006) and is presumably due to the heated dust by the central AGN. The observed dust continuum images (Fig. 2) show a strongly peaked structure in the centre as well as a discy structure with spiral- and bar-like features in the outer part, suggesting the presence of a distinct central unresolved component and a resolved extended component. The former is likely to be associated with the warm dust heated by the central AGN and the latter with the cold dust heated by star formation in the host galaxy.

For these reasons, we decompose the continuum image into an unresolved AGN component and an extended host galaxy component using the observed [C ii] moment 0 map as a template tracing the dust distribution heated by the star formation in the galaxy. This approach has several advantages: (i) there is only one parameter (total flux), (ii) the [C ii] emission captures disc substructures like arms and a bar that are also seen in the dust continuum map, and that are not easily representable by an analytic model, and (iii) there is independent evidence that [C ii] is a good overall neutral gas tracer (e.g. Herrera-Camus et al. 2018a, b); we expect this gas to form stars at a rate roughly proportional to the density (Krumholz, Dekel & McKee 2011), and there is good evidence that [C ii] directly correlates with star formation (De Looze et al. 2011). [C ii] in this galaxy also has rotating disc kinematics and a nearly exponential profile typical of discs (Sersic index n ∼ 0.9; Tsukui & Iguchi 2021). Given these considerations, our decomposition model includes the following components: (i) a component proportional to the PSF of the observation I_PSF to reproduce the unresolved emission, (ii) a component with the same spatial distribution as the observed [C ii] intensity map representing a discy, star-forming component that exhibits spiral and bar-like structure, I_disc, and (iii) a circular Gaussian component, convolved by the PSF of the observation, I_Gauss, which we include in case the central compact region includes emission that is extended rather than truly point-like (e.g. a nuclear starburst), and thus is partially resolved by the ALMA beam. The model with these three components is described as

$$\begin{eqnarray} I_{\mathrm{model}}&= & I_{\mathrm{PSF}}(cx,cy, f_{\nu ,\mathrm{point}}) \\&& +\,I_{\mathrm{Gauss}}(cx,cy, \sigma _{\mathrm{Gauss}}, f_{\nu ,\mathrm{Gauss}}) \\&& +\,I_{\mathrm{disc}}(f_{\nu ,\mathrm{disc}}), \end{eqnarray}$$

(4)

where the model has only six parameters: the central coordinates of the PSF and Gaussian components (cx, cy); the intrinsic size of the Gaussian σ_Gauss (standard deviation); and the total flux of each component f_ν,point, f_ν,Gauss, and f_ν,disc. We create the spatial template for the disc component I_disc from the [C ii] moment map as follows: First, we convolve the [C ii] cube to the same resolution as the observed continuum image we are fitting, and then we make a [C ii] moment map by the masked moment method (Dame 2011; convolution with a kernel of 2 times the beam size and 2 times the spectral channel width, then masking with a threshold set to the RMS of the original cube) to reduce the bias introduced by the presence of noise in the [C ii] moment map. This bias pushes the best-fitting model towards lower f_ν,disc values because larger f_ν,disc values amplify the noise present in the [C ii] moment map, giving an additional penalty in χ²; convolution and masking suppress this noise in the spatial template and therefore suppress the bias.

We fit I_model to the observed Band 7 and Band 9 continuum images with simple χ² minimization. For both the Band 7 image |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and the Band 9 image |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠, evaluation of the Bayesian information criterion (BIC) indicates that this three-component fit is preferred over any alternatives omitting one of the components. Full details are provided in Appendix [see Fig. A3 and Fig. A4 for BIC and chi-square values, see Fig. A5 and Fig. A6 for fitting residuals]. Figs 8 and 9 show the original image (a: see Fig. 2 for the same images in log scale), the best-fitting PSF model (b), the Gaussian model (c) and the [C ii] model (d), and the best-fitting residual (e) for Band 7 and Band 9 continuum images, respectively.

Figure 8.

Decomposition of the rest frame 161 μm continuum image: (a) the original image on a linear scale; (b) the best-fitting point source component; (c) the best-fitting circular Gaussian component; (d) the best-fitting [C ii] component; and (e) the absolute residual |data − model|. The intensity scale for all maps is in units of the RMS noise σ of the continuum image (a). We show the size of the synthesized beam (FWHM) in the lower-left corner of each panel.

Figure 9.

Decomposition of the rest frame 90 μm continuum image. Same as in Fig. 8 except showing Band 9 results.

In Table 2, we show the best-fitting parameters we derive. The uncertainties on these parameters are computed by Monte Carlo resampling with noise which has realistic spatial correlations. Our procedure to make this estimate is to: (i) measure the autocorrelation function (ACF) of the noise map in the primary beam uncorrected Band 7 and Band 9 continuum images, (ii) randomly generate noise maps with the same correlation properties characterized by the noise ACF using the ESSENCE package (Tsukui et al. 2022, 2023), and (iii) repeat the fitting procedure 300 times with different random realizations of the noise. We then take the 68 per cent confidence interval on the resulting distribution of parameters as our confidence interval in Table 2.

Table 2.

The best-fitting result of the image decomposition. The x and y offset of the centre of the point source and Gaussian is denoted relative to the image centre coordinate (RA, Dec) = (204.514231 deg, −4.543055 deg), which corresponds to the highest peak pixel of the temperature. The offset of the quasar position from the centre is RA = 14.5 mas and Dec = 4.5 mas, consistent with that of the point source component within the ALMA’s positional uncertainty. Model total flux f_ν,model is the sum of each component flux f_ν,point + f_ν,Gauss + f_ν,disc. Measured total flux f_ν,total are measured from the images, same as listed in Table A1.

Parameters	Band 7 (⁠|$\lambda _{\mathrm{rest}}\ =\ 161\ {\mu\rm{m}}$|⁠)	Band 9 (⁠|$\lambda _{\mathrm{rest}}\ =\ 90\ {\mu\rm{m}}$|⁠)
x offset of the centre (RA), cx (mas)	19.6 ± 0.5	9.6 ± 1.6
y offset of the centre (Dec), cy (mas)	13.6 ± 0.6	1.3 ± 1.4
Point source flux, fν,point (mJy)	3.09|$^{+0.09}_{-0.08}$|	10.87|$^{+0.85}_{-1.18}$|
Gaussian flux fν,Gauss (mJy)	7.4|$^{+0.18}_{-0.14}$|	16.83|$^{+1.08}_{-1.06}$|
[C ii] template flux fν,disc (mJy)	10.26|$^{+0.25}_{-0.36}$|	27.62|$^{+2.32}_{-2.19}$|
Gaussian size σGauss (pc)	787|$^{+18}_{-17}$|	777|$^{+66}_{-75}$|
Model total flux fν,model (mJy)	20.76 ± 0.36	55.32 ± 2.69
Derived parameters
Measured total flux fν,total (mJy)	20.54 ± 0.27	52.92 ± 2.13

Our fits pass several consistency checks. First, the measured total flux f_ν,total from the observed images and the model total flux f_ν,model = f_ν,point + f_ν,Gauss + f_ν,disc agree within the statistical uncertainty for both images. Second, the centre position of the point source and Gaussian components coincides with the location of the optical quasar position and the highest peak in the continuum images and temperature map within the uncertainties. Third, the sizes of the additional Gaussian components that we derive independently by fitting the two continuum images are consistent within the uncertainty. We speculate that this component is likely needed to compensate for the central decrease in the [C ii]/FIR mentioned in Section 3. We do see some correlated residuals in the residual map of |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| (Fig. 8), which shows a positive residual at the location of the southern arm but not the northern arm. This indicates that dust continuum emission is enhanced relative to [C ii] emission in the southern arm but not in the northern arm. The residual of |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| also shows a positive residual in a similar position, but without a clear spiral structure (Fig. 9). This result may be interpreted as the two arms having different interstellar medium (ISM) properties (opacity, temperature, metallicity, etc.). However, despite this feature, our simple model with only six parameters describes well the overall morphology of the continuum emission without any other significant structures in the residual.

Fig. 10 shows the surface brightness profile of the Band 7 and Band 9 continuum images and their best-fitting decomposed components. The contribution of the point sources to the integrated flux over the image is small, 15 per cent and 21 per cent in rest frame 161 and 90 μm, respectively (Table 2). However, the point source dominates the flux in the central pixel (flux in the central resolution element, or beam area), contributing 49.8|$_{-1.3}^{+1.5}$| per cent and 56.7|$^{-6.1}_{+4.5}$| per cent in |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠, respectively. The best-fitting model profile is slightly higher than the dust continuum profile at the outskirts of the galaxy (∼0.6 arcsec; Fig. 10, left-hand panel). This indicates the [C ii] model component is more extended than the FIR emission, even after the central compact point source and the compact Gaussian component are removed. This is consistent with the presence of extended [C ii] halo structures extending farther than the underlying FIR emission (out to ∼10 kpc) in other star-forming galaxies at a similar redshift (e.g. Fujimoto et al. 2019, 2020; however, for a contrary view see Novak et al. 2020, who find no evidence for halo structures in quasar host galaxies at z > 6 using stacking analysis). The relatively short (1 h) Band 7 observation relying on the brightness of the source is not sensitive enough to probe the faint structure out to 10 kpc in BRI1335–0417. Note that in what follows, we use only the point source fluxes f_ν,point and the measured total fluxes f_ν,total from the images, not the fluxes of the individual Gaussian f_ν,Gauss or [C ii] disc components f_ν,disc, or the total model fluxes f_ν,model. Therefore, the slight model-data offset in the surface brightness profile at the outer part of the galaxy does not affect our subsequent results.

Figure 10.

The surface brightness profiles for rest frame 161 μm (left-hand panel) and 90 μm (right-hand panel) continuum images. Blue points with an error bar show the azimuthally averaged surface brightness profile of the continuum image, while lines show the best-fitting decomposed model, including the point source (dotted grey line), Gaussian (dot–dashed grey line), and [C ii] (dashed grey line) components, and their sum (black solid line). These profiles are derived using Photutils code (Bradley et al. 2019) as follows. We first fit elliptical isophotes to the dust continuum images, where the four parameters describing each ellipse (2D centre position, PA, and ellipticity) are left free. We then measure azimuthally averaged surface brightness along each of the best-fitting ellipses. The error bars on the data points include uncertainties in the ellipse fitting due to spatially correlated noise in the images.

4.2 SED fitting using the image decomposition result

In the previous section, we found that an unresolved component is required to reproduce the dust continuum images, and we decompose the images into a sum of this component and the rest of the galaxy. In this subsection, we model the SEDs of the unresolved component (AGN-heated dust), with flux f_ν,point, and the rest of the galaxy, with flux f_ν,total − f_ν,point (colder dust heated by the star formation in the disc). We use the fluxes of the point-like components |$f_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}, \mathrm{point}}$| and |$f_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}, \mathrm{point}}$| to constrain the AGN-heated warm dust contribution to the FIR part of SED, along with total flux measurements at various wavelength bands listed in Table A1. Note that we do not use the fluxes of the rest of the galaxy f_ν,total − f_ν,point as independent data points, since they are already constrained implicitly by the total flux measurements f_ν,total and point source flux measurements f_ν,point. We consider several possible approaches to fitting the resolved and unresolved components described below.

4.2.1 Two greybody components SED modelling

We first model the FIR part of the SED with two greybody functions in the optically thin limit (equation 1, τ ≪ 1), one for the AGN-heated dust and the other for the cold dust heated by star formation. Although detailed AGN dust torus models exist to describe the warm dust component heated by AGN (Stalevski et al. 2016; Hönig & Kishimoto 2017), the spatial extent and geometry of tori are highly uncertain for high redshift quasars. This motivates us to consider a simpler model where we fit the AGN-dust component with a greybody spectrum in the optically thin limit but with a free dust emissivity index β. This model has six free parameters (T_dust, M_dust, and β_dust for each component). The sum of the two greybody functions (AGN-heated dust and colder dust heated by star formation) is fitted to the spatially integrated fluxes (from rest frame 36–472 |${\mu\rm{m}}$| bands in Table A1), and the AGN-heated dust greybody is fitted to point fluxes, |$f_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}, \mathrm{point}}$| and |$f_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}, \mathrm{point}}$| (Table 2). The fitting is done simultaneously to constrain both the AGN-heated dust component and the colder dust component.

We compute the distribution of their posteriors for a uniform prior distribution using the emcee package (Foreman-Mackey et al. 2013). Fig. 11 shows the SED decomposed into the AGN-heated dust component and the host galaxy. Fig. A7 shows the posterior distribution of the model parameters, indicating that all model parameters are well-constrained. From this fit, we find that the SFR of the host galaxy is 1.5|$^{+0.3}_{-0.2}\times 10^{3}\ \mathrm{M}_\odot$| yr⁻¹ based on the FIR luminosity of the cold dust component, which is more than 3 times smaller than the previous estimate 5040 ± 1300 M_⊙ yr⁻¹ (Wagg et al. 2014). This is due to the high-AGN contribution, 63 per cent, to the FIR luminosity, which is neglected in the previous study. The AGN-heated dust component has a temperature |$T_{\mathrm{dust}}=95.60^{+37.33}_{-25.81}$| K, which is consistent with the high temperature reported for the centre of high redshift quasars found in earlier high-resolution observations (∼200 pc, Shao et al. 2022; Walter et al. 2022). However, from this fit, we also find that the temperature of the cold dust component is |$T_{\mathrm{dust}}=29.3^{+4.5}_{-3.6}$| K, which is lower than the temperature 29–47 K with a median value of ∼38 K at the outer part of the galaxy (Fig. 5). This is expected due to the optically thin assumption we adopted in the SED modelling, which yields systematically lower temperatures than fits that allow the dust to have finite optical depth (Cortzen et al. 2020).

Figure 11.

The FIR SED of BRI1335−0417 (black points, see Table A1) with the best-fitting model (grey) composed of two greybody spectra: a warm dust component heated by the AGN (red) and a cold dust component associated with the host galaxy (green). Both components are assumed to be optically thin. In addition to the total photometric points, we used point source fluxes, |$f_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}, \mathrm{point}}$| and |$f_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}, \mathrm{point}}$| (red points) measured by the image decomposition to constrain the flux contributed by the AGN-heated warm dust component (Table 2).

Improving this cour model by introducing an additional free parameter for the dust optical depth τ_ν at |$161\ {\mu\rm{m}}$|⁠, which we denote 𝜏161 μm, for the cold dust component. We fit the model again with the same procedure for the above optically thin case, but with a prior constraint that the cold dust temperature cannot exceed the warm dust temperature. Fig. A8 shows the posterior distribution of the model parameters, which are summarised in Table 3.⁵ Including a finite optical depth increases our best-fitting cold dust temperature to |$52.6_{-11.0}^{+10.3}$| K. This confirms that relaxing the optically thin assumption provides a dust temperature consistent with the spatially resolved result, while at the same time not significantly changing the overall functional shape of two greybody functions or the parameters derived from them such as the SFR and the AGN luminosity. The updated SED fit is shown in Fig. A9. The derived range of optical depth 𝜏161 μm from 0 to 0.0146 is much smaller than the spatially resolved result with the median value of 0.1556. The optical depth derived by the spatially integrated SED assuming single temperature T_dust may have much less physical meaning than the optical depth derived in the spatially resolved image (Fig. 4b) because the galaxy shows a significant variation of non-linear parameters such as temperature and optical depth over the galaxy. This may illustrate that the optical depth τ_ν measurement becomes inaccurate if the dust temperature and opacity structure are smaller than the integration aperture.

Table 3.

The derived parameter with SED modelling. ^aThe stardust code also returns the average radiation field intensity <U > = 66.0 ± 6.4 (=L_dust/M_dust/125), which is a proxy of the luminosity-weighted dust temperature T_dust. The temperature is converted by using the <U > = (T_dust/18.9 K)^6.04 following Magdis et al. (2017).

Parameters	Unit	Two greybodies	Two greybodies	Stardust
Cold dust opacity assumption		Thin	Thick	Thin
log10Mdust (warm dust)	(M⊙)	|$7.98_{-0.13}^{+0.13}$|	|$8.02_{-0.13}^{+0.13}$|
Tdust (warm dust)	(K)	|$96.0_{-25.9}^{+38.9}$|	|$87.1_{-18.3}^{+34.1}$|
βdust (warm dust)		|$0.85_{-0.49}^{+0.56}$|	|$0.91_{-0.49}^{+0.49}$|
log10Mdust (cold dust)	(M⊙)	|$9.48_{-0.15}^{+0.14}$|	|$9.00_{-0.13}^{+0.16}$|	9.39 ± 0.03
Tdust (cold dust)	(K)	|$29.3_{-3.6}^{+4.5}$|	|$52.6_{-11.0}^{+10.3}$|	37.8 ± 0.6a
βdust (cold dust)		|$2.75_{-0.40}^{+0.46}$|	|$2.56_{-0.30}^{+0.37}$|
|$\tau _{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|			0.0072 (<0.0146)
Derived parameters
Total LIR	(1013 L⊙)	|$3.9_{-0.8}^{+1.3}$|	|$3.8_{-0.6}^{+0.9}$|	5.4 ± 0.1
AGN LIR	(1013 L⊙)	|$2.4_{-0.9}^{+1.3}$|	|$2.0_{-0.7}^{+1.1}$|	3.6 ± 0.1
Host galaxy LIR	1013 L⊙	|$1.5_{-0.2}^{+0.3}$|	|$1.7_{-0.4}^{+0.5}$|	1.77 ± 0.02
SFR	(103 M⊙ yr−1)	|$1.5_{-0.2}^{+0.3}$|	|$1.7_{-0.4}^{+0.5}$|	1.77 ± 0.02
AGN fraction in LIR	–	|$0.62_{-0.13}^{+0.11}$|	|$0.53_{-0.15}^{+0.14}$|	0.67 ± 0.01
χ2 (DOF)	–	0.98 (4)	0.66 (3)	284.19 (20)

From this fit, we also obtain dust emissivity indices |$\beta _{\mathrm{dust}}=0.91^{+0.56}_{-0.49}$| and |$\beta _{\mathrm{dust}}=2.56^{+0.37}_{-0.30}$| for the warm and cold components, respectively. Our finding of a shallow index β ∼ 1 for the warm component and a steeper index β ∼ 2.5 for the cool component is consistent with the large range of 1 ≲ β ≲ 2.5 suggested by the observations and models of infrared to submillimetre dust emission and laboratory experiments (e.g. laboratory experiments: Agladze et al. 1996; Mennella et al. 1998, models: Pollack et al. 1994, and observations: Dupac et al. 2003; Galametz et al. 2012). In agreement with our results, these studies suggest that β_dust anticorrelates with the dust temperature T_dust; β ∼ 1 for small grains, primarily radiating in MIR bands, and β ∼ 2 for large grains, which can achieve lower equilibrium temperature and thus radiate their energy in FIR bands (da Cunha, Charlot & Elbaz 2008).

4.2.2 Stardust panchromatic SED modelling

Our third fitting approach is to model the panchromatic SED from the rest UV to FIR bands using the more realistic galaxy SED fitting code stardust (Kokorev et al. 2021), which assumes three components; a QSO template dominating in the rest-frame UV to optical bands (Shen 2016), an AGN-heated dust component (Mullaney et al. 2011) dominating in the rest-frame MIR bands, and a cold dust component of the host galaxy dominating in the FIR bands (Draine & Li 2007). The original stardust code fits the combination of the above three templates to the spatially integrated photometric data points given as input. In addition to the spatially integrated data, we modified the code to accept the point source fluxes as independent constraints to be fitted with the AGN-heated dust component. Using the modified code, we simultaneously fit (i) the sum of three templates to the spatially integrated data and (ii) the AGN-heated dust component to point source fluxes, |$f_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}, \mathrm{point}}$| and |$f_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}, \mathrm{point}}$|⁠.⁶ Fig. 12 shows the best-fitting SED decomposed into the QSO component (AGN-heated dust) and the host galaxy component (cold dust associated with the star-forming host galaxy). The SFR and AGN contribution to L_TIR we derive from this fit are consistent with those we obtain from the two-component greybody fit (Table 3). Our stardust calculation also finds that a QSO template without dust extinction (A_V = 0) fits the data better than a template including extinction, which may indicate that the QSO has already cleared dust at least from the line of sight between us and the central accretion disc. We also consider a fit from stardust using a stellar spectrum template instead of a quasar template for the UV-to-optical part of SED, but this results in a χ² value that is larger by 536, with 16 degrees of freedom. This difference in χ² is highly significant and strongly favours the quasar template. This result is also supported by the HST optical image (Fig. 1) consistent with the PSF. Moreover, the stellar template fit produces an estimated stellar mass of ∼10¹² M_⊙, which easily exceeds the dynamical mass. This confirms that the UV to optical part of the SED is dominated by AGN emission, supporting our interpretation that the compact dust component identified at the position of the highest temperature peak is predominantly heated by AGN.

Figure 12.

The ultraviolet (UV) to radio SED of BRI1335−0417 (black points, see Table A1) with the best-fitting stardust model (Kokorev et al. 2021) composed of a quasar template (blue solid line: UV to optical), AGN-heated hot and warm dust template [red solid line: near-IR (NIR) to MIR], and cold dust associated with the host galaxy (green solid line: FIR). In addition to the total photometric points, we used point source fluxes, |$f_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}, \mathrm{point}}$| and |$f_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}, \mathrm{point}}$| (red points) measured by the image decomposition to constrain the flux contributed by the AGN-heated warm dust component (Table 2). The dotted vertical line denotes the Ly α line wavelength.

For completeness, we also perform stardust SED fitting without using the image decomposition results, |$f_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}, \mathrm{point}}$| and |$f_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}, \mathrm{point}}$|⁠, to constrain the AGN-heated dust contribution to the FIR SED. In this case, stardust provides a best-fitting model with a AGN subtracted SFR of 5260 ± 31 M_⊙, much higher than that derived when including the spatially resolved constraints; moreover, the χ² value of this fit (χ² = 282, DOF = 18) is no worse than that derived including the spatially resolved data. We show the fit in Fig. A10, and visual inspection confirms that the fit to the global SED is no worse than the one shown in Fig. 12 derived using the spatially decomposed constraints. This demonstrates that spatially resolved information at rest frame ∼90–|$161\ {\mu\rm{m}}$| band is required to constrain and remove the AGN-heated dust component in order to estimate SFR accurately; it is not possible to make this correction with unresolved photometric data alone.

4.2.3 Conclusion of the SED fitting results

The most important conclusion from the discussion above is that all three methods of SED modelling that include the spatially resolved constraints provide consistent estimates of the AGN contribution to L_TIR, and therefore for the host galaxy SFR (Table 3). As final values for the remainder of this paper, we adopt the best-fitting model parameters derived from the two greybodies fit with a free optical thickness |$\tau _{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| for the cold component (fourth column of Table 3); this is a conservative choice since it includes the widest range of possible model uncertainties. Our fiducial estimates are therefore |$T_{\mathrm{dust}}=87.1_{-18.3}^{+34.1}$|⁠, |$\beta _{\mathrm{dust}}=0.91_{-0.49}^{+0.49}$|⁠, and |$T_{\mathrm{dust}}=52.6_{-11.0}^{+10.3}$| K, |$\beta _{\mathrm{dust}}=2.56_{-0.30}^{+0.37}$| for the warm and cold components, respectively. The former is consistent with expectations for small dust primarily radiating in the MIR (Cunha et al. 2008), while the latter is consistent with the typical temperatures and dust opacity indices of high-redshift starburst galaxies (Magnelli et al. 2012), bolstering our interpretation of these two components as primarily AGN-heated and primarily star formation-heated.

4.3 AGN-subtracted dust properties

Armed with our successful decomposition of the emission into a point-like AGN component and an extended star-formation component, we are now also in a position to construct a map of the spatially resolved dust properties in the galaxy with the AGN contamination removed. To this end, we start from the rest frame 161 and 90 μm continuum images, and in each pixel, we subtract our best-fitting model estimate of the contribution from the point-like AGN-heated dust. We then repeat the fitting procedure presented in Section 3 on the point source-subtracted images. This procedure yields a map of properties for the cold dust component alone, rather than the unknown admixture of cold and hot components we obtain in Section 3. Figs 13a–c show the temperature map, dust optical depth, and surface density of the SFR, respectively, derived in this manner. We again use Monte Carlo simulations to propagate the uncertainties, including not only the noise in the image and the uncertainty in the fiducial dust emissivity index β = 2.14 ± 0.17 as in Section 3, but also our uncertainties on the properties of the point source component (see Table 2).

Figure 13.

The best-fitting dust temperature (a), optical depth at Band 7 continuum λ_rest = 161 μm or dust mass surface density (b), and surface density of the integrated TIR luminosity over 8–1000 μm, Σ_TIR (L_⊙ kpc⁻²) (c). These are derived by fitting a greybody function at each pixel of the dust continuum images |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| after subtracting the point source (AGN-heated dust) contribution. The colour scales and contour levels are the same as in Fig. 4 in order to aid in comparison. The size of the synthesized beam (FWHM) is shown in the lower-left corner of panels (A–C).

Fig. 14 shows the 1D dust temperature profile of the cold dust and the star-formation density estimated after the warm dust component is removed. The temperature still shows a slight increase towards the centre, but does not show the steep temperature gradient seen in the results without the subtraction. The flat dust temperature indicates that the steep gradient in Fig. 5 is entirely due to the warm dust heated by the AGN, suggested by the compact nature and the two-component dust SED fitting (Figs 11 and 12), providing further confidence for the successful subtraction of the warm dust component heated by the central AGN.

Figure 14.

Radial profile of the temperature of BRI1335–0417, derived after subtracting the point source component (AGN-heated dust) shown in blue points and blue shade extracted from Fig. 13a in the same way as Fig. 5. The grey shade shows the radial profile derived before subtracting the point source component (same as the grey shade Fig. 5).

Fig. 15 shows a comparison of the SFR surface density Σ_SFR estimated from the FIR luminosity derived with and without removing the warm dust component heated by AGN (Figs 4c and 13c). Σ_SFR is derived by (Σ_TIR/L_⊙ kpc⁻²) × 10⁻¹⁰ M_⊙ kpc⁻² yr⁻¹ assuming the Chabrier (2003) initial mass function, where Σ_TIR is measured for the cold dust component after subtracting the unresolved warm dust component from the images (Fig. 13c). We see that if one does not remove the component, the central SFR surface density can be overestimated by over a factor of 3 at currently available resolutions, and that the point source AGN-heated dust contaminates the derived dust physical properties out to ∼1 synthesized beam FWHM in radius (∼1.3 kpc). In future higher resolution observations, the contaminated region would presumably be smaller, but the contribution of the AGN-heated warm dust would be correspondingly larger in the central resolution elements, leading to an even larger overestimation of the central SFR surface density.

Figure 15.

Radial profile of SFR surface density Σ_SFR of BRI1335–0417, derived after subtracting the point source component (AGN-heated dust) shown in blue points and blue shade. The points and shades are extracted from Fig. 13c in the same way as Fig. 14 after converting the surface TIR luminosity density to SFR density (see text). The grey shade shows the radial profile of Σ_SFR derived without subtracting the point source component, derived by converting Σ_TIR as shown in Fig. 4c to Σ_SFR.

Fig. 16 compares the SFR surface density estimated from the point source-subtracted FIR luminosity with that estimated from the [C ii] luminosity using the calibration of De Looze et al. (2011). As expected from the central decrease in [C ii]/L_TIR, the [C ii] luminosity underpredicts the SFR surface density in the central region. In the outer parts of the galaxy, the predicted SFR density from the [C ii] line luminosity agrees well with that derived from L_TIR. As a consistency check, we also compute the total SFR from the resolved map by integrating over it, using SFRs derived from [C ii] in pixels where dust-based estimates are unavailable because the 68 per cent confident interval on T_dust is >10 K (white regions in the image). The resulting integrated SFR is ∼1476 M_⊙ yr⁻¹, consistent with the SED modelling result |$1700_{-400}^{+500}\ \mathrm{M}_\odot$| yr⁻¹.

Figure 16.

Radial profile of SFR surface density Σ_SFR of BRI1335–0417 derived from L_TIR shown in Fig. 13c, compared with the SFR surface density derived from [C ii] line luminosity |$L_{\mathrm{[C~\small {II}]}}$| using the calibration of De Looze et al. (2011).

5 DISCUSSION

The decomposition of BRI1335–0417’s emission into AGN- and star formation-powered components has implications both for the interpretation of the source itself and future studies of similar systems. We explore these in turn.

5.1 Implications for the dynamics and history of BRI1335–0417

Our analysis allows us to draw a number of conclusions about the structure formation of BRI1335–0417. First, our results suggest that the galaxy has an extended disc and a compact bulge both actively forming stars. Part of the evidence is dynamical: Tsukui & Iguchi (2021) conclude based on [C ii] kinematics that a compact mass <1.3 kpc in size resides in the centre of the galaxy. Our work here adds morphological evidence on top of this: Both the |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| images require a compact Gaussian component with an effective radius R_e of ∼0.9 kpc (calculated as 0.5 × FWHM from the sizes obtained in Table 2) in addition to a nearly exponential [C ii] disc template. The most natural explanation for the combined dynamical and morphological evidence is the presence of a compact, massive nuclear region that is actively star-forming. However, the fact that the rest-frame UV to optical part of the SED is best fit by a QSO template without dust extinction also suggests that the AGN of BRI1335–0417 has cleared dust along our line of sight to the central accretion disc (e.g. Fujimoto et al. 2022), so gas is actively being ejected, and the central region may be in the last throes of star formation. Thus, we conclude that BRI1335–0417 may be in the transition from a QSO phase to being a passive bulge-dominated galaxy (Hopkins et al. 2008a, b). With the compact star formation and the extended disc star formation, the galaxy may evolve into the passive bulged-disc system (Fudamoto et al. 2022) unless the galaxy is affected by external perturbation such as mergers.

Our results also explain the organized disc rotation in the galaxy. In Fig. 17, we show pixel-by-pixel measurements of Σ_gas versus Σ_SFR (i.e. the Kennicutt–Schmidt diagram) for BRI1335–0417. The gas mass is derived by converting the spatially resolved dust mass surface density in Fig. 4b to a total mass surface density assuming the same gas-to-dust ratio of 54.2 ± 9.3 estimated in Section 3. We find that most pixels are located in the starburst regime characterized by gas depletion times of 50–200 Myr (Genzel et al. 2010), much less than the ∼1 Gyr values typically found for main-sequence star-forming galaxies at redshift of 1–3 (Tacconi et al. 2013). If the starburst is triggered by a gas-rich major merger, which violently disturbs the gas kinematics, it may take at least one orbital period for the gas to settle into organized disc rotation. The 50–200 Myr depletion time is comparable to the ≈120 Myr orbital period at the disc effective radius derived from [C ii] kinematics (Tsukui & Iguchi 2021). This explains why the galaxy has had time to settle into relaxed, disc-like kinematics, at least in the centre where the orbital period is shortest; the outer disc, where the orbital time is longer, is also rotating, but shows significant morphological disturbance (Tsukui & Iguchi 2021), consistent with a system age near the upper end of the depletion time range.

Figure 17.

The SFR surface density Σ_SFR and the gas surface density Σ_gas relation (Kennicutt–Schmidt diagram). Left-hand panel: The black points indicate our spatially resolved measurements of BRI1335–0417, where three pixels are extracted per beam, compared with the star-forming main sequence at redshift of 1–3 (Tacconi et al. 2013) and starburst galaxies at redshift of 1–3 (Genzel et al. 2010). The red cross shows the previous estimate for BRI1335–0417 (Jones et al. 2016), refined by our new spatially resolved SFR estimate with the AGN-heated dust component removed. Right-hand panel: Same as the black points in the left, but colour-coded with the radial distance from the centre. The diagonal lines in both panels indicate linear relations between Σ_SFR and Σ_gas corresponding different gas depletion time-scale Σ_gas/Σ_SFR.

It is worth noting that our estimated gas depletion times are much larger than the previous estimates by Jones et al. (2016), who not only found a much higher SFR due to AGN contribution, but also adopted a size of 1.8 kpc² for the area of the star-forming disc based on the 44 GHz continuum, much smaller than the size revealed in our spatially resolved maps. Both factors contribute to an underestimate of the depletion time. Using this underestimate leads to the conclusion that the system must be less than one dynamical time old, which is hard to reconcile with the kinematics as noted in Section 1. Our increased depletion time resolves this discrepancy.

Finally, we note that Fig. 17 also shows that the central region of the galaxy has shorter gas depletion time on average than the disc (outer) part of the galaxy, while the outer part has a larger scatter in depletion time. The large scatter reflects the high-/low-temperature regions aligned with the minor/major axis (Fig. 7), and suggests that the SFR in part of the outer disc may be locally enhanced by cold gas/satellite accretion from a preferential direction, or by the AGN wind to the perpendicular direction to the disc as discussed in Section 3.6. This combined with the shorter depletion time in the centre suggests that the galaxy may quench star formation inside out to form a bulge-disc system (Tacchella et al. 2015; van Dokkum et al. 2015) due to the combination of the removal of the gas by the AGN wind, earlier consumption of the gas in the centre, and/or gas supply from ongoing gas/satellite accretion supporting the star formation in the disc.

5.2 Implications for future studies of high-redshift starbursts and QSO hosts

In BRI1335–0417, the warm dust component heated by AGN contributes only a small fraction of the total observed flux, ∼15 per cent and ∼21 per cent at λ_rest = 161 μm and λ_rest = 90 μm, respectively. However, ignoring this component leads to a factor of three overestimate of the total SFR (see Figs 3 and 11), because the warm dust component has high temperature, |$T_\mathrm{dust}=87.1^{+34.1}_{-18.3}$|⁠, and thus contributes a substantially larger fraction of the inferred total TIR luminosity (⁠|$L_{\mathrm{TIR}}\propto T_{\mathrm{dust}}^{4+\beta }$|⁠). The warm dust component also contributes a significant fraction of the flux within the central resolution element, leading to a factor of ∼4 overestimate of the local SFR density (see Fig. 15), and contaminates the SFR density estimate out to roughly the FWHM of the beam in radius. Thus the point source subtraction is crucial for accurate estimates of both the total SFR and its spatial distribution, even when the AGN contribution to the observed bands is relatively small.

Nor does higher resolution remove this need. In a higher-resolution observation, the region contaminated by the AGN is smaller, but the AGN contribution is correspondingly greater because the relative fraction of the extended cold dust component becomes smaller in the central beam. As an extreme example of this, Walter et al. (2022) reported 200 pc resolution observation towards the quasar J234833.34–305410.0 at redshift 6.9 and estimated the dust temperature of the central resolution element (radius of 110 pc) to be >132 K, and the SFR density to be 25 500 M_⊙ yr⁻¹ kpc⁻² with the assumption that the AGN is negligible. The effective radius of the resolution element 110 pc is much smaller than our study (r = 720 pc). As the authors discussed, the compact emission from the AGN-heated dust may dominate in the central resolution element. Our findings here strongly support that conjecture.

In this paper, we decompose the AGN and star-formation contributions using ALMA observations of the FIR part of SED, where the ratio of AGN to galactic luminosity is less extreme than at rest-frame UV to optical, where the BH accretion disc greatly outshines the host galaxy (Marshall et al. 2020). Recently Ding, Silverman & Onoue (2022) demonstrated that JWST data in rest-frame optical also enable the decomposition of the surface brightness distribution into the point source AGN component and extended host galaxy stellar component. The separation of AGN-host galaxies with JWST, coupled with decomposition in the FIR and SED modelling such as that we have demonstrated here, will provide measurements of the already-formed stellar masses in addition to the SFR. These complementary measurements will provide the information crucial to understanding how the galaxy–BH relation is set in the early universe.

6 CONCLUSION AND SUMMARY

We present ∼1 kpc ALMA imaging of the dust continuum, [C ii] emission, and CO(7-6) emission of a z = 4.4 quasar host galaxy, BRI1335−0417. Due to the unique brightness of the galaxy, the observations provide a number of resolution elements across the galaxy, allowing us to study the spatially resolved ISM properties of the host galaxy, and to disentangle the light coming from the unresolved, AGN-powered region from that produced by star formation in the surrounding galaxy. Our main findings and their implications are the following.

• Using the spatially resolved continuum images |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠, but without first separating the AGN and galaxy components, we constrain the shape of the greybody for individual pixels and derive the dust temperature, optical thickness, and dust mass surface density. The central resolution element is found to be optically thick τ_ν ∼ 0.6 at |$\lambda _{\mathrm{rest}}\ =\ 161\ {\mu\rm{m}}$| and τ_ν ∼ 1.3 at |$\lambda _{\mathrm{rest}}\ =\ 90\ {\mu\rm{m}}$|⁠. The temperature shows a steep radial gradient towards the centre, reaching 57.7 ± 0.4 K, which is higher than the typical temperatures of 47 and 40 K for quasar host galaxies (Beelen et al. 2006) and star-forming galaxies (Magnelli et al. 2012), respectively.

• The pixel-by-pixel temperature distribution image shows that high-temperature peaks are preferentially aligned with the disc minor axis. The anisotropic distribution coincides with the high-velocity dispersion region of [C ii], both of which show a conical shape aligned with the disc minor axis. With current data, we cannot conclude whether the feature is due to the AGN wind heating the dust (Saito et al. 2018) or to the presence of a star-forming region driven by anisotropic cold gas supply from gas/satellite accretion (Dekel et al. 2009).

• Image decomposition analysis reveals the presence of a point source in the two dust continuum images |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠, whose positions coincide with the highest temperature, the peak of the dust continuum image, and the optical quasar position. The point source contribution to the total flux is small, ∼15 per cent and ∼21 per cent for |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠, respectively. However, in the central resolution element, the contribution is much larger, 49.8|$_{-1.3}^{+1.5}$| per cent and 56.7|$^{-6.1}_{+4.5}$| per cent for the respective images to the flux.

• We model the FIR SED assuming that the point-source flux comes from a warm dust component heated by the AGN, and the remainder comes from cold dust component heated by the star-forming host galaxy. The decomposed fluxes constrain the AGN-heated warm dust contribution to the FIR part of the SED, which is usually assumed to be dominated by the cold dust heated by star formation. We estimate temperatures of |$T_\mathrm{dust}=87.1^{34.1}_{-18.3}$| K and |$T_{\mathrm{dust}}=52.6^{+10.3}_{-11.0}$| K for the warm and cold components, respectively. We estimate the SFR from the FIR luminosity of the cold component, finding a SFR of |$1700_{-400}^{+500}\ \mathrm{M}_\odot$| yr⁻¹. This is a factor of three less than the previously estimated value 5040 ± 1300 M_⊙ yr⁻¹ due to the high AGN fraction in the FIR luminosity |$53^{+14}_{-15}$| per cent.

• The SED fit suggests that there are two dust components with different temperatures in the central resolution element in the images. The central temperature of 57 ± 0.3 K we obtain without decomposing the AGN and star-formation components should be interpreted as the luminosity-averaged solution within the beam, which results from fitting these two dust components with a single greybody function. After removing the point source (warm dust component) component from the images, we remeasure the dust properties for the star formation-heated cold dust component only. This fit shows a nearly flat temperature profile with a slight increase towards the centre, suggesting that the steep temperature gradient found in the one-component fit is entirely due to the unresolved warm dust component. The single-component fit overestimates the central SFR density by a factor of ∼4 as a result.

• Our estimates of SFR surface density Σ_SFR after AGN subtraction and gas surface density Σ_gas for individual pixels show a roughly linear sequence with a gas depletion time of 50–200 Myr. This places the galaxy in the starburst regime, clearly separated from main-sequence galaxies at z ∼ 1–3 with gas depletion times of ∼1 Gyr, but makes it a typical starburst rather than an extreme outlier.

Through this study, we demonstrate a method to constrain the SED shape of an unresolved warm dust component heated by AGN and an extended cold dust component in the host galaxy by combining spatially resolved information in two images |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$| and |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$| with integrated UV to FIR SED analysis. This method provides the total spatially resolved SFR, with contamination from warm dust heated by the AGN removed allowing us to study how quasar activity affects stellar mass build-up. Our method and recently demonstrated AGN-host galaxy decomposition with JWST (Ding et al. 2022) provide complementary measurements crucial to understanding how the galaxy–BH relation is set in the early universe.

ACKNOWLEDGEMENTS

TT is grateful for the helpful discussions with Satoru Iguchi, Takuya Hashimoto, Kentaro Nagamine, and Yuichi Matsuda. This research was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. MRK acknowledges support from the Australian Research Council through Laureate Fellowship FL220100020. Data analysis was carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00394.S, and #2018.1.01103.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC 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 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(s) GO 8572. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This paper makes use of the Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This work is based in part on archival data obtained with the Spitzer space telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by an award issued by JPL/Caltech. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

DATA AVAILABILITY

The ALMA data we use in this work are publicly available at: https://almascience.nrao.edu/aq/. The HST data we use in this work can be obtained from: https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html. The photometric data from Spitzer and Herschel we use in this work are publicly available at doi: 10.26131/IRSA3 and Marton et al. (2017), respectively.

Footnotes

1

Starburst galaxies are commonly referred to galaxies which exhibit the SFR three or four times above the tight SFR—stellar mass relation which majority of star-forming galaxies form (e.g. Elbaz et al. 2011; Rodighiero et al. 2011; Schreiber et al. 2015).

2

Strictly speaking, only two bands are not sufficient to constrain the three physical parameters of the greybody function, M_dust, temperature T_dust, and opacity index β_dust. However, tight constraints can be obtained by assuming the typical opacity index β_dust or choosing two bands near the peak of the spectrum which are sensitive to the dust temperature T_dust and insensitive to the β_dust.

3

We used the same casa version as used in the quality assurance at the ALMA Observatory, which is 5.1.1–5 for Band 7, and 5.4.0–70 for Band 4 and Band 9 data.

4

ALMA Technical Handbook: https://almascience.nrao.edu/documents-and-tools/cycle9/alma-technical-handbook

5

Note that, since the posterior distribution for τ_161μm remains significantly above zero all the way down to τ_161μm = 0, for this parameter report only a best-fitting value and an upper limit τ_161μm < 0.0146 (84th percentile of the posterior distribution), not a lower limit.

6

The modified stardust code will be available at the GitHub repository: https://github.com/takafumi291, upon publication.

APPENDIX A: SUPPLEMENTAL FIGURES

Figure A1.

The radial distribution of HST STIS/50CCD image (black line) and the PSF (blue line: Medallon, Welty & et al. 2023). The profile is measured by fitting a Moffat function to the image (Fig. 1). The grey-shaded region is the 1σ uncertainty of the measured profile due to statistical noise. This comparison demonstrates that the optical emission in the STIS/50CCD image is consistent with a point source within the uncertainty.

Figure A2.

Herschel SPIRE images at 250 |${\mu\rm{m}}$| (a), 350 |${\mu\rm{m}}$| (b), and 500 |${\mu\rm{m}}$| (c) bands for BRI1335−0417. The blue cross shows the source position of BRI1335−0417. The blue ellipse at the left corner of each panel is the size of the PSF (FWHM). At 250 |${\mu\rm{m}}$| with the highest resolution, the galaxy is separated from the nearby bright source at East. At 350 and 500 |${\mu\rm{m}}$| bands, the galaxy is not separated from and can be contaminated by the neighbouring bright source. Therefore, we did not include the fluxes of 350 and 500 |${\mu\rm{m}}$| for our SED modelling.

Figure A3.

The reduced chi-square values (left-hand panel) and BIC (right-hand panel) for all combinations of model components used to fit Band 7 continuum image, |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|⁠.

Figure A4.

Same as Fig. A3 for Band 9 continuum image, |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠.

Figure A5.

The residual of observed data minus best-fitting model for all combinations of model components used to fit the Band 7 continuum image, |$I_{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|⁠.

Figure A6.

The same as Fig. A5 for the Band 9 continuum image, |$I_{\lambda _{\mathrm{rest}}\ =\ 90\ \mu \mathrm{m}}$|⁠.

Figure A7.

Posterior probability distributions of the dust mass M_dust, dust temperature T_dust, and dust emissivity β_dust of two greybody functions fitted to FIR part of the SED, sampled using emcee code (Foreman-Mackey et al. 2013). The colour in the sampled distribution indicates the relative log-likelihood of the sample; blue shows the least likely to white the most, while black shows the even less likely points with Δχ² > 9, corresponding to 3σ confidence interval in χ² statistics. The derived parameters are summarised in Table 3.

Figure A8.

Same as Fig. A7 for two greybody functions with free optical depth at rest frame 161 |${\mu\rm{m}}$|⁠, τ_{161 μm} for the cold dust component. We used prior on the temperature of the warm dust component does not exceed the cold dust component. The derived parameters are summarised in Table 3.

Figure A9.

Same as Fig. 11, but for the fit relaxing the optically thin assumption with the free optical depth at rest frame 161 |${\mu\rm{m}}$|⁠, |$\tau _{\lambda _{\mathrm{rest}}\ =\ 161\ \mu \mathrm{m}}$|⁠. Optically thin case is included in the range of the confidence interval of the model fitting. The overall shape of the model does not change and thus provides consistent parameters such as SFR, and AGN fraction with the optically thin case (see Table 3).

Figure A10.

The UV to radio SED of BRI1335−0417 (black points, see Table A1) with the best-fitting stardust model (Kokorev et al. 2021) without using point source fluxes (red points) at rest frame 161 and 90 μm as AGN-heated dust component. The best-fitting model returns the SFR estimate of 5260 ± 31 M_⊙ yr⁻¹ much higher than 1770 ± 20 M_⊙ yr⁻¹ estimated using point source fluxes at rest frame 161 and 90 μm bands as emission from AGN-heated dust.

Table A1.

The data points used in this paper. From the left column, the name of the instrument, central wavelength in microns, measured flux, flux uncertainty, and sources. SEIP denotes Spitzer Enhanced Imaging Products Source List (doi: 10.26131/IRSA3) and HPPSC denotes the Herschel/PACS Point Source Catalogue (Marton et al. 2017). For ALMA photometric points, the absolute flux uncertainty (10 per cent of the flux) is added in the quadrature with the statistical uncertainty.

Instrument	Wavelength (micron)	Flux (mJy)	Flux uncertainty (mJy)	Sources
Pan-STARRS1-g	0.477	0.0099	3.100e−04	(Chambers et al. 2016)
Pan-STARRS1-r	0.612	0.053	1.300e−03	(Chambers et al. 2016)
Pan-STARRS1-i	0.747	0.102	1.000e−03	(Chambers et al. 2016)
Pan-STARRS1-z	0.865	0.118	1.000e−03	(Chambers et al. 2016)
Pan-STARRS1-y	0.960	0.123	3.000e−03	(Chambers et al. 2016)
VISTA-Y	1.018	0.143	4.000e−03	(Lasker et al. 2021)
2MASSJ	1.239	0.133	5.000e−03	(Lasker et al. 2021)
2MASSH	1.649	0.141	8.000e−03	(Lasker et al. 2021)
2MASSKs	2.163	0.135	8.000e−03	(Lasker et al. 2021)
WISE W1	3.350	0.140	7.000e−03	(Lasker et al. 2021)
Spitzer IRAC3.6	3.600	0.141	3.759e−04	SEIP
Spitzer IRAC4.5	4.500	0.114	4.042e−04	SEIP
WISE W2	4.600	0.116	0.011	(Lasker et al. 2021)
Spitzer IRAC5.8	5.800	0.151	1.487e−03	SEIP
Spitzer IRAC8.0	8.000	0.280	2.782e−03	SEIP
Spitzer MIPS24	24.00	1.484	0.215	SEIP
Herschel PACS 70	70	6.768	3.606	HPPSC
Herschel PACS 160	160	20.976	7.838	HPPSC
Herschel SPIRE 250	250	37.6	7.900	HPPSC
ALMA	484	52.9	5.7	Our work
ALMA	869	20.5	2.1	Our work
ALMA	1110	9.0	0.9	(Lu et al. 2018)
IRAM interferometer	1350	5.6	1.1	(Guilloteau et al. 1997)
ALMA	2012	1.17	0.12	(Lu et al. 2018)
ALMA	2080	1.00	0.10	Our work
VLA	61 685	0.076	0.011	(Carilli, Menten & Yun 1999)
VLA	203 940	0.220	0.043	(Carilli et al. 1999)