# The interstellar medium distribution, gas kinematics, and system dynamics of the far-infrared luminous quasar SDSS J2310+1855 at z = 6.0

Yali Shao, Ran Wang, Axel Weiss, Jeff Wagg, Chris L Carilli, Michael A

Strauss, Fabian Walter, Pierre Cox, Xiaohui Fan, Karl M Menten, et al.

► To cite this version:

Yali Shao, Ran Wang, Axel Weiss, Jeff Wagg, Chris L Carilli, et al.. The interstellar medium distribution, gas kinematics, and system dynamics of the far-infrared luminous quasar SDSS J2310+1855 at z = 6.0. Astronomy & Astrophysics - A&A, 2022, 668, pp.A121. 10.1051/0004-6361/202244610. hal-04265432

HAL Id: hal-04265432

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

Submitted on 30 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# The interstellar medium distribution, gas kinematics, and system dynamics of the far-infrared luminous quasar SDSS J2310+1855 at z = 6.0

Yali Shao<sup>1,2</sup>, Ran Wang<sup>2</sup>, Axel Weiss<sup>1</sup>, Jeff Wagg<sup>3</sup>, Chris L. Carilli<sup>4</sup>, Michael A. Strauss<sup>5</sup>, Fabian Walter<sup>6,4</sup>, Pierre Cox<sup>7</sup>, Xiaohui Fan<sup>8</sup>, Karl M. Menten<sup>1</sup>, Desika Narayanan<sup>9</sup>, Dominik Riechers<sup>10</sup>, Frank Bertoldi<sup>11</sup>, Alain Omont<sup>7</sup>, and Linhua Jiang<sup>2</sup>

1. Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
e-mail: yshao@mpifr-bonn.mpg.de

2. Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, PR China

3. SKA Observatory, Lower Withington Macclesfield, Cheshire SK11 9FT, UK

4. National Radio Astronomy Observatory, Socorro, NM 87801-0387, USA

5. Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

6. Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

7. Sorbonne Université, CNRS UMR 7095, Institut d'Astrophysique de Paris, 98bis bvd Arago, 75014 Paris, France

8. Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

9. Department of Astronomy, University of Florida, 211 Bryant Space Science Center, Gainesville, FL 32611, USA

10. I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany

11. Argelander-Institut für Astronomie, University at Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

Received 27 July 2022 / Accepted 19 October 2022

## ABSTRACT

We present Atacama Large Millimeter/submillimeter Array (ALMA) sub-kiloparsec- to kiloparsec-scale resolution observations of the [C II], CO (9–8), and OH$^+$ (1$_{1}$–0$_{1}$) lines along with their dust continuum emission toward the far-infrared (FIR) luminous quasar SDSS J231038.88+185519.7 at z = 6.0031, to study the interstellar medium distribution, the gas kinematics, and the quasar-host system dynamics. We decompose the intensity maps of the [C II] and CO (9–8) lines and the dust continuum with two-dimensional elliptical Sérsic models. The [C II] brightness follows a flat distribution with a Sérsic index of 0.59. The CO (9–8) line and the dust continuum can be fit with an unresolved nuclear component and an extended Sérsic component with a Sérsic index of ~1, which may correspond to the emission from an active galactic nucleus dusty molecular torus and a quasar host galaxy, respectively. The different [C II] spatial distribution may be due to the effect of the high dust opacity, which increases the FIR background radiation on the [C II] line, especially in the galaxy center, significantly suppressing the [C II] emission profile. The dust temperature drops with distance from the center. The effective radius of the dust continuum is smaller than that of the line emission and the dust mass surface density, but is consistent with that of the star formation rate surface density. This may indicate that the dust emission is a less robust tracer of the dust and gas distribution but is a decent tracer of the obscured star formation activity. The OH$^+$ (1$_{1}$–0$_{1}$) line shows a P-Cygni profile with an absorption at ~400 km s$^{-1}$, which may indicate an outflow with a neutral gas mass of (6.2 ± 1.2) × 10$^8$ M$_{\odot}$ along the line of sight. We employed a three-dimensional tilted ring model to fit the [C II] and CO (9–8) data cubes. The two lines are both rotation dominated and trace identical disk geometries and gas motions. This suggest that the [C II] and CO (9–8) gas are coplanar and corotating in this quasar host galaxy. The consistent circular velocities measured with [C II] and CO (9–8) lines indicate that these two lines trace a similar gravitational potential. We decompose the circular rotation curve measured from the kinematic model fit to the [C II] line into four matter components (black hole, stars, gas, and dark matter). The quasar-starburst system is dominated by baryonic matter inside the central few kiloparsecs. We constrain the black hole mass to be 2.97$^{+0.51}_{-0.77}$ × 10$^9$ M$_{\odot}$; this is the first time that the dynamical mass of a black hole has been measured at z ~ 6. This mass is consistent with that determined using the scaling relations from quasar emission lines. A massive stellar component (on the order of 10$^9$ M$_{\odot}$) may have already existed when the Universe was only ~0.93 Gyr old. The relations between the black hole mass and the baryonic mass of this quasar indicate that the central supermassive black hole may have formed before its host galaxy.

**Key words.** galaxies: high-redshift – quasars: general – submillimeter: galaxies

## 1. Introduction

Almost 400 quasars at $z \ge 5.7$ have been discovered in the past ~20 yr, mostly from optical wide-field multiband surveys (e.g., Fan et al. 2000; Matsuoka et al. 2022). These quasars provide a unique opportunity to study a number of key issues, for example the formation of young luminous quasars, the evolving impact of the central black holes on the host galaxies, and the typical

interstellar medium (ISM) conditions in the quasar host galaxies, during the epoch at which the intergalactic medium was being reionized by the first luminous sources.

Toward some of these $z \sim 6$ quasars, bright millimeter dust continuum emission has been detected at millijansky or sub-millijansky levels (e.g., Bertoldi et al. 2003a; Wang et al. 2007, 2008, 2011a; Omont et al. 2013; Venemans et al. 2018; Li et al. 2020c) using the IRAM facilities, the *James Clerk Maxwell*
Telescope (JCMT), and the Atacama Large Millimeter/submillimeter Array (ALMA), indicating far-infrared (FIR) luminosities of $10^{11}$ to $10^{13}$ L$_{\odot}$ and star formation at rates of hundreds to thousands $M_{\odot}$yr$^{-1}$ in these young quasar host galaxies. A few of these $z \sim 6$ FIR luminous quasars have been detected in multi-$J$ CO transitions from $J=2-1$ to $17-16$ (e.g., Bertoldi et al. 2003b; Walter et al. 2003, 2004; Carilli et al. 2007; Riechers et al. 2009; Wang et al. 2010, 2011a,b; Gallerani et al. 2014; Shao et al. 2019; Li et al. 2020b; Pensabene et al. 2021) using the IRAM facilities and ALMA, indicating abundant gas reservoirs on the order of $10^{10}$ M$_{\odot}$ and highly excited molecular gas in the quasar host galaxies. Most of these high-$z$ FIR-luminous quasars have been detected in [C II] 158 µm fine-structure line emission down to arcsecond- and subarcsecond-scale resolution using the IRAM facilities and ALMA (e.g., Maiolino et al. 2005; Walter et al. 2009, 2022; Wang et al. 2013; Shao et al. 2017; Decarli et al. 2018; Izumi et al. 2018; Venemans et al. 2020; Yue et al. 2021; Meyer et al. 2022), suggesting that high-luminosity quasar host galaxies have lower dynamical masses than local galaxies with similar black hole masses, although the ratios of the black hole mass to the host galaxy dynamical mass of most of the low-luminosity quasars are consistent with the local value. In addition, high-resolution CO, [C II], and ultraviolet (UV) continuum observations reveal close companions for some $z \sim 6$ quasars (e.g., Wang et al. 2011b, 2019; McGreer et al. 2014; Decarli et al. 2017; Miller et al. 2020; Izumi et al. 2021; Pensabene et al. 2021). These studies suggest an overall scenario of the coevolution of quasars and their host galaxies: the central supermassive black holes (SMBHs) may grow rapidly through major galaxy mergers and the accretion of large amounts of gas, which may trigger high rates of star formation in the quasar host galaxies and thus influence the relationship between the black hole mass and the dynamical mass of these quasar-host systems.

However, only a few of these studies are based on high-resolution (i.e., sub-kiloparsec scale) and high sensitivity observations in the millimeter (e.g., Yue et al. 2021; Walter et al. 2022). The quasar host galaxies at $z \sim 6$ are at most a few beam sizes across. Thus, the inferred ISM properties and the dynamical information of the quasar-host systems are mostly based on unresolved or slightly resolved spatially integrated flux stacking along the velocity or frequency axes (i.e., the intensity maps) and on the distribution of spatially integrated flux for each velocity interval as a function of apparent radial Doppler velocity (i.e., the integrated spectra). High-resolution (i.e., a few hundred parsecs) ALMA observations of the ISM in these quasar host galaxies at $z \sim 6$ can provide insight into these issues, for example the spatial distribution of different ISM tracers in addition to the star formation activity, the gas kinematics, and the overall dynamics probed by different ISM tracers in the SMBH-host system at the earliest epochs.

Measuring the sizes and morphologies of galaxies in the early Universe is critical for understanding the initial stage of galaxy formation and evolution. The long-wavelength FIR dust continuum traces regions of young and compact star formation that are severely obscured by dust at optical wavelengths. In both observed galaxies (e.g., Riechers et al. 2013, 2014; Ikarashi et al. 2015; Hodge et al. 2016; Wang et al. 2019; Tadaki et al. 2020) and simulated galaxies (e.g., Cochrane et al. 2019; Popping et al. 2022), the distribution of the dust continuum emission generally is more compact than both the cold gas and the dust mass, but is more extended (more compact) than the stellar component at $z \le 3$ ($z \ge 4$). The [C II] fine-structure

transition at 157.74 µm is one of the primary coolants of the star-forming ISM in galaxies (e.g., Stacey et al. 1991), and it traces both the neutral atomic and ionized gas. The different rotational transitions of CO trace cold and warm gas that will form stars in the future. Their sizes and morphologies are strong probes of ISM properties. Thus, comparing the size and morphology between the dust emission and the [C II] and CO lines will provide insights into the evolutionary process of the ISM in which star formation takes place.

The kinematics of the ISM can yield important information on the dynamical structure of the quasar host, as well as the distribution and fraction of baryonic and dark matter in the system (e.g., van de Hulst et al. 1957; Rubin 1983; de Blok et al. 2008). The rotation velocities of different phases of the ISM frequently show a complex picture and differ from source to source. The rotation velocities measured from the molecular components (i.e., CO) and atomic components (i.e., H I) are generally in agreement (e.g., Wong et al. 2004; Frank et al. 2016). Übler et al. (2018) find that the kinematics of CO and H α are in good agreement at a $z = 1.4$ star-forming galaxy. However, Simon et al. (2005) find that the H α line in NGC 4605 shows systematically slower rotation than the CO. de Blok et al. (2016) find that the velocity of the [C II] line is systematically larger than that of the CO or H I lines in a few nearby galaxies, which they attribute to systematics in the data reduction and the low velocity resolution of the [C II] data. It is still not clear why discrepancies between different kinematic tracers are observed in some galaxies.

In the nearby Universe, the stellar bulge can be observed directly through optical and near-infrared (NIR) imaging, and its dynamics can be probed via investigations of the rotation curve (e.g., Begeman et al. 1991; de Blok & Bosma 2002; Sofue et al. 2009; Gao et al. 2018). At high redshifts, the stellar bulge may already exist long before the peak of cosmic star formation, as demonstrated by the [C II] gas kinematics of $4 \le z \le 5$ dusty star-forming galaxies and quasars (e.g., Rizzo et al. 2020, 2021; Tsukui & Iguchi 2021). However, at higher redshifts, into the reionization epoch, it remains unknown whether the stellar bulge is already present.

In this paper we report on ALMA sub-kiloparsec- to kiloparsec-resolution observations of the [C II], CO (9–8), and OH$^+$ ($1_{150\AA} - 0_1$) lines and their underlying dust continuum toward the FIR-luminous quasar SDSS J231038.88+185519.7 (hereafter J2310+1855) at $z = 6.0031$, to study the ISM distribution, the gas kinematics, and the quasar-host system dynamics. This quasar is first reported by Wang et al. (2013) with 250 GHz dust continuum and CO (6–5) line observations using the IRAM facilities and [C II] line observations using the ALMA. Jiang et al. (2016) are credited with its discovery; they discovered this optically bright quasar with $m_{1450\AA} = 19.30$ mag in the Sloan Digital Sky Survey (SDSS) imaging data. The black hole mass has been measured to be $(4.17 \pm 1.02) \times 10^9$ M$_{\odot}$ and $(3.92 \pm 0.48) \times 10^9$ M$_{\odot}$ from the GEMINI/GNIRS spectra of Mg II and C IV lines, respectively (Jiang et al. 2016). A smaller black hole mass of $(1.8 \pm 0.5) \times 10^9$ M$_{\odot}$ based on the C IV emission line detected in the X-shooter/VLT spectrum is reported by Feruglio et al. (2018). All of these estimates of the black hole mass are based on the local scaling relations (e.g., Shen 2013). However, it is still under debate if the local relationship is suitable at high redshift.

J2310+1855 has a well-measured dust content and a detailed rest-frame NIR-to-FIR spectral energy distribution (SED) observed by SDSS, WISE, Herschel, IRAM, and ALMA. (Shao et al. 2019). These observations reveal a dust temperature
of ~40 K and a star formation rate (SFR) of ~2000 $M_{\odot}$ yr$^{-1}$ under the optically thin approximation. However, the dust and star formation distribution within the host are unknown, and the origins of different dust components have not yet been identified. The multi-$J$ CO transitions (from $J = 2-1$ to $13-12$) have been detected using the Karl G. Jansky Very Large Array (VLA), the IRAM facilities, and ALMA (Wang et al. 2013; Feruglio et al. 2018; Shao et al. 2019; Carniani et al. 2019; Li et al. 2020b; Riechers et al., in prep.). A detailed CO spectral line energy distribution analysis reveals that, in addition to the far-UV radiation from young and massive stars, another gas heating mechanism (e.g., X-ray radiation and/or shocks) may be needed to explain the observed CO luminosities (Carniani et al. 2019; Li et al. 2020b). The low excitation water para-H$_2$O ($20_2-11_1$) and para-H$_2$O ($32_2-31_3$), $OH^+$ ($1_1-0_1$) lines, and lines from ionized gas such as [O I] 146 µm, [O III] 88 µm, and [N II] 122 µm have been detected with ALMA toward J2310+1855, and a solar-level metallicity is proposed based on the [O III]/[N II] ratio (Hashimoto et al. 2019b; Li et al. 2020a; Tripodi et al. 2022). We should note that [N II] 122 µm is an upper [N II] line (i.e., the transition of $^3P_2-^3P_1$), so there are potential excitation effects. In addition, [O III] 88 µm and [N II] 122 µm lines have significantly different ionization potentials (~14 versus ~35 eV), so they do not necessarily trace the same parts of the H II regions. The bright [C II] line has been detected with ALMA at low angular resolution (>0.7; Wang et al. 2013; Feruglio et al. 2018). A velocity gradient is obvious, which may indicate that the disk is dominated by rotating gas. The ~4 kpc resolution [C II] data reveal a dynamical mass of $9.6 \times 10^{10} M_{\odot}$ with an approximate estimate of the inclination angle (46°, determined from the ratio between the minor and major axis), suggesting a $M_{BH}/M_{bulge}$ value that is higher than the local value (Wang et al. 2013). Tripodi et al. (2022) used ~1 kpc resolution [C II] data to find a dynamical mass of $5.2 \times 10^{10} M_{\odot}$ within a 1.7 kpc region, based on a kinematic modeling on the data cube with a rotating disk inclination angle of 25°. However, the limited spatial resolution introduces large uncertainties in the determination of the gas kinematics, making it difficult to perform a dynamical decomposition of the quasar-host system. In this paper we use ALMA high-resolution observations of the [C II] and CO(9-8) lines (~0.6 and ~1 kpc, respectively) and their underlying dust continuum in this quasar to explore its dynamics in detail.

The outline of this paper is as follows. In Sect. 2 we describe our ALMA observations and the data reduction. In Sect. 3 we present the measurements of the [C II], CO(9-8), and $OH^+$ ($1_1-0_1$) lines and the underlying dust continuum. In Sect. 4 we fit the line and continuum intensity maps to 2D Sérsic functions, constrain the dust properties, apply a 3D tilted ring model to the [C II] and CO(9-8) data cubes, and decompose the circular rotation curve measured from the high-resolution [C II] line into multiple components. In Sect. 5 we discuss the spatial distribution and extent of the ISM, the surface density of the gas and the star formation, the ionized and molecular gas kinematics, the gas outflow, and the dynamics of the quasar-host system. In Sect. 6 we summarize our results. Finally, in Appendices A-C, we present the channel maps of [C II] and CO(9-8) lines and describe the 2D Sérsic function and the tilted ring model. Throughout the paper we adopt a $\Lambda$CDM cosmology with $H_0 = 67.8$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.3089$, and $\Omega_\Lambda = 0.6911$ (Planck Collaboration XIII 2016). Under this cosmological assumption and at $z = 6.0031$, 1" on the sky corresponds to a physical size of 5.84 kpc, the luminosity distance is $D_L = 59.0$ Gpc, and the age was 0.9297 Gyr since the Big Bang.

## 2. ALMA observations and data reduction

We conducted ALMA band-6 observations of the [C II] line ($v_{\text{rest}} = 1900.5369$ GHz; redshifted to $v_{\text{obs}} = 271.3851$ GHz), along with band-4 observations of the CO(9-8) line ($v_{\text{rest}} = 1036.9124$ GHz; redshifted to $v_{\text{obs}} = 148.0648$ GHz) and the $OH^+$ ($1_1-0_1$) line ($v_{\text{rest}} = 1033.0582$ GHz; redshifted to $v_{\text{obs}} = 147.5144$ GHz; hereafter $OH^+$), toward J2310+1855 at $z = 6.0031$ from 2019 August 08 to 20 (PI: Yali Shao; Project code: 2018.1.00597.S). We used 40-47 12-m antennas in the C43-7 configuration with a maximum projected baseline of 3.6 km for observations with both bands. We centered one of the 2 GHz spectral windows (four in total) on the redshifted [C II]/CO(9-8) line frequency, and used the rest of the spectral windows to observe the dust continuum. The total observing times are 1.60 and 1.44 h, resulting in on-source integration times of 0.66 and 0.80 h for the [C II] and CO(9-8)/$OH^+$ observations, respectively. We established the flux density scale using scans of the standard ALMA calibrator J2253+1608. The flux calibration uncertainties are ~5-10% for ALMA Cycle 6 bands-4 and -6 observations (Warmels et al. 2018). We checked the phase and water vapor by observing nearby calibrators of J2316+1618 and J2307+1450. The data were calibrated using the ALMA standard pipeline using the Common Astronomy Software Application (CASA<fnref n="1" />). The original channel width is 15.625 MHz, corresponding to ~17 and 32 km s$^{-1}$ for the band-6 and -4 observations, respectively, which can sample the intrinsic line widths well. The underlying dust continuum emission was subtracted in the uv-plane for both data sets.

In order to improve the uv coverage and thus the final image and spectrum sensitivity, we also included archival data from projects – 2011.0.00206.S (PI: Ran Wang), 2015.1.00997.S (PI: Roberto Maiolino), and 2015.1.01265.S (PI: Ran Wang), which observed the [C II] line in the 12-m array and 12-m + 7-m arrays, and the CO(9-8)/$OH^+$ ($1_1-0_1$) line in the 12-m array, toward J2310+1855, respectively. We only consider the frequency range overlapping with our science goals. As we used ALMA Cycle 0 observations for which the data weights are not correct, before the final combination of multi-epoch data, we re-weighted the calibrated target data with the STATWT task in CASA, which attempts to assess the sensitivity per visibility and adjust the weights accordingly with line-free data. Finally, we made the line data cube from the combined calibrated data using the TCLEAN task in CASA with robust weighting (robust = 0.5) for the [C II] and CO(9-8) line in order to optimize the sensitivity per frequency bin and the resolution of the final maps, and natural weighting (robust = 2.0) for the $OH^+$ line in order to improve the sensitivity. For the continuum images, we used robust weighting (robust = 0.5), which allows us to compare the emission of the continuum with that of the covered lines ([C II] and CO(9-8)) at similar angular resolution. As the $OH^+$ line is very weak, we TCLEAN the data cube deeply with a threshold of $1\sigma$ using a 1" × 1" square centered on the quasar position. For the rest of the lines and continuum, we TCLEAN to a level of $3\sigma$. In addition, during TCLEAN for the 12-m and 7-m combined data, we used mosaic gridder mode, which can correctly image data with different antenna sizes. The synthesized beam size of the final [C II] and CO(9-8) images are 0.111 × 0.092 and 0.187 × 0.153, corresponding to 0.65 × 0.54 kpc and 1.09 × 0.89 kpc, respectively, at the quasar redshift. The noise levels in a 15.625 MHz channel are 0.17 and 0.10 mJy beam$^{-1}$ for the [C II] and CO(9-8) lines, respectively. The root mean

<footnote n="1">¹ https://casa.nrao.edu/</footnote>
<table><caption>**Table 1.** Measurements from ALMA observations.</caption><thead><tr><th></th><th>[C II]</th><th>CO (9–8)</th><th>OH<sup>+</sup></th></tr></thead><tbody><tr><td colspan="4"><strong>Line Intensity Map</strong></td></tr><tr><td>Weighting</td><td>0.5</td><td>0.5</td><td>2.0</td></tr><tr><td>Size<sub>beam</sub></td><td>0.111 × 0.092</td><td>0.187 × 0.153</td><td>0.243 × 0.198</td></tr><tr><td colspan="4"><strong>Line Spectrum</strong></td></tr><tr><td>z</td><td>6.0032 ± 0.0001 <sup>(a)</sup></td><td>6.0028 ± 0.0002 <sup>(a)</sup></td><td>—</td></tr><tr><td>FWHM (km s<sup>-1</sup>)</td><td>384 ± 5 <sup>(a)</sup></td><td>388 ± 16 <sup>(a)</sup></td><td>225 ± 117 <sup>(c)</sup>; 240 ± 74 <sup>(d)</sup>; 330 ± 198 <sup>(e)</sup></td></tr><tr><td>v<sub>S spectrum</sub> (Jy km s<sup>-1</sup>)</td><td>6.50 ± 0.17 <sup>(a)</sup>; 6.31 ± 0.18 <sup>(b)</sup></td><td>1.44 ± 0.11 <sup>(a)</sup>; 1.39 ± 0.11 <sup>(b)</sup></td><td>0.072 ± 0.049 <sup>(c)</sup>; 0.072 ± 0.023 <sup>(d)</sup>; -0.065 ± 0.044 <sup>(e)</sup></td></tr><tr><td>L<sub>spectrum</sub> (10<sup>9</sup> L<sub>☉</sub>)</td><td>6.39 ± 0.17 <sup>(a)</sup>; 6.21 ± 0.18 <sup>(b)</sup></td><td>0.77 ± 0.06 <sup>(a)</sup>; 0.74 ± 0.06 <sup>(b)</sup></td><td>0.038 ± 0.026 <sup>(c)</sup>; 0.038 ± 0.013 <sup>(d)</sup>; -0.035 ± 0.023 <sup>(e)</sup></td></tr><tr><td colspan="4"><strong>Continuum Map</strong></td></tr><tr><td>Size<sub>beam_con</sub></td><td>0.113 × 0.092</td><td>0.192 × 0.156</td><td>—</td></tr><tr><td>γ<sub>con</sub> (GHz)</td><td>262</td><td>147</td><td>—</td></tr></tbody></table>

**Notes.** The parameters for the line spectra are measured using a single and double Gaussian fit for [C II] and CO (9–8) lines, and a triple-Gaussian fit for the OH⁺ line. The adopted [C II], CO (9–8) and OH⁺ line spectra are shown in the right panels of Figs. 1–3. For the [C II] and CO (9–8) lines: (a) single Gaussian fit results; (b) double-Gaussian fit results. For the OH⁺ line: (c) blue part emission component; (d) red part emission component; (e) absorption component.

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

**Fig. 1.** ALMA observed [C II] line. Left panel: [C II] velocity-integrated map. The white plus sign is the Hubble Space Telescope (HST) quasar position (RA of 23:10:38.90; Dec of 18:55:19.85). The shape of the restoring beam with a FWHM size of 0.′′111 × 0.′′092 is plotted in the bottom-left corner. The contour levels are [-3, 3, 6, 9, 12, 15] × rms$_{[C II]}$, where rms$_{[C II]}$ = 0.035 Jy beam⁻¹ km s⁻¹. Middle panel: [C II] velocity map produced using pixels above 4σ in the channel maps. The black plus sign is the HST quasar position. A clear velocity gradient can be seen. Right panel: [C II] line spectrum (black histogram) overplotted with the best-fit Gaussian profiles. The spectrum is extracted from the 2σ region of the source emitting area in the [C II] intensity map. The solid blue line presents the best-fit single Gaussian model. The solid red line represents the best-fit double-Gaussian model, and the dashed and dotted red lines are for each of the double Gaussians. The kinematic local standard of rest (LSRK) velocity scale is relative to the [C II] redshift from our ALMA Cycle 0 observations (Wang et al. 2013). The spectral resolution is 15.625 MHz, corresponding to 17 km s⁻¹. The rms of the spectrum measured using the line-free spectrum is 0.65 mJy and is shown as a black bar. The [C II] spectral profile is asymmetric, with enhancement on the red side. The green bar indicates the velocity range of the OH⁺ absorption.

square (rms) noise in the continuum maps at 262 and 147 GHz are 0.02 and 0.01 mJy beam⁻¹, respectively.

### 3. Results

The [C II] and CO (9–8) lines and their underlying dust continuum are all spatially resolved. The intensity peak of OH⁺ (1₁–0₁) is detected at >5σ significance. We list the observational results in Table 1.

#### 3.1. The [C II] line

The velocity-integrated intensity map, the intensity-weighted velocity map and the spectrum of the [C II] line are presented in Fig. 1. A clear velocity gradient can be seen in the velocity map. In addition, the [C II] emission peak moves in a circular, counterclockwise path from 219 to -281 km s⁻¹ shown in the [C II] line channel maps (Fig. A.1). These are

the main characteristics of a disk with rotating gas. Similar rotating disks have been widely detected in the local and high-z Universe in both quasar hosts and galaxies without an active galactic nucleus (AGN; Wang et al. 2013; Lucero et al. 2015; Shao et al. 2017; Jones et al. 2017, 2021; Banerji et al. 2021; Alonso-Herrero et al. 2018; Bewketu Belete et al. 2021). We measured the line spectrum by integrating the intensity of each channel from the [C II] line data cube, including pixels determined in the line-emitting region above 2σ in the [C II] intensity map. The line spectrum with original spectral resolution (15.625 MHz) is shown in the right panel of Fig. 1, which reveals that this line has an asymmetric profile with an enhancement to positive velocities.

The single Gaussian fit to the [C II] spectrum shows that the line width and the source redshift are consistent with our previous Cycle 0 observations (393 ± 21 km s⁻¹ and 6.0031 ± 0.0002, respectively; Wang et al. 2013). The [C II] line flux calculated from the Gaussian fit to the line spectrum agrees with our previous ALMA observations at 0.′′7 resolution
<table><caption>**Table 2.** Parameters derived from best-fit of the image decomposition.</caption><thead><tr><th>Species</th><th>Model</th><th>Center<sub>x</sub><br>(hh:mm:ss)</th><th>Center<sub>y</sub><br>(dd:mm:ss)</th><th>n<sub>ISM</sub></th><th>r<sub>s_x</sub><br>(")</th><th>r<sub>s_y</sub><br>(")</th><th>θ<br>(°)</th><th>r<sub>e_x</sub><br>(kpc)</th><th>r<sub>e_y</sub><br>(kpc)</th><th>v<sub>S</sub><sub>v_point</sub><br>(Jy km s<sup>-1</sup>)</th><th>v<sub>S</sub><sub>v_ext</sub><br>(Jy km s<sup>-1</sup>)</th><th>S<sub>v_point</sub><br>(mJy)</th><th>S<sub>v_ext</sub><br>(mJy)</th><th>f<sub>point_all</sub><br>(%)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th><th>(13)</th><th>(14)</th><th>(15)</th></tr></thead><tbody><tr><td>[C II]</td><td>S</td><td>23:10:38.8994 ± 0.0001s</td><td>18:55:19.7925 ± 0"0011</td><td>0.59<sup>+0.01</sup><sub>-0.01</sub></td><td>0.236<sup>+0.003</sup><sub>-0.003</sub></td><td>0.165<sup>+0.002</sup><sub>-0.003</sub></td><td>43<sup>+1</sup><sub>-1</sub></td><td>1.262<sup>+0.012</sup><sub>-0.011</sub></td><td><b>0.883<sup>+0.008</sup><sub>-0.007</sub></b></td><td>—</td><td>6.85<sup>+0.05</sup><sub>-0.06</sub></td><td>—</td><td>—</td><td>—</td></tr><tr><td>CO (9–8)</td><td>S</td><td>23:10:38.9016 ± 0.0001s</td><td>18:55:19.7781 ± 0"0019</td><td>2.01<sup>+0.07</sup><sub>-0.06</sub></td><td>0.013<sup>+0.002</sup><sub>-0.002</sub></td><td>0.010<sup>+0.002</sup><sub>-0.001</sub></td><td>109<sup>+3</sup><sub>-4</sub></td><td><b>1.088<sup>+0.040</sup><sub>-0.035</sub></b></td><td>0.82<sup>+0.03</sup><sub>-0.03</sub></td><td>—</td><td>1.51<sup>+0.03</sup><sub>-0.03</sub></td><td>—</td><td>—</td><td>—</td></tr><tr><td>CO (9–8)</td><td>P+S</td><td>23:10:38.9016 ± 0.0001s</td><td>18:55:19.7787 ± 0"0019</td><td>1.21<sup>+0.09</sup><sub>-0.07</sub></td><td>0.087<sup>+0.012</sup><sub>-0.015</sub></td><td>0.066<sup>+0.010</sup><sub>-0.011</sub></td><td>107<sup>+4</sup><sub>-4</sub></td><td><b>1.259<sup>+0.054</sup><sub>-0.045</sub></b></td><td>0.949<sup>+0.040</sup><sub>-0.040</sub></td><td>0.16<sup>+0.02</sup><sub>-0.02</sub></td><td>1.31<sup>+0.04</sup><sub>-0.04</sub></td><td>—</td><td>—</td><td>11<sup>+1</sup><sub>-1</sub></td></tr><tr><td>[C II]con</td><td>S</td><td>23:10:38.9003 ± 0.0001s</td><td>18:55:19.8033 ± 0"0017</td><td>1.57<sup>+0.02</sup><sub>-0.02</sub></td><td>0.022<sup>+0.001</sup><sub>-0.001</sub></td><td>0.020<sup>+0.001</sup><sub>-0.001</sub></td><td>33<sup>+1</sup><sub>-1</sub></td><td>0.652<sup>+0.003</sup><sub>-0.003</sub></td><td><b>0.587<sup>+0.003</sup><sub>-0.003</sub></b></td><td>—</td><td>—</td><td>—</td><td>9.06<sup>+0.03</sup><sub>-0.03</sub></td><td>—</td></tr><tr><td>[C II]con</td><td>P+S</td><td>23:10:38.9003 ± 0.0001s</td><td>18:55:19.8032 ± 0"0016</td><td>1.16<sup>+0.02</sup><sub>-0.02</sub></td><td>0.053<sup>+0.002</sup><sub>-0.002</sub></td><td>0.047<sup>+0.002</sup><sub>-0.002</sub></td><td>33<sup>+1</sup><sub>-1</sub></td><td>0.687<sup>+0.004</sup><sub>-0.004</sub></td><td><b>0.621<sup>+0.004</sup><sub>-0.004</sub></b></td><td>—</td><td>—</td><td>0.51<sup>+0.02</sup><sub>-0.02</sub></td><td>8.31<sup>+0.04</sup><sub>-0.04</sub></td><td>6<sup>+1</sup><sub>-1</sub></td></tr><tr><td>CO (9–8) con</td><td>S</td><td>23:10:38.9013 ± 0.0001s</td><td>18:55:19.7740 ± 0"0005</td><td>1.18<sup>+0.03</sup><sub>-0.03</sub></td><td>0.039<sup>+0.003</sup><sub>-0.002</sub></td><td>0.034<sup>+0.002</sup><sub>-0.002</sub></td><td>79<sup>+3</sup><sub>-3</sub></td><td><b>0.522<sup>+0.006</sup><sub>-0.006</sub></b></td><td>0.452<sup>+0.005</sup><sub>-0.005</sub></td><td>—</td><td>—</td><td>—</td><td>1.52<sup>+0.01</sup><sub>-0.01</sub></td><td>—</td></tr><tr><td>CO (9–8) con</td><td>P+S</td><td>23:10:38.9013 ± 0.0001s</td><td>18:55:19.7741 ± 0"0005</td><td>1.17<sup>+0.05</sup><sub>-0.04</sub></td><td>0.042<sup>+0.003</sup><sub>-0.004</sub></td><td>0.037<sup>+0.003</sup><sub>-0.003</sub></td><td>80<sup>+2</sup><sub>-3</sub></td><td><b>0.558<sup>+0.012</sup><sub>-0.012</sub></b></td><td>0.484<sup>+0.011</sup><sub>-0.010</sub></td><td>—</td><td>—</td><td>0.08<sup>+0.02</sup><sub>-0.02</sub></td><td>1.45<sup>+0.02</sup><sub>-0.02</sub></td><td>5<sup>+2</sup><sub>-1</sub></td></tr></tbody></table>

**Notes.** Column 1: the name of the ISM tracer. Column 2: the best-fit model. “S” and “P” present the Sérsic and point component, respectively. Columns 3–4: source position – RA and Dec, respectively. Column 5: the Sérsic index of the ISM (the gas and the dust emission). Columns 6–7: the scale lengths along the major and minor axes, respectively. Column 8: the rotation angle defined in Appendix B. Columns 9–10: the half-light radii. The values in boldface are measurements near the kinematic major axes. Columns 11–12: the line flux for the point component and extended component, respectively. Columns 13–14: the dust continuum flux density for the point component and extended component, respectively. Column 15: the fraction of the nuclear component in the whole emission. Note that the errors of these parameters are all fitting-type errors. The complex correlated noise is not included. As claimed by the task reference of the CASA IMFIT, the correlated noise can bring 4% effect when 5σ detection of the emission for a 2D elliptical Gaussian fitting.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Similar to Fig. 1, but for the CO (9–8) line. The FWHM size of the restoring beam is 0".187 × 0".153. The contour levels are [-3, 3, 6, 9, 12, 15, 18, 21]×rms<sub>co</sub>, where rms<sub>co</sub> = 0.020 Jy beam<sup>-1</sup> km s<sup>-1</sup>. The spectral resolution is 15.625 MHz, corresponding to 32 km s<sup>-1</sup>. The rms of the spectrum is 0.31 mJy. The CO (9–8) spectral profile is also asymmetric, with enhancement on the red side.</figcaption></figure>

(8.83 ± 0.44 Jy km s<sup>-1</sup>; Wang et al. 2013) within the ~15% calibration uncertainty. Our double-Gaussian fit results are consistent with that of the single Gaussian fit. We also measured a consistent value by modeling the observed intensity map with the 2D elliptical Sérsic model as described in Sect. 4.1 and shown in Table 2.

## 3.2. The CO (9–8), OH<sup>+</sup> lines and the possible detection of para-H<sub>2</sub>O ground-state emission line

The velocity-integrated intensity map, the intensity-weighted velocity map and the spectrum of the CO (9–8) line are presented in Fig. 2. A similar asymmetric component at ~–50 km s<sup>-1</sup> is seen as in the [C II] line.

The Gaussian fit of the CO (9–8) line reveals a line width and source redshift that are consistent with the previous ALMA Cycle 3 observations (376 ± 18 km s<sup>-1</sup> and 6.0031 ± 0.0002, respectively; Li et al. 2020b). The CO (9–8) line flux calculated from the Gaussian fit on the line spectrum agrees with the previous ALMA observations at ~0".7 resolution (1.31 ± 0.06 Jy km s<sup>-1</sup>; Li et al. 2020b). We also obtained a consistent value by modeling the observed intensity map with the 2D elliptical Sérsic model shown in Table 2 (Sect. 4.1).

The OH<sup>+</sup> velocity-integrated intensity map and spectrum are shown in Fig. 3. The peak flux density is 0.066 ± 0.012 Jy beam<sup>-1</sup> km s<sup>-1</sup> (>5σ) in the intensity map averaged only over the emission component. There appears a double-horn emission profile and an absorption at ~–400 km s<sup>-1</sup>. We used a triple-Gaussian model (red line in Fig. 3) to fit the spectrum, and we present the flux and line parameters in Table 1. The P-Cygni profile (with an absorption at –384 ± 128 km s<sup>-1</sup> from the Gaussian fit) of the OH<sup>+</sup> spectrum may indicate an outflow.

Our new ALMA band-4 observations serendipitously partially cover the para-H<sub>2</sub>O (1<sub>11</sub>–0<sub>00</sub>) line at v<sub>rest</sub> = 1113.3430 GHz (redshifted to 158.9786 GHz), which is detected in emission in J2310+1855. The spectrum is shown in Fig. 4, where we overlaid a spectrum of the para-H<sub>2</sub>O (2<sub>02</sub>–1<sub>11</sub>) line from the ALMA archive. The para-H<sub>2</sub>O (1<sub>11</sub>–0<sub>00</sub>) line is usually detected as an absorption line in most galaxies (e.g., Weiß et al. 2010; Yang et al. 2013); however it shows up as a conspicuous but weak emission feature in a few galaxies (e.g., González-Alfonso et al. 2010; Spinoglio et al. 2012; Liu et al. 2017). The cold (T<sub>dust</sub> ~ 20–30 K) and widespread diffuse ISM component gives rise to the ground-state para-H<sub>2</sub>O line (Weiß et al. 2010; Liu et al. 2017). Liu et al. (2017) find the ratio of para-H<sub>2</sub>O (1<sub>11</sub>–0<sub>00</sub>) and para-H<sub>2</sub>O (2<sub>02</sub>–1<sub>11</sub>) to be <1 in
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** ALMA observed OH⁺ line. Left panel: OH⁺ emission line velocity-integrated map. The white plus sign is the HST quasar position. The shape of the restoring beam with a FWHM size of 0.243 × 0.198 is plotted in the bottom-left corner. The contour levels in black are [3, 4, 5] × rms<sub>OH⁺ em</sub>, where rms<sub>OH⁺ em</sub> = 0.012 Jy beam⁻¹ km s⁻¹. The peak flux density is 0.066 ± 0.012 Jy beam⁻¹ km s⁻¹. The overplotted dashed white contours are from the OH⁺ absorption line, with contour levels of [-3, -2] × rms<sub>OH⁺ ab</sub>, where rms<sub>OH⁺ ab</sub> = 0.008 Jy beam⁻¹ km s⁻¹. Right panel: OH⁺ line spectrum (black histogram) overplotted with the best-fit triple-Gaussian profile (red line). The dashed and dotted blue lines present the blueward and redward emission components, respectively, and the dash-dotted blue line represent the absorption component. The spectrum is extracted from the 2σ region of the source emitting area in the OH⁺ intensity map. The LSRK velocity scale is relative to the [C II] redshift from our ALMA Cycle 0 observations (Wang et al. 2013). The spectral resolution is 31.25 MHz, corresponding to 64 km s⁻¹. The rms of the spectrum is 0.05 mJy and is shown as a black bar. The green bar indicates the full width at zero intensity of the [C II] line, which is roughly consistent with that of the best-fit model for the OH⁺ emission line.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Spectra of the para-H₂O (1₁₁–0₀₀) line (black histogram) and the scaled para-H₂O (2₀₂–1₁₁) line (blue histogram). The LSRK velocity scale is relative to the [C II] redshift from our ALMA Cycle 0 observations (Wang et al. 2013). The spectral resolution for both lines is 15.625 MHz, corresponding to ~32 km s⁻¹. The observations may only cover the red tail of the para-H₂O (1₁₁–0₀₀) emission line.</figcaption></figure>

galaxies in which both lines are detected as emission lines. As the ground-state para-H₂O line is very close to the edge of the spectral window in our observations, we only detect it partially and may not cover the peak of its spectrum. Thus, we only simply report the detection of this line and will not discuss it further in what follows.

### 3.3. The dust continuum emission

The dust continuum distribution, overplotted with that of the [C II] and CO (9–8) lines, is shown in the left and middle panels of Fig. 5. We study the brightness distribution and the effective radius (half-light radius) for both lines and their underlying dust continuum emission in Sect. 4.1. The measured intensities are shown in Table 2. The *z* ~ 6 source is spatially resolved with our data, so we investigate the brightness distributions with 2D elliptical Gaussian functions. However, for the low-resolution data, the observed brightness distribution is dominated by the large beam, so the CASA IMFIT tool, which fits a 2D elliptical Gaussian, can be used effectively. The flux densities (listed in Table 2) of the [C II] and CO (9–8) underlying dust continua measured using the fits in Sect. 4.1 are consistent with the ALMA Cycle 0 data (8.91 ± 0.08 mJy; Wang et al. 2013) and ALMA Cycle 3 data (1.59 ± 0.04 mJy; Li et al. 2020b), respectively, measured using the CASA IMFIT tool including additional calibration uncertainties. As shown in Table 2, the effective radii of the [C II] and CO (9–8) lines are larger than that of their underlying dust continuum emission, as we discuss in Sect. 5.1.1.

## 4. Analysis

### 4.1. Image decomposition

In order to investigate the ISM distribution, we used a 2D elliptical Sérsic function (Eq. (B.3)) to reproduce the observed intensity maps of the [C II] and CO (9–8) lines and their underlying dust continuum. The 2D elliptical models represent the inclined brightness distribution. We first convolved the model with the restoring beam kernel, and then determined the best fit Sérsic model by comparing the convolved model and the
<figure><img src="image_6.png" /><figcaption>Fig. 5. Dust continuum images at different frequencies and their radial profile comparison. Left and right panels: [C II] and CO (9–8) underlying dust continuum maps at 262 and 147 GHz, respectively. The black plus sign is the quasar position from HST snapshot observations. The bottom-left ellipse in each panel shows the restoring beam with FWHM sizes of 0.113 × 0.092 and 0.192 × 0.156, respectively. The black contours are [−3, 3, 9, 27, 81] × rms<sub>con</sub> for both maps. The rms<sub>con</sub> at 262 and 147 GHz are 0.02 and 0.01 mJy beam<sup>−1</sup>, respectively. The overplotted white contours are from [C II] and CO (9–8) emission lines, and the contour levels are the same as those of the black contours in the left panels of Figs. 1 and 2. Right panel: comparison of the surface brightness of the continuum at two different frequencies. The blue and red symbols with error bars are for the scaled 262 and 147 GHz dust continuum, respectively. The error bar represents the deviation of the values of all pixels in each ring used to do the aperture photometry.</figcaption></figure>

data pixel-by-pixel using Markov chain Monte Carlo (MCMC) method with the emcee<fnref n="2" /> package (Foreman-Mackey et al. 2013). The best-fit models are shown in the left panels of Figs. 6 (for the [C II] line), 7 (for the CO (9–8) line), B.2 (for the dust continuum underlying the [C II] emission), and B.3 (for the dust continuum underlying the CO (9–8) emission). During the fitting, we considered uniform error maps as we only used the central regions of the observed intensity maps. Motivated by the possible existence of a torus surrounding the AGN in the AGN unification model and from observations of several molecular lines and dust continuum (e.g., Hönig 2019; Combes et al. 2019; García-Burillo et al. 2021), we also add an additional unresolved nuclear component (which is concentric with the extended component) to represent the dusty and molecular torus for the CO (9–8) line and the dust continuum emission. For the [C II] line, we do not include this component, as the observed [C II] emission is suppressed by the high dust opacity, as we discuss in Sect. 5.1.1. We note that in the following analysis and discussion on the spatially resolved gas and dust content in the quasar host galaxy, we remove the emission associated with the dusty and molecular AGN torus (i.e., subtract the unresolved nuclear component model from the observed image). To determine which model (Sérsic or point source+Sérsic) can better fit the line and/or the dust continuum data, we first measured the pixel-to-pixel rms within a 1'' × 1'' square at the center of the residual map (the measured values are shown in the captions of Figs. 7, B.2 and B.3), and then calculate the distribution of residual values inside the square (the top-right histogram of Figs. 7, B.2, and B.3). In addition, we compared the surface brightness profiles of both the modeled and observed intensity maps with projected rings (corrected by the inclination angle and position angle given in Table 3 determined from the kinematic models in Sect. 4.3).

The parameters of the best-fit image decomposition for each emission line and its associated continuum are listed in Table 2. For the [C II] line (shown in Fig. 6), the best-fit Sérsic index is 0.59, which is smaller than a typical exponential disk with a Sérsic index of 1. It is also smaller than that of the CO (9–8) line

and the dust continuum of our quasar J2310+1855. This discrepancy may be due to the high dust opacity, rather than the intrinsic distribution of the [C II] line, as we discuss in Sect. 5.1.1. For the CO (9–8) line, the best-fit Sérsic index is 2.01 in the case of a single Sérsic model. When we add a point component in the center, the Sérsic index of the extended component decreases to 1.21. The statistics on the residual maps for the one- and two-component scenarios appear identical. But the two-component model better matches the observed CO (9–8) within the central ~1 kpc region, as shown by the surface brightness difference distribution in Fig. 7. The [C II] and CO (9–8) underlying dust continua, similar to the CO (9–8) line, seem better described by a combination of a point component and an extended Sérsic component, as shown in Figs. B.2 and B.3. And the Sérsic index for the dust continuum emission is a little bit larger than 1. The dust continuum is more concentrated than the gas emission, as measured by the half-light radius. This is consistent with what is found in observational and theoretical studies of high-redshift galaxies (e.g., Strandet et al. 2017; Tadaki et al. 2018; Cochrane et al. 2019). We discuss the spatial distribution and extent of the ISM in Sect. 5.1. The central positions (see Cols. 3–4 in Table 2) of the gas and the dust continuum are identical.

## 4.2. Dust diagnostic

With the dust emission size from the image decomposition, and following Weiß et al. (2007) and Walter et al. (2022), we are able to constrain the dust temperature and the dust mass with a general gray-body formula instead of using the optically thin approximation as done in our previous work (Shao et al. 2019):

$$S_{\nu} = \Omega_{\text{app}}[B_{\nu}(T_{\text{dust}}) - B_{\nu}(T_{\text{CMB}})][1 - \exp(-\tau_{\nu})](1 + z)^{-3}, \quad (1)$$

and

$$\tau_{\nu} = \kappa_{0}(\nu/\nu_{\text{ref}})^{\beta} M_{\text{dust, app}}/(D_{A}^{2}\Omega_{\text{app}}), \quad (2)$$

where $S_{\nu}$ is the observed flux density for a target at redshift z. $B_{\nu}(T_{\text{dust}})$ and $B_{\nu}(T_{\text{CMB}})$ are the black-body functions with dust

<footnote n="2">² http://dfm.io/emcee/current/</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Image decomposition of the [C II] line. Left and middle panels: 2D elliptical Sérsic model and the residual between the observed and modeled [C II] intensity maps. The shape of the [C II] synthesized beam with a FWHM size of 0.111 × 0.092 is plotted in the bottom-left corner of each panel. We measured the residual rms to be 0.037 Jy beam⁻¹ km s⁻¹ inside the dashed gray square with a side length of 1". Top-right panel: [C II] luminosity surface brightness (black squares and red diamonds with error bars are measurements from the observed intensity map and the 2D elliptical Sérsic model, respectively) at different radii, measured using elliptical rings with the ring width along the major axis half (0.05) that of the [C II] clean beam size, the rotation angle equal to PA (=199°), and the ratio of semiminor and semimajor axis – *b*/*a* of cos(î) (î = 42°), where PA and î come from the [C II] line kinematic modeling. The error bar represents the standard deviation of the values of all pixels in each ring. The solid and dashed gray lines are three times the [C II] luminosity surface brightness limit measured in the emission-free region with an elliptical annulus that is the same with the ones used to measure the [C II] luminosity surface brightness but with a larger radius, and its corresponding rms. Bottom-right panel: luminosity surface brightness difference (red diamonds) between the measurements from the observed intensity map and the 2D elliptical Sérsic model. Note that the vertical scale is linear and in units of 10⁷ L⊙ kpc⁻².</figcaption></figure>

temperature $T_{\text{dust}}$ and the cosmic microwave background (CMB) temperature $T_{\text{CMB}}$, respectively. $\Omega_{\text{app}}$ is the apparent solid angle. In principle, $\Omega_{\text{app}}$ can be different for different wavelengths. Due to the dust continuum sizes remaining constant within ~20% at observed frame wavelengths from 500 µm to 2 mm for $z \sim 1-5$ main-sequence galaxies from simulations (Popping et al. 2022), and our findings for the similar dust sizes (within ~10%) at wavelengths of 158 and 290 µm in the rest frame, we do not consider changes in this quantity below. $\tau_{\nu}$ is the optical depth, which is a function of frequency $\nu$ (rest-frame) and dust mass surface density $M_{\text{dust, app}}/(D_A^2 \Omega_{\text{app}})$. $M_{\text{dust, app}}$ is the apparent dust mass. $D_A$ is the angular distance. $\beta$ is the emissivity index. We here adopted the absorption coefficient per unit dust mass $\kappa_0$ of 13.9 cm² g⁻¹ at the reference frequency $\nu_{\text{ref}}$ of 2141 GHz (Draine 2003; Walter et al. 2022).

Motivated by the AGN contribution to the strength of the dust continuum emission of the quasar-host system (e.g., Di Mascia et al. 2021), we decomposed the UV-to-FIR SED of J2310+1855 shown in Fig. 8. With the best-fitting UV/optical power law (red lines) and NIR to mid-infrared (MIR) CAT3D AGN torus model (brown lines; Hönig & Kishimoto 2017) in Shao et al. (2019), and the source size of twice the effective radius of the dust continuum underlying CO(9–8) listed in Table 2, we derived an average dust temperature $T_{\text{dust}}$ of 53$_{-4}^{+4}$ K, a total dust mass $M_{\text{dust, app}}$ of 1.27$_{-0.35}^{+0.62}$ × 10⁹ M☉ and an overall emissivity index $\beta$ of 1.90$_{-0.27}^{+0.33}$. Our dust temperature is smaller than the value of 71 ± 4 K derived by Tripodi et al. (2022), and our dust mass is ~3 times higher. This is due to the differences in the AGN dust torus model used, which has a significant contribution in the rest frame wavelength range ≲100 µm in our case, whereas the impact of the AGN dust torus on the total dust continuum emission in

Tripodi et al. (2022) is limited ≲50 µm. As a result, the peak of our gray-body model is at a longer wavelength, which leads to a lower dust temperature and a higher dust mass according to Eqs. (1) and (2), as well as a smaller gas-to-dust ratio (GDR) compared with the dust property results in Tripodi et al. (2022). The FIR and IR luminosities are 8.84$_{-2.19}^{+2.68}$ × 10¹² L☉ and 1.19$_{-0.29}^{+0.40}$ × 10¹³ L☉, respectively, by integrating the gray-body model (which can be taken as purely star-forming heated dust emission) from 42.5–122.5 µm and from 8–1000 µm. The average SFR is 2055$_{-499}^{+691}$ M☉ yr⁻¹ calculated from the above-mentioned IR luminosity and the formula in Kennicutt (1998).

With the resolved dust continuum at two different frequencies, we are able to investigate the radial distributions of the dust temperature, dust mass and SFR. We first used circular annuli to do aperture photometry on both the [C II] and CO(9–8) underlying dust continuum maps (after removing the central point components that may come from the AGN dust torus) with ring width of 0.1". The right panel of Fig. 5 shows the comparison of the surface brightness of the continuum at observed frame 262 and 147 GHz. The slight differences in the dust continuum profiles at different wavelengths may indicate different distribution of the dust temperature along with distance. For each ring we fitted the measured flux densities with Eqs. (1) and (2) using the $\beta$ value of 1.90 from our SED fitting, to get the average dust temperature, dust mass and dust optical depth at each radius. Finally, with fixed gray-body formula for each ring, we calculated the IR luminosity and SFR. The results are shown in Fig. 9. The dust temperature drops with increasing radius as seen in the top-left panel. A similar dependence is found in a star-forming galaxy at $z = 7.13$ (Akins et al. 2022). We fit a power law, $T_{\text{dust}} \propto r^{-0.41 \pm 0.20}$. The radial dependence of $\Sigma_{\text{dust\_mass}}$ and $\Sigma_{\text{SFR}}$ are consistent with
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Image decomposition of the CO (9–8) line. **Top-left panels:** 2D elliptical Sérsic model and the residual between the observed and modeled CO (9–8) intensity maps, from left to right. We measured the residual rms inside the dashed gray square with a side length of 1". **Bottom-left panels:** same as in the top panels but for the Point+2D elliptical Sérsic modeling. The peak values and rms within the dashed gray square are 0.061 and 0.021, and 0.059 and 0.020 Jy beam⁻¹ km s⁻¹ for the Sérsic and P+Sérsic residual maps, respectively. The shape of the CO (9–8) synthesized beam with a FWHM size of 0.187 × 0.153 is plotted in the bottom-left corner of each panel. **Top-right panel:** distributions of the residuals for pixels inside the dashed gray squares (Sérsic modeling: blue histogram; P+Sérsic modeling: red histogram). **Middle-right panel:** CO (9–8) luminosity surface brightness (black squares, blue circles, and red diamonds with error bars are measurements from the observed intensity map, the 2D elliptical Sérsic model, and the Point+2D elliptical Sérsic model, respectively) at different radii, measured using elliptical rings with the ring width along the major axis half (0.1') that of the CO (9–8) clean beam size, the rotation angle equal to PA (=198°), and the ratio of semimajor and semimajor axis – *b/a* of cos(̃i) (̃i = 43°), where PA and ̃i come from the CO (9–8) line kinematic modeling (listed in Table 3). The error bar represents the deviation of the values of all pixels in each ring. The solid and dashed gray lines are three times the CO (9–8) luminosity surface brightness limit measured in the emission-free region with an elliptical annulus that is the same with the ones used to measure the CO (9–8) luminosity surface brightness but has a larger radius, and its corresponding rms. **Bottom-right panel:** luminosity surface brightness difference between the measurements from the observed intensity map and the ones from the 2D elliptical Sérsic model (blue circles) and from the Point+2D elliptical Sérsic model (red diamonds). Note that the vertical scale is linear and in units of 10⁶ L⊙ kpc⁻².</figcaption></figure>

<table><caption>**Table 3.** Kinematic parameters derived from the ³D BAROLO modeling.</caption><thead><tr><th>Data cube</th><th><i>i</i><br>(°)</th><th>PA<br>(°)</th><th>̃i<br>(°)</th><th>PÃ<br>(°)</th><th><i>V</i><sub>flat</sub><br>(km s<sup>-1</sup>)</th><th><i>V</i><sub>max</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>ext</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>med</sub><br>(km s<sup>-1</sup>)</th><th><i>V</i><sub>flat</sub>/σ<sub>ext</sub></th><th><i>V</i><sub>max</sub>/σ<sub>med</sub></th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th></tr></thead><tbody><tr><td>[C II]</td><td>(37–45)<sup>+1(6)</sup><sub>−2(8)</sub></td><td>(191–214)<sup>+3(26)</sup><sub>−3(23)</sub></td><td>42 ± 3</td><td>199 ± 8</td><td>250<sup>+22</sup><sub>−13</sub></td><td>283<sup>+45</sup><sub>−47</sub></td><td>37<sup>+13</sup><sub>−11</sub></td><td>54<sup>+6</sup><sub>−6</sub></td><td>7<sup>+2</sup><sub>−2</sub></td><td>5<sup>+1</sup><sub>−1</sub></td></tr><tr><td>CO (9–8)</td><td>(41–43)<sup>+2(9)</sup><sub>−3(9)</sub></td><td>(195–202)<sup>+7(9)</sup><sub>−5(9)</sub></td><td>43 ± 1</td><td>198 ± 3</td><td>–</td><td>215<sup>+13</sup><sub>−14</sub></td><td>–</td><td>64<sup>+8</sup><sub>−9</sub></td><td>–</td><td>4<sup>+1</sup><sub>−1</sub></td></tr></tbody></table>

**Notes.** Column 1: the data cube name. Columns 2–3: the inclination angle and position angle for each ring. Columns 4–5: the mean values of the inclination angle and position angle shown in Cols. 2–3. Note that the errors only represent the standard deviations of the measured values for all rings. Column 6: the rotation velocity in the flat part. Column 7: the maximum rotation velocity. Column 8: the velocity dispersion in the flat part. Column 9: the median velocity dispersion. Columns 10–11: the ratios between the rotation velocity and the velocity dispersion.

an exponential distribution $\Sigma(r) = \Sigma_0 \exp^{-r/r_s}$, where $r_s$ is the exponential scale length and $\Sigma_0$ is the central surface density. The best-fit relation for $\Sigma_{\text{dust\_mass}}$ is shown in the top-right panel of Fig. 9 with $r_{s, \text{dust\_mass}} = 0.77 \pm 0.27$ kpc, corresponding to an effective radius of $1.29 \pm 0.45$ kpc. It is about two times larger than the dust continuum effective radius (i.e., ~0.6 kpc) shown in

Table 2. This may indicate that the single-band dust continuum emission is not a good tracer of the dust mass, which is consistent with the finding from numerical simulations (Popping et al. 2022), and that the dust emission depends more strongly on dust temperature than on dust mass. Thus, one requires at least two bands of high-resolution imaging to map the dust temperature
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** SED decomposition toward J2310+1855. **Left panel:** rest-frame UV-to-FIR SED fitting. The red points with error bars or downward arrows are observed data (see Table 3 in Shao et al. 2019 for details). The pink lines represent the UV/optical power law from the accretion disk. The brown lines correspond to the CAT3D AGN torus model (Hönig & Kishimoto 2017). The physical properties derived for these two components are detailed in Table 4 in Shao et al. (2019). The green lines correspond to a gray-body profile (Eqs. (1) and (2)) associated with star formation activity in the quasar host galaxy. The black lines are the sum of all components. The fit employed the MCMC method with the emcee package (Foreman-Mackey et al. 2013). We visualized the model uncertainties with shaded areas by randomly selecting 100 models from the parameter space. **Right panel:** corner map for the parameters (dust temperature, dust mass, and dust emissivity index) associated with the gray-body model for the FIR region SED fitting (green lines in the **left panel**). The contours are drawn at $1 - \exp(-m^2/2)$ ($m = 0.5, 1, 1.5, 2$) of the volume. The vertical dashed lines show 16th, 50th (median), and 84th percentiles.</figcaption></figure>

across the galaxy disks when using the dust continuum emission of galaxies as a reliable tracer of the dust mass distribution. The total dust mass summing over all rings is $(1.84 \pm 0.69) \times 10^9 M_{\odot}$, which is consistent with the value measured from the SED fitting. An exponential fit for $\Sigma_{SFR}$ is shown in the bottom-right panel of Fig. 9 with $r_{s, SFR} = 0.38 \pm 0.09$ kpc, corresponding to an effective radius of $0.64 \pm 0.15$ kpc. It is consistent with the dust effective radius, which means that the resolved dust continuum emission is very closely linked to the SFR distribution. These are consistent with the conclusions from TNG50 star-forming galaxy simulations that the single-band dust emission is a less robust tracer of the dust distribution, but is a decent tracer of the obscured star formation activity in galaxies (Popping et al. 2022). The total SFR summing over all rings is $3415_{-1777}^{+5746} M_{\odot} \text{ yr}^{-1}$, which is consistent with the value measured from the SED fitting. We suggest adopting the total SFR derived from the SED fitting, as the SFR value for each ring is measured from only two data points at different wavelengths and thus has a lot of uncertainty. As shown in the bottom-left panel of Fig. 9, as the dust mass surface density decreases with the galactocentric radius, the dust emission in both bands becomes less optically thick. The rest-frame wavelength of both dust bands (~290 µm for the CO (9–8) underlying dust continuum, and ~158 µm for the [C II] underlying dust continuum) is shortward of the Rayleigh-Jeans tail (i.e., a rest-frame wavelength of ~350 µm). The dust mass surface density is very high (i.e., $\Sigma_{\text{dust\_mass}} > 10^8 M_{\odot} \text{ kpc}^{-2}$), and the optical depth for the dust continuum emission at both frequencies is above 0.1 inside ~1 kpc.

### 4.3. Gas kinematic modeling

The [C II] and CO (9–8) lines can be used to trace the kinematic properties of atomic, ionized, and molecular gas in the quasar host galaxy. As shown in Fig. A.1, the [C II] emission peak

moves in a circular path with frequency. As shown in Figs. 10 and C.2, the velocity fields traced by both [C II] and CO (9–8) show clear velocity gradients. In addition, the position-velocity diagrams along the kinematic major axes have an “S” shape (especially apparent for the [C II] line; however, the resolution of the CO data is roughly two times lower than that of the [C II] data, so the “S” shape is not obvious). These are consistent with a rotating gas disk. We are able to construct a 3D disk model for these data cubes with the package called 3D-Based Analysis of Rotating Object via Line Observations ($^{3D}$BAROLO$^3$; Di Teodoro & Fraternali 2015).

The $^{3D}$BAROLO software fits 3D tilted ring models to spectroscopic data cubes. For each ring, the algorithm builds a 3D disk model of the gas distribution in both the spatial and velocity axes, and then convolves it with the restoring beam of the observed data cube. Finally, it compares the convolved data set with the observed one. The geometry of the tilted ring model can be seen in Fig. C.1 and the corresponding geometrical and kinematic parameters (e.g., the inclination angle $i$; the position angle $\phi$; the rotation velocity $V_{\text{rot}}$) are described in Appendix C. The 3D tilted ring model performs better than the standard 2D modeling on the velocity fields. For example, a common problem when deriving the kinematic properties from the velocity fields is the beam smearing effect, which smears the steep velocity gradient especially in the central region. And there exists differences among the velocity fields measured using different methods (e.g., an intensity-weighted velocity field from CASA, versus a mean gas velocity map based on Gaussian fitting using AIPS$^4$). The $^{3D}$BAROLO algorithm avoids the beam smearing with the convolution step and directly conducts the modeling on the data cube. In addition to the rotation velocity for each

<footnote>$^3$ http://editeodoro.github.io/Bbarolo/</footnote>

<footnote>$^4$ http://www.aips.nrao.edu/</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Radial distributions of the dust temperature (top left), dust mass surface density (top right), optical depth (bottom left) at [C II] and CO (9–8) underlying dust continuum frequencies (red dots and blue squares with error bars, respectively), and the SFR surface density (bottom right). The purple, gold, and green lines are best-fit lines. Note that the uncertainties of the optical depth at different wavelengths are same in the bottom-left panel but appear different due to different ranges of the left and right axes. The error bars of the SFR surface density (which is associated with the dust temperature and dust mass) in the bottom-right panel are measured with a Monte Carlo method: we calculate samples by changing the parameter values with random Gaussian draws centered on their best-fit values and deviated by their errors shown in the top-left and top-right panels, and the lower and upper values are taken as the 16th and 84th percentiles.</figcaption></figure>

ring, 3D BAROLO can also measure the asymmetric-drift correction ($V_A$), which is caused by the random motions and can be directly measured with the velocity dispersion and the density profile (e.g., Iorio et al. 2017). There is no direct way to estimate parameter errors, and 3D BAROLO adopts a Monte Carlo method: it calculates models by changing the parameter values with random Gaussian draws centered on the minimum of the function, once the minimization algorithm has converged.

We took the ring width $W_{\text{ring}}$ as half of the clean beam size (i.e., Nyquist sampling): 0.05 and 0.1 for the modeling of the [C II] and CO (9–8) lines, respectively. The central radius of each ring is $W_{\text{ring}} \times t$ where t is the t-th ring. Eight rings are adopted for the kinematic modeling on our high-resolution [C II] data. The [C II] emission extends to a radius of ~2.8 kpc as shown by the [C II] surface brightness distribution (top-right panel in Fig. 6). We chose to use data within a radius of ~2.4 kpc, taking advantage of higher S/N data. We should note that as the projected [C II] kinematic minor axis (~2.5 kpc) is longer than the kine-

matic major axis, we will lack some information near the kinematic major axis for the outer two rings. However, since we have plenty of information from near the kinematic minor axis and in the region between these two axes, we still have enough data to carry out the kinematic modeling in the outer two rings. The initial guess of the position angle and the inclination angle are 200° and 40°. The initial position angle are roughly measured from the velocity maps based on the definition of the position angle (see Appendix C). One method for roughly determining the inclination angle, assuming an intrinsic round disk, uses the photometric minor and major source size ratio ($r_{\text{size}}$: $i = \cos^{-1}(r_{\text{size}})$). As the kinematic major axis of the [C II] line is far from its photometric major axis and [C II] kinematic major axis is shorter than its kinematic minor axis, but the kinematic major axis of the CO (9–8) line is close to its photometric major axis (see Table 2) and CO (9–8) kinematic major axis is longer than its kinematic minor axis, we derived the initial inclination angle from the CO (9–8) minor and major source size ratio. In addition,
during the fitting, we adopt pixel-by-pixel normalization, which allows the code to exclude the parameter of the surface density of the gas from the fit. It means that we force the value of each spatial pixel along the spectral dimension in the model to equal that in the observations, which allows a non-axisymmetric model in density and avoids untypical regions (Lelli et al. 2012). We present the fitted parameters in Table 3. The average inclination angle and position angle derived from the [C II] and CO (9–8) lines are consistent. We plot the intensity, velocity, velocity dispersion, and position-velocity maps and the line spectra measured from the modeled line data cubes in Figs. 10 and C.2 for the [C II] and CO (9–8) lines, respectively. The rotation velocity, gas velocity dispersion, and asymmetric-drift correction for each ring are shown in Fig. 11. The circular rotation velocities (the rotation velocities corrected by the asymmetric drifts; see Eq. (3)) presented in Fig. 12 for [C II] and CO (9–8) lines are consistent with each other.

### 4.4. Rotation curve decomposition

In order to quantify the dynamical contribution of each matter component, we performed a decomposition of the circular rotation curve measured from the high-resolution [C II] line.

The circular velocity ($V_c$) directly traces the galactic gravitational potential ($\Phi$) and can be measured by correcting the rotation curve for the asymmetric drift:

$$R \frac{\partial \Phi}{\partial R} = V_c^2 = V_{\text{rot}}^2 + V_A^2, \qquad (3)$$

where $V_{\text{rot}}$ is the rotation velocity, $V_A$ is the asymmetric drift correction caused by random motions and can be modeled given the velocity dispersion and the density profile (e.g., Iorio et al. 2017) by $^{3D}$BAROLO, and $\Phi$ is the sum of the potentials of the different mass components.

We consider four matter components (black hole, stellar, gas, and dark matter) that contribute to the total gravitational potential of the quasar-host system. Thus, the circular velocity can be expressed as

$$V_c = \sqrt{V_{\text{BH}}^2 + V_{\text{star}}^2 + V_{\text{gas}}^2 + V_{\text{DM}}^2}, \qquad (4)$$

where $V_{\text{BH}}$, $V_{\text{star}}$, $V_{\text{gas}}$ and $V_{\text{DM}}$ are the contributions of the black hole, stellar, gas and dark matter components to the circular velocity. The following applies for each component at a radius $R$:

First, the Keplerian velocity due to the central black hole – $V_{\text{BH}}$ is

$$V_{\text{BH}} = \sqrt{\frac{GM_{\text{BH}}}{R}}, \qquad (5)$$

where $G$ and $M_{\text{BH}}$ are the gravitational constant and the black hole mass, respectively. We measured a circular velocity of $249_{-12}^{+9} \text{ km s}^{-1}$ at the innermost ring with a central radius of 0.29 kpc. Only considering the gravitational potential caused by the central black hole at the radius of 0.29 kpc, we calculated a black hole mass upper limit of $4.18_{-0.39}^{+0.31} \times 10^9 M_{\odot}$. The black hole mass for our target J2310+1855 measured from the Mg II and C IV lines is $(3.15-5.19) \times 10^9 M_{\odot}$ (Jiang et al. 2016), which corresponds to a $V_{\text{BH}}$ range of $216-277 \text{ km s}^{-1}$ at $R = 0.29 \text{ kpc}$. The black hole contribution to the total potential is significant in the innermost region, but some values of $V_{\text{BH}}$ exceed the measured value at 0.29 kpc. We set the black hole mass as a free parameter in the dynamical modeling. The black hole sphere of influence ($r_h$) is defined as the region within which the gravitational

potential of the black hole dominates over that of the surrounding stars, which can be measured with $r_h = GM_{\text{BH}}/\sigma_{\text{star}}^2$ where $\sigma_{\text{star}}$ is the velocity dispersion of the surrounding stellar population. When $M_{\text{BH}} = 2.97 \times 10^9 M_{\odot}$ (from this dynamical modeling) and $\sigma_{\text{star}} = 250 \text{ km s}^{-1}$ (see Sect. 5.5.2), $r_h = 0.21 \text{ kpc}$. When $M_{\text{BH}} = 10^9 M_{\odot}$ (the order value of the black hole mass measured from NIR spectral lines with the local scaling relations; Jiang et al. 2016; Feruglio et al. 2018) and $\sigma_{\text{star}} = 100 \text{ km s}^{-1}$ (assuming it is equal to the maximum velocity dispersion measured from our kinematic modeling of the gas), $r_h = 0.43 \text{ kpc}$. Our high-resolution ALMA [C II] observations can zoom into the sphere of influence considering that our innermost point is at 0.29 kpc and, thus, can be used to derive the dynamical mass of the central black hole. Beyond the sphere of influence, the gravitational dominance of the black hole quickly vanishes, which is shown as the yellow line of the left and middle panels of Fig. 12. We are not spatially resolving the sphere of influence considering the ratio $r_h/r_{\text{res}}$ (i.e., 0.4–0.8 in this work) between the radius of the black hole sphere of influence and the spatial resolution of the data. However, gas dynamical studies (e.g., Ferrarese & Merritt 2000; Graham et al. 2001; Marconi & Hunt 2003; Valluri et al. 2004) have addressed that resolving the sphere of influence is an important but not sufficient factor to dynamically estimate the black hole mass, and the ratio $r_h/r_{\text{res}}$ can be taken as a rough indicator of the quality of the black hole mass estimate. Ferrarese & Ford (2005) have summarized a list of black hole mass detection based on resolved dynamical studies in their Table II, which shows a $r_h/r_{\text{res}}$ range of 0.4–1700. Our $r_h/r_{\text{res}}$ ratio is just within the above-mentioned range. In summary, we are able to derive a dynamical mass of the black hole from our [C II] data.

Second, the stellar component can be described by a Sérsic profile (e.g., Terzić & Graham 2005; Rizzo et al. 2021), giving rise to

$$V_{\text{star}} = \sqrt{\frac{GM_{\text{star}} \gamma(n_{\text{star}}(3-p), b(R/R_{e; \text{star}})^{1/n_{\text{star}}})}{R \Gamma(n_{\text{star}}(3-p))}}, \qquad (6)$$

where $M_{\text{star}}$, $R_{e; \text{star}}$ and $n_{\text{star}}$ are the total stellar mass, the effective radius and the Sérsic index, $\gamma$ and $\Gamma$ are the incomplete and complete gamma functions. $p$ and $b$ are related to the Sérsic index by $p = 1.0 - 0.6097/n_{\text{star}} + 0.05563/n_{\text{star}}^2$ (when $0.6 < n_{\text{star}} < 10$, and $10^{-2} \le R/R_{e; \text{star}} \le 10^3$) and $b = 2n_{\text{star}} - 1/3 + 0.009876/n_{\text{star}}$ (when $0.5 < n_{\text{star}} < 10$). Due to the limited resolution (~600 pc) of our ALMA data (giving us few data points), the lack of resolved rest-frame optical or NIR data to constrain the two stellar components (i.e., bulge and disk), and the strong degeracies between the two subcomponents, following the fitting method from Rizzo et al. (2021), we adopted a single Sérsic element to globally describe the stellar component.

Third, the gas component can be represented by a thin exponential disk (Binney & Tremaine 1987). Thus,

$$V_{\text{gas}} = \sqrt{\frac{2GM_{\text{gas}}}{R_{\text{gas}}}y^2[I_0(y)K_0(y) - I_1(y)K_1(y)]}, \qquad (7)$$

where $M_{\text{gas}}$ and $R_{\text{gas}}$ are the total gas mass and the gas scale radius, $I_{ii}(y)$ and $K_{ii}(y)$ are the first and second kind modified Bessel functions of zeroth ($ii = 0$) and first ($ii = 1$) orders, and $y$ is given by $y = R/(2R_{\text{gas}})$. The CO (9–8) image decomposition in Sect. 4.1 suggests an approximately exponential distribution of the molecular gas. We fixed the Sérsic index to be 1 for the CO (9–8) gas distribution, and got a scale
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Observed and modeled [C II] line of J2310+1855. **Upper-left panels:** intensity, velocity, and velocity dispersion maps for the observed data from top to bottom, which are labeled “D.” **Upper-central panels:** Same as the left panels but for the best-fit modeling from $^{3D}$BAROLO, labeled “M”. **Upper-right panels:** same as the left panels but for the residuals between the observed data (“D”) and the modeled ones (“M”), which are labeled “R” (note that the dynamical range of these panels are different from others). The black cross in each panel marks the center of the rotating gas disk. The shape of the synthesized beam with a FWHM size of 0.111 × 0.092 is plotted in the bottom-left corner of each panel. In the velocity field panels, the dashed black lines are the kinematic major axis of the gas disk, and the plotted solid gray dots represent the ring positions. **Lower-left and lower-middle panels:** position-velocity maps extracted along the kinematic major and minor axes from the observed data cube (blue contours) and the $^{3D}$BAROLO model cube (red contours). The contour levels are [-2, 2, 6, 10] × 0.15 mJy beam$^{-1}$ for both the data and the model. **Lower-right panel:** [C II] spectra extracted from the observed data cube (blue histogram) and the $^{3D}$BAROLO model cube (red histogram). The residual spectrum between the data and the model is shown as the black histogram. The spectral resolution is 62.5 MHz, corresponding to 68 km s$^{-1}$. In the $^{3D}$BAROLO modeling, we used a re-binned (four original channels) data cube in order to optimize the data S/N per frequency and velocity bins and the sampling of the FWHM of the [C II] line.</figcaption></figure>

radius of ~0.78 kpc, which can be taken as the overall molecular gas scale radius under the assumption of an exponential gas disk. With the detected CO (2–1) line toward J2310+1855 ($vS_v = 0.18 \pm 0.02$ Jy km s$^{-1}$; Shao et al. 2019), and assuming

$L'_{CO(2-1)} \approx L'_{CO(1-0)}$ (Carilli & Walter 2013), we are able to measure the gas mass with an assumed CO-to-gas mass conversion factor ($\alpha_{CO}$), namely: $M_{gas} = \alpha_{CO} L'_{CO(1-0)}$. The only free parameter for the gas component is $\alpha_{CO}$.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Measured rotation velocity (left panel), gas velocity dispersion (middle panel), and asymmetric-drift velocity (right panel) as a function of radius (the central radius of each ring) measured from the 3DBAROLO modeling for the [C II] (red points with error bars) and CO (9–8) (blue points with error bars) lines, respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Best-fit of the decomposition of the circular rotation curve traced by the [C II] line when allowing the black hole component to be free (left), fixed to $1.8 \times 10^9 M_{\odot}$ (Feruglio et al. 2018; middle) and none (right). Black points with error bars are the circular velocities, these are the rotation velocities corrected for the asymmetric drift (caused by gas random motions). The solid yellow, green, brown and blue lines represent the black hole, stellar, gas and dark matter components, respectively. The dashed red lines are the sum of these four components. The blue diamonds with error bars are the circular velocities measured from the CO (9–8) line. These two lines trace the identical gravitational potential within the errors.</figcaption></figure>

Finally, the dark matter component contributes much less to the total gravitational potential than the baryonic components in the innermost regions. However, it is required to fit a flat rotation curve in the outer parts in the Galaxy and other local galaxies (e.g., Carignan et al. 2006; Sofue et al. 2009). We consider the dark matter halo as a Navarro-Frenk-White (NFW; Navarro et al. 1996) spherical halo. These assumptions lead to

$$V_{\text{DM}} = \sqrt{\frac{GM_{\text{DM}}}{R_{\text{DM}}} \frac{1}{x} \frac{\ln(1+cx) - cx/(1+cx)}{\ln(1+c) - c/(1+c)}}, \quad (8)$$

where $M_{\text{DM}}$ and $R_{\text{DM}}$ are the virial mass and radius ($M_{200}$ and $R_{200}$ in Navarro et al. 1996), respectively, c is the concentration parameter, and $x = R/R_{200}$ is the radius in units of the virial radius. The virial mass and radius are correlated with the critical density ($\rho_{\text{crit}}$): $M_{200} = 200\rho_{\text{crit}}(4\pi/3)R_{200}^3$, where $\rho_{\text{crit}} = 3H_z^2/(8\pi G)$ and $H_z = H_0 \sqrt{\Omega_m(1+z)^3 + \Omega_{\Lambda}}$ (where z is the redshift). We calculate the concentration parameter (c) from the mass-concentration relation estimated in N-body cosmological simulations by Ishiyama et al. (2021). The only free parameter for the dark matter component is $M_{\text{DM}}$.

In summary, for the dynamical modeling, we have six free parameters ($M_{\text{BH}}$, $M_{\text{star}}$, $n_{\text{star}}$, $R_{\text{e; star}}$, $\alpha_{\text{CO}}$, $M_{\text{DM}}$) and eight data points (the circular velocities corrected by the asymmetric drift shown in Fig. 12, which have already been corrected for the inclination angle from 3DBAROLO). During the fitting

with the emcee package (Foreman-Mackey et al. 2013), a loose prior constraint of $[10^8, 10^{10}] M_{\odot}$ is adopted for $M_{\text{BH}}$. The quantities of $M_{\text{star}}$, $n_{\text{star}}$ and $R_{\text{e; star}}$ are tightly coupled together, and we use prior ranges of $[10^7, 10^{11}] M_{\odot}$, [0.5, 10], and [0, 2.5] kpc, respectively. As for $\alpha_{\text{CO}}$, we used a range of [0.2, 14] $M_{\odot}/(\text{K km s}^{-1} \text{ pc}^2)$, which is measured from a sample of nearby AGN, ultra-luminous infrared galaxies, and starburst galaxies (Mashian et al. 2015). For $M_{\text{DM}}$, we consider a uniform prior of $10^{10}-10^{13} M_{\odot}$ dark matter halo.

We present the best-fit decomposition for the circular rotation curve in the left panel of Fig. 12, and the physical parameters measured from the dynamical modeling in Col. (a) of Table 4. The black hole mass of $2.97^{+0.51}_{-0.77} \times 10^9 M_{\odot}$ is consistent with the measurements – (4.17 ± 1.02) × 10⁹ and (3.92 ± 0.48) × 10⁹ $M_{\odot}$ from the Mg II and C IV lines (Jiang et al. 2016), and (1.8 ± 0.5) × 10⁹ $M_{\odot}$ from the C IV line (Feruglio et al. 2018). The [C II] emission can be taken as a molecular gas mass tracer in galaxies (e.g., Zanella et al. 2018; Madden et al. 2020; Vizgan et al. 2022) with a conversion factor ($\alpha_{[C II]}$) by $M_{\text{gas}} = \alpha_{[C II]} L_{[C II]}$. The derived gas mass corresponds to a $\alpha_{[C II]}$ of $3.2^{+0.6}_{-0.8} M_{\odot}/L_{\odot}$. Given that there are only two degrees of freedom, and the strong degeneracies among some parameters (i.e., $M_{\text{star}}$, $n_{\text{star}}$ and $R_{\text{e; star}}$), we cannot constrain most of these parameters well with the ~600 pc resolution data we have. We tried another experiment by fixing the black hole mass to a smaller value of $1.8 \times 10^9 M_{\odot}$ (Feruglio et al. 2018), and present the dynamical results in the middle panel of Fig. 12 and Col. (b)
<table><caption>**Table 4.** Derived parameters from the [C II] dynamical modeling.</caption><thead><tr><th>Parameter</th><th></th><th>(a)</th><th>(b)</th><th>(c)</th></tr></thead><tbody><tr><td>M<sub>BH</sub> (10<sup>9</sup> M<sub>☉</sub>)</td><td>(1)</td><td>2.97<sup>+0.51</sup><sub>-0.77</sub></td><td>(<sup>aa</sup>)1.8</td><td>—</td></tr><tr><td>M<sub>star</sub> (10<sup>9</sup> M<sub>☉</sub>)</td><td>(2)</td><td>1.16<sup>+6.51</sup><sub>-1.08</sub></td><td>5.76<sup>+5.89</sup><sub>-3.53</sub></td><td>6.30<sup>+6.25</sup><sub>-2.62</sub></td></tr><tr><td>R<sub>e, star</sub> (kpc)</td><td>(3)</td><td>1.31<sup>+0.74</sup><sub>-0.70</sub></td><td>1.03<sup>+0.94</sup><sub>-0.77</sub></td><td>0.21<sup>+0.55</sup><sub>-0.17</sub></td></tr><tr><td>n<sub>star</sub></td><td>(4)</td><td>4.88<sup>+3.47</sup><sub>-3.18</sub></td><td>6.66<sup>+2.34</sup><sub>-2.64</sub></td><td>7.17<sup>+2.02</sup><sub>-3.25</sub></td></tr><tr><td>M<sub>DM</sub> (10<sup>10</sup> M<sub>☉</sub>)</td><td>(5)</td><td>4.55<sup>+32.28</sup><sub>-4.21</sub></td><td>1.63<sup>+26.66</sup><sub>-1.58</sub></td><td>3.42<sup>+30.62</sup><sub>-3.34</sub></td></tr><tr><td>α<sub>CO</sub> [M<sub>☉</sub>/(K km s<sup>-1</sup> pc<sup>2</sup>)]</td><td>(6)</td><td>0.37<sup>+0.07</sup><sub>-0.09</sub></td><td>0.37<sup>+0.09</sup><sub>-0.10</sub></td><td>0.34<sup>+0.09</sup><sub>-0.09</sub></td></tr><tr><td>M<sub>gas</sub> (10<sup>10</sup> M<sub>☉</sub>)</td><td>(7)</td><td>2.05<sup>+0.40</sup><sub>-0.52</sub></td><td>2.02<sup>+0.47</sup><sub>-0.55</sub></td><td>1.84<sup>+0.47</sup><sub>-0.49</sub></td></tr><tr><td>M<sub>bary</sub> (10<sup>10</sup> M<sub>☉</sub>)</td><td>(8)</td><td>2.37<sup>+0.50</sup><sub>-0.52</sub></td><td>2.74<sup>+0.45</sup><sub>-0.52</sub></td><td>2.65<sup>+0.42</sup><sub>-0.52</sub></td></tr><tr><td>f<sub>gas</sub></td><td>(9)</td><td>0.95<sup>+0.05</sup><sub>-0.25</sub></td><td>0.77<sup>+0.13</sup><sub>-0.18</sub></td><td>0.74<sup>+0.11</sup><sub>-0.18</sub></td></tr><tr><td>V<sub>200</sub> (km s<sup>-1</sup>)</td><td>(10)</td><td>111<sup>+112</sup><sub>-65</sub></td><td>79<sup>+125</sup><sub>-54</sub></td><td>101<sup>+116</sup><sub>-72</sub></td></tr><tr><td>R<sub>200</sub> (kpc)</td><td>(11)</td><td>15.85<sup>+15.99</sup><sub>-9.21</sub></td><td>11.25<sup>+17.91</sup><sub>-7.77</sub></td><td>14.42<sup>+16.59</sup><sub>-10.24</sub></td></tr></tbody></table>

**Notes.** Columns (a)–(c) are for the dynamical models when allowing the black hole component to be free, fixed and none, respectively. Rows 1–6: fitted parameters from the dynamical modeling in Sect. 4.4. They are the black hole mass, the stellar mass, the effective radius of the stellar component, the Sérsic index of the stellar component, the dark matter mass, and the CO-to-gas mass conversion factor, respectively. Row 7: the gas mass ($M_{gas} = \alpha_{CO}L'_{CO(1-0)}$). Row 8: the baryonic mass ($M_{bary} = M_{star} + M_{gas}$). Row 9: gas mass fraction ($f_{gas} = M_{gas}/M_{bary}$). Rows 10–11: the virial velocity and virial radius of the dark matter halo. (aa) The black hole mass is fixed to the value of $1.8 \times 10^9 M_{\odot}$ (Feruglio et al. 2018).

of Table 4. In both situations, the derived stellar mass is on the order of $10^9 M_{\odot}$ (however with large uncertainty). And when fixing the black hole component, the stellar component (green line in the middle panel of Fig. 12) has a bump in the inner region. This may indicate that a massive stellar bulge already formed at $z = 6$.

# 5. Discussion

## 5.1. The spatial distribution and extent of the ISM

### 5.1.1. Comparisons among different ISM tracers

The different distribution and extent of various ISM tracers may be due to dissociation regions heated by starbursts and/or AGN in the quasar-host system of J2310+1855. We model the spatial distribution and extent (e.g., quantified by the half-light radii and the Sérsic index) of the [C II] and the CO (9–8) lines and the dust continuum emission, using a 2D elliptical Sérsic function in Sect. 4.1.

The extended FIR dust continuum and the CO high-J line emission follow a Sérsic distribution with the Sérsic index $n_{ISM}$ a bit larger than 1. Exponential dust-continuum profiles are also observed for nearby (e.g., Haas et al. 1998; Bianchi & Xilouris 2011) and high-redshift galaxies (Hodge et al. 2016; Calistro Rivera et al. 2018; Wang et al. 2019; Novak et al. 2020). However, the [C II] emission profile is slightly flatter, with a Sérsic index of 0.59.

Walter et al. (2022) used simple RADEX (van der Tak et al. 2007) modeling on a $z \sim 7$ quasar host galaxy, and find that the [C II] surface brightness in the very central part (i.e., $<0.4$ kpc) is significantly suppressed by the high dust opacity, which increases the FIR background radiation field on the [C II] line, thereby making the [C II] fainter than one would observe without the FIR background. As shown in Fig. 9, both the dust mass surface density ($\Sigma_{dust\ mass} > 10^8 M_{\odot} \text{ kpc}^{-2}$) and the [C II] underlying dust optical depth ($\tau \sim 0.5$) are very high inside the central 1 kpc region in our quasar host galaxy. In order to test the

high dust opacity effect on the spatial distribution of [C II], following Walter et al. (2022), we model the [C II] line intensities as a function of radius with and without the FIR background determined using the dust emission in the host galaxy (i.e., a gray-body radiation of $B_{V}(T_{dust})(1 - \exp^{-\tau_{V}})$). For both situations, we also consider the CMB as a background (i.e., a black-body radiation of $B_{V}(T_{CMB})$; $T_{CMB} = 2.73 \times (1 + z)$ K). It should be noted that we only model collisional excitation with molecular hydrogen $H_2$, and ignore collisions with electrons and neutral hydrogen $H I$. We model the [C II] line intensities for the whole galaxy by dividing the galaxy disk into nine concentric rings with ring width of 0.05. The number density of $H_2$ is fixed to $10^5 \text{ cm}^{-3}$. The column density of [C II] for each ring can be calculated with a fixed gas-to-dust mass ratio of $\sim 10$ (from our measured total dust mass in Sect. 4.2 and total gas mass in Sect. 4.4) and a fixed [C II] abundance relative to hydrogen ($7 \times 10^{-6}$ is the best value in our fitting). We further assume that the average kinetic temperature is equal to the dust temperature at the central radius of each ring, and the dust temperature can be predicted from the $T_{dust} - r$ relation measured in Sect. 4.2. The average optical depth in each ring is associated with the dust mass surface density, which can be predicted from the $\Sigma_{dust\ mass} - r$ relation measured in Sect. 4.2. In order to match the observed surface brightness of the [C II] line, we reduce the scale length of the $\Sigma_{dust\ mass} - r$ relation by a factor of 2 for the inner 0.5 kpc region, which is consistent with the observed surface density profiles of the dust mass for some nearby galaxies (e.g., Casasola et al. 2017). The full width at half maximum (FWHM) of the [C II] line emission for each ring is measured from our observed [C II] line spectrum, which is extracted from each circular annulus.

The RADEX modeling results are shown in Fig. 13. The high FIR background associated with the dust radiation, which has dust optical depth $\sim 0.5$ inside 1 kpc, indeed lowers the [C II] surface brightness, thus changing the observed spatial distribution of [C II], reducing the [C II] equivalent width, and lowering the surface brightness ratio between [C II] and FIR emission. The effect becomes weaker in the outer parts of the galaxy, as the dust optical depth decreases. As shown in the right panel of Fig. 9, the CO (9–8) underlying dust continuum optical depths are below 0.2. And importantly, the CO (9–8) emits at longer wavelength than [C II], further away from the dust emission peak. Thus, the FIR background radiation is not strong enough to significantly modify the observed CO (9–8) line intensity. In summary, the difference of the observed surface brightness distribution between [C II] (with a Sérsic index of 0.59) and CO (9–8) (with a Sérsic index of 1.21 or 2.01) may be just the effect of high dust opacity (see also Riechers et al. 2014), which diminishes the [C II] emission much more strongly in the center than in the outer region.

We next compare the rotation angles (Col. 8 in Table 2) of the major axes of the tilted brightness distribution of the ISM tracers with the position angles (Col. 5 in Table 3) from the kinematic modeling. The photometric major axis of the [C II] line is close to that of its underlying dust continuum, but is almost perpendicular to that of the CO (9–8) line. The photometric major axis of the CO (9–8) line is coincident with the kinematic major axis of the [C II]/CO (9–8) shown in Sect. 4.3, but makes an angle of $\sim 30^{\circ}$ with respect to that of its underlying dust continuum. Such a misalignment between kinematic and photometric position angles is also reported for the hosts of Palomar-Green quasars by Molina et al. (2021), which may indicate kinematic twisting. We compare the effective radii (in boldface in Cols. 10–11 of Table 2) of the emission near the kinematic major axis, as they are less influenced by projection effects and can be
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** RADEX modeling of the [C II] line emission with (red lines) and without (blue lines) FIR background caused by high dust opacity in the quasar host galaxy, and compared with the observed [C II] line emission (black circles with error bars): the [C II] surface brightness profile (**left panel**), the [C II] line equivalent width (**middle panel**), and the ratio between $\Sigma_{[CII]}$ and $\Sigma_{FIR}$ (**right panel**) as a function of radius. The dashed black line in the **right panel** represents the radial distribution of the optical depth at the frequency of the [C II] underlying dust continuum.</figcaption></figure>

interpreted as the matter distribution along the kinematic major axis. The half-light radius of the [C II] underlying dust continuum at observed frame 262 GHz is larger by a factor of 1.1 than that of the CO (9–8) underlying dust continuum at observed frame 147 GHz. The half-light radii of the dust continuum at both frequencies are smaller than those of their nearby lines, with ratios of ~0.7 and ~0.5 of [C II] underlying dust continuum to the [C II] line, and CO (9–8) underlying dust continuum to the CO (9–8) line, respectively. The half-light radius of the [C II] line is smaller than that of the CO (9–8) line. The more concentrated dust continuum than the [C II] and the CO (9–8) emission is also found in other $z \sim 6$ quasars (e.g., Wang et al. 2019; Novak et al. 2020), high-redshift galaxies (e.g., Riechers et al. 2013, 2014; Chen et al. 2017; Gullberg et al. 2019), and numerical simulations (e.g., Cochrane et al. 2019; Popping et al. 2022).

The higher compactness for the continuum than the lines reflects the temperature dependence on radius; the dust surface brightness depends on both column density and temperature (Eq. (1)), while the lines mainly scale with temperature (e.g., Calistro Rivera et al. 2018). The dust continuum sizes at the rest frame wavelengths of 158 and 290 µm are identical within ~10%, which is consistent with simulations conducted by Popping et al. (2022), who found that the dust continuum sizes remain constant within ~20% at observed frame wavelengths from 500 µm to 2 mm for $z \sim 1–5$ main-sequence galaxies. As these two wavelengths are very close, they are both dominated by dust at similar temperatures.

Our RADEX modeling of the effect of the dust opacity suggests that the actual half-light radius of the [C II] gas is much smaller than that measured from the observed *n* = 0.59 Sérsic distribution listed in Table 2. The intensity of the CO (9–8) line depends strongly on the ISM density that drives the excitation of high-*J* CO lines. The [C II] line is mainly collisionally excited by electron and hydrogen, which are heated by UV photons. The different effective radii of the [C II] and CO (9–8) lines may be a reflection of the radial dependence of the gas excitation and the radiation field.

### 5.1.2. The nuclear and extended components of the ISM

The long-wavelength dust continuum emission appears to trace the bulk of the star formation activity in quasar-starburst systems at $z \sim 6$ (e.g., Venemans et al. 2018; Li et al. 2020c). Detailed SED decomposition from rest-frame UV to FIR allows one to constrain the relative contribution to the continuum of differ-

ent components (i.e., accretion disk, torus, host galaxy; e.g., Leipski et al. 2013, 2014; Shao et al. 2019). We find (Fig. 8) that the AGN dust torus contributes ~40% of the integrated FIR luminosity of J2310+1855, and the star-formation heated dust continuum emission dominates the FIR luminosity above rest-frame 50 µm (Shao et al. 2019; Sect. 4.2 in this paper). The dust continuum decomposition results presented in Table 2 and Figs. B.2 and B.3 reveal both a compact central component and an extended component, which may stem from the AGN dust torus and the quasar host galaxy associated with star formation, respectively. The AGN dusty molecular torus is a key element in the AGN unification model, which obscures the accretion disk and the broad line region in type-2 AGN (e.g., Hönig 2019). The AGN torus might be unstable, warped and unaligned with the host galaxy orientation, and has been revealed by observations of several molecular lines and dust continuum in low-$z$ AGN (i.e., 10 pc scale; e.g., García-Burillo et al. 2014, 2021; Combes et al. 2019). Our image decomposition assigns ~5% of the dust continuum near both [C II] and CO (9–8) to a point source, which we interpret as emission from the AGN torus.

The diameters of the dusty molecular tori measured by high-resolution ALMA observations of nearby Seyfert galaxies range from ~25 pc to ~130 pc, with a median value of ~42 pc (García-Burillo et al. 2021). We here consider a radius of 25 pc of the torus in our quasar, and with the point emission at the two wavelengths listed in Table 2 and Eqs. (1) and (2) ($\beta$ fixed to be 1.90), we calculate a torus dust temperature of 1600 ± 500 K and a torus dust mass of (1.46 ± 1.14) ×10⁶ $M_\odot$. Esparza-Arredondo et al. (2021) explored the torus GDR for a sample of 36 nearby AGN with NuSTAR and Spitzer spectra, and found that it lies between 1 and 100 times the Galactic ratio (i.e., 100; Bohlin et al. 1978). From our dynamical modeling, we get an average GDR of ~10 for the quasar host galaxy. We assume a GDR value of 100 for the torus alone, as this is consistent with the existence of gas located in the dust-free inner region of the torus. Thus, the molecular gas mass derived from the dust emission in the torus is (1.46 ± 1.14) ×10⁸ $M_\odot$.

We present the best-fit CO (9–8) image decomposition results in Fig. 7 and Table 2 with both a nuclear and extended components, representing the CO (9–8) emission from the AGN molecular torus and the AGN host galaxy, respectively. The nuclear component comprises ~11% of the whole CO (9–8) emission. The nuclear CO (9–8) emission can be converted to a molecular gas mass with a flux ratio between CO (9–8) and CO (2–1) lines (i.e., ~8; Shao et al. 2019) and a CO-to-H₂
conversion factor – $\alpha_{CO}$, assuming $L'_{CO(2-1)} \approx L'_{CO(1-0)}$ (Carilli & Walter 2013). We derive an average $\alpha_{CO}$ value of ~0.37 from our dynamical study, but this may be smaller in the core. For example the conversion factor in the circumnuclear disk can be 3–10 lower than that of the ISM (i.e., Usero et al. 2004; García-Burillo et al. 2014). Assuming a $\alpha_{CO}$ that is 10 times smaller than the overall value, the molecular gas mass is derived to be $(2.44 \pm 0.28) \times 10^8 M_{\odot}$, which is roughly consistent with that measured from the dust emission from the AGN dusty molecular torus. However, the decomposed unresolved emission from the AGN dusty molecular torus is quite uncertain given our current angular resolution and sensitivity. And the estimate of the molecular gas mass in the dusty molecular torus depends on many uncertain quantities: the torus size, the dust emissivity index, the GDR, the CO line ratio and the $\alpha_{CO}$ in the unresolved central region.

## 5.2. The surface density of the gas and the star formation

Next, we investigated the surface density of the gas and star formation. We measured the surface densities of gas and dust as a function of distance from the quasar center from the intensity maps using circular rings. These rings have widths of 0.′05 and 0.′1 for the [C II] and CO (9–8) lines, respectively, about half the restoring beam sizes. For each ring, we used Eqs. (1) and (2) with a fixed dust emissivity index of 1.90 (from the UV to FIR SED decomposition in Sect. 4.2) and with the dust temperature and dust mass surface density following the relations from our dust diagnostics in Sect. 4.2. Then we measured the FIR luminosity by integrating from 42.5 to 122.5 µm and the IR luminosity by integrating from 8 to 1000 µm.

### 5.2.1. $\Sigma_{[C\text{II}]}/\Sigma_{\text{FIR}}$

Local spiral galaxies have a typical [C II]-FIR luminosity ratio of ~5 × 10⁻³ when $\Sigma_{\text{FIR}} < 10^{11} L_{\odot} \text{ kpc}^{-2}$. However, the ratio is substantially smaller at higher FIR luminosity surface densities (e.g., Díaz-Santos et al. 2013, 2017). This so-called [C II]-FIR deficit is shown in the left panel of Fig. 14, together with results from other spatially resolved z ~ 6 quasars (e.g., Wang et al. 2019) and z ~ 3 submillimeter galaxies (SMGs; e.g., Rybak et al. 2019). As shown in the right panel of Fig. 13, the ratio of $\Sigma_{[C\text{II}]}$ and $\Sigma_{\text{FIR}}$ of our target increases with increasing radii, appearing to asymptote at larger radii. Our [C II]-FIR ratios are a few times higher than those of spatially resolved star-forming associations at similar FIR surface brightness. Our measurements push to higher $\Sigma_{\text{FIR}}$ than do the comparison samples, but do not probe deeper into the faint end (i.e., $\Sigma_{\text{FIR}} < 10^{10} L_{\odot} \text{ kpc}^{-2}$) occupied by Lyman break galaxies (LBGs; e.g., Hashimoto et al. 2019a) and the local luminous infrared galaxies in the GOALS sample (Díaz-Santos et al. 2013, 2017).

The origin of the overall [C II]-FIR deficit of galaxies has been investigated in detail (e.g., Luhman et al. 1998; Stacey et al. 2010; Díaz-Santos et al. 2013; Narayanan & Krumholz 2017; Lagache et al. 2018), yet remains under debate. The [C II] line is mainly excited through collisions with electrons, neutral hydrogen H I, and molecular hydrogen H₂. Close to the quasar host galaxy center (corresponding to the region with the highest $\Sigma_{\text{FIR}}$), the UV photons are absorbed by large dust grains (which may be reemitted at FIR wavelengths, thus boosting the FIR luminosity). Therefore, the hydrogen and electrons will be less heated, giving rise to few collisions to excite ionized carbon. In addition, the dust absorption reduces

the number of ionizing photons, decreasing the extent to which carbon becomes ionized (e.g., Luhman et al. 2003). Other popular explanations for the [C II]-FIR deficit are the saturation of [C II] emission at very high temperatures (e.g., $T_{\text{gas}} > T_{[\text{C II}]}$; e.g., Muñoz & Oh 2016) or in extremely dense environments (in which more carbon is in the form of CO than C⁺; e.g., Narayanan & Krumholz 2017).

The observed larger [C II]-FIR deficit in the center of our quasar host galaxy may be a reflection of the high dust temperature at small radius, and the different dependence of $\Sigma_{\text{FIR}}$ and $\Sigma_{[\text{CII}]}$ on $T_{\text{dust}}$. As shown in the right panel of Fig. 14, $\Sigma_{\text{FIR}}$ is an integral over the gray-body relating to the dust temperature with a power-law index (~$T_{\text{dust}}^{5.0 \pm 0.4}$; black line), larger than that (~$T^4$) of the Stefan-Boltzmann law. However, $\Sigma_{[\text{CII}]}$ scales as ~$T^{2.4 \pm 0.7}$ (green line). Thus, higher dust temperatures in the galaxy center produce a larger [C II]-FIR deficit (see also Gullberg et al. 2015; Walter et al. 2022). In addition, as demonstrated by Walter et al. (2022) and our RADEX test in Sect. 5.1.1, the high dust opacity, which increases the FIR background radiation field on the [C II] line (see also Riechers et al. 2014), suppresses the observed [C II] line, especially in the central region, thus enhancing the [C II]-FIR deficit.

### 5.2.2. $\Sigma_{M_{\text{gas}}}$ and $\Sigma_{\text{SFR}}$

The Kennicutt-Schmidt (KS) power-law relation between the gas mass and the SFR surface density, $\Sigma_{\text{SFR}} \sim \Sigma_{M_{\text{gas}}}^N$, describes how efficiently galaxies turn their gas into stars, which enables us to understand galaxy formation and evolution across cosmic time. The KS relation is nearly linear ($N \sim 1-1.5$) in galaxies ranging from the local Universe to high redshifts ($z < 4$; e.g., Kennicutt 1998; Bothwell et al. 2010; Genzel et al. 2010, 2013), indicating that the star formation processes may be independent of cosmic time. It is still unknown whether it holds for higher-redshift galaxies up to the reionization epoch.

In order to determine whether the spatially resolved ISM of J2310+1855 at $z = 6.0$ deviates from the star formation law observed in local galaxies and high-z star-forming galaxies, Fig. 15 plots the SFR surface density as a function of the molecular gas mass surface density of the resolved ISM of J2310+1855 (black open circles and diamonds with error bars), and of other populations from local to high-z galaxies in the literature (references in the caption). The molecular clouds where stars form are primarily composed of molecular hydrogen (H₂) and helium. Given cosmological abundances, we consider the mass of the star-forming gas to be the mass of the molecular hydrogen times 1.36. We measure the gas mass surface density $\Sigma_{M_{\text{gas}}}$ in each ring from the CO (9–8) intensity map with a flux ratio between CO (9–8) and CO (2–1) lines of ~8 (Shao et al. 2019) and assume $L'_{\text{CO(2-1)}} \approx L'_{\text{CO(1-0)}}$ (Carilli & Walter 2013) and a conversion factor $\alpha_{\text{CO}} \sim 0.37 M_{\odot}/(\text{K km s}^{-1} \text{ pc}^2)$ from the dynamical modeling in Sect. 4.4. In addition, we also measure $\Sigma_{M_{\text{gas}}}$ from the [C II] intensity map with $\alpha_{[\text{C II}]} \sim 3.2 M_{\odot}/L_{\odot}$ from the dynamical modeling in Sect. 4.4. We similarly measure the SFR surface density $\Sigma_{\text{SFR}}$ from the resolved dust continuum by converting the IR luminosity to SFR using the formula in Kennicutt (1998). Note that we subtracted the IR luminosity contributed by the AGN torus as a point source (see Sects. 4.1 and 5.1.2).

Our SFR and gas surface densities are at the high end of the samples shown in Fig. 15, comparable to the SMGs at $2 < z < 4$. The best-fit KS relation with $\Sigma_{M_{\text{gas}}}$ measured from the CO (9–8)
<figure><img src="image_15.png" /><figcaption>Fig. 14. [C II]-FIR deficit. Left panel: [C II]/FIR luminosity ratio as a function of FIR surface brightness. The red open diamonds with error bars connected by the solid red line represent measurements for J2310+1855. The comparison samples are spatially resolved $z \sim 6$ quasars (green triangles connected by dashed lines; Wang et al. 2019) and $z \sim 3$ SMGs (purple triangles connected by dashed lines; Rybak et al. 2019), other $z > 6$ quasars (yellow stars; Wang et al. 2013; Venemans et al. 2016, 2017; Decarli et al. 2018), other high-$z$ SMGs (pink triangles; Riechers et al. 2013; Neri et al. 2014; Gullberg et al. 2018; Rybak et al. 2019), the GOALS sample of local luminous infrared galaxies (gray circles; Díaz-Santos et al. 2013, 2017), and LBGs (green triangles; Capak et al. 2015; Jones et al. 2017; Hashimoto et al. 2019a). Right panel: Comparison between the trends of the surface brightnesses of [C II] (green circles with error bars) and the scaled FIR (black circles with error bars) luminosities, as a function of dust temperature. The green and black lines are the power-law fits to the $\Sigma_{[C II]}$ and $\Sigma_{FIR}$ as a function of $T_{dust}$, respectively.</figcaption></figure>

line for J2310+1855 is

$$ \frac{\log_{10} \Sigma_{SFR}}{M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}} = -20.50(\pm 1.87) + 2.50(\pm 0.21) \times \frac{\log_{10} \Sigma_{\text{gas}}}{M_{\odot} \text{kpc}^{-2}}, \quad (9) $$

and the best-fit KS relation with $\Sigma_{\text{gas}}$ measured from the [C II] line is

$$ \frac{\log_{10} \Sigma_{SFR}}{M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}} = -10.59(\pm 0.64) + 1.37(\pm 0.07) \times \frac{\log_{10} \Sigma_{\text{gas}}}{M_{\odot} \text{kpc}^{-2}}, \quad (10) $$

for which we do not fit with the inner most measurement, as the observed [C II] emission in the central of the galaxy is highly suppressed due to the high dust opacity (see Sect. 5.1.1).

The power-law index $N = 2.50 \pm 0.21$ measured from the CO(9–8) line is higher than that of local disk galaxies (e.g., Kennicutt 1998; Kennicutt et al. 2007) and high-redshift star-forming populations (e.g., Genzel et al. 2013; Sharon et al. 2019). This may indicate that this quasar host galaxy at $z = 6.00$ is undergoing much faster star formation than do galaxies in the local Universe and at lower redshifts. However, we should note that, as mentioned in other studies (e.g., Daddi et al. 2010; Thomson et al. 2015) the CO-to-H$_2$ conversion factor $\alpha_{CO}$ can vary for different galaxy types, position scales, and metallicities. We adopted the $\alpha_{CO}$ value from our dynamical studies in Sect. 4.4, which only represents the global value, and ignored any possible variation within the galaxy. This will bias the slope of the KS law we found in Fig. 15. In addition, the indirect measurements of the low-$J$ CO line transitions requires assuming a CO excitation model, which brings additional uncertainties when comparing KS relations from different studies. We simply use a bulk ratio between the ALMA CO(9–8) emission and the VLA CO(2–1) emission, which in reality will decrease with radius. This leads to an over-estimation of the slope of the KS law. Low-$J$ CO observations with higher spatial resolution will be needed to improve on this analysis. And there appears to be an excess in CO(9–8) luminosities in distant starbursts

due to for example shock excitation (e.g., Riechers et al. 2021a), which would make the CO(9–8) line a poorer tracer of gas mass.

The power-law index $N = 1.37 \pm 0.07$ measured from the [C II] line is consistent with both the local and other high-redshift samples (e.g., Bothwell et al. 2010; Carilli et al. 2010; Genzel et al. 2010; Thomson et al. 2015). This may indicate that this quasar host galaxy at $z = 6$ in our work follows similar star formation process with that of the local Universe and lower redshifts. However, the value of the $\alpha_{[C II]}$ might be variable within the galaxy. The depletion times defined by $\Sigma_{\text{gas}}/\Sigma_{SFR}$ for the ISM of our quasar host galaxy are 1–100 Myr, increasing from the inner region (corresponding to the highest $\Sigma_{\text{gas}}$ and $\Sigma_{SFR}$) to outer. The short gas consumption timescales suggest a rapid starburst mode for our quasar host galaxy.

## 5.3. The asymmetric line profile

The integrated spectra of the gaseous disks of many galaxies both in the local Universe and at high redshift often have a characteristic double-horned shape (e.g., Roberts 1978; Walter et al. 2008; Shao et al. 2017). Higher velocity dispersion and steeper density profiles of the gas tend to decrease the depth of the valley between the horns, and can even reduce the global profile to flat-topped Gaussian shapes (e.g., de Blok & Walter 2014; Stewart et al. 2014). In addition, the gaseous disk can be disturbed by interactions of galaxies with their neighbors and environments (e.g., van Eymeren et al. 2011; Bok et al. 2019; Watts et al. 2020; Reynolds et al. 2020), ram-pressure stripping (e.g., Gunn et al. 1972; Kenney et al. 2015) and gas accretion from the cosmic web (e.g., Bournaud et al. 2005). These distortions may drive morphological (e.g., the gas distribution) and kinematic (e.g., the velocity field, and the global spectral profile) asymmetries. We now turn to the spectral profile of J2310+1855 to understand it in the context of these ideas.

The [C II] and CO(9–8) spectral profiles of J2310+1855 are asymmetric as shown in Figs. 1 and 2. The spectra
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Spatially resolved KS plot for J2310+1855 (black open circles with error bars, and its best-fit power-law shown as black solid line when calculating Σ<sub>Mgas</sub> from the CO (9–8) line; black open diamonds with error bars, and its best-fit power-law shown as dashed black line when calculating Σ<sub>Mgas</sub> from the [C II] line), compared with the surface density KS relations of other galaxy samples including local starburst galaxies (brown filled circles; Kennicutt 1998), local spiral galaxies (blue filled circles; Kennicutt 1998; Bigiel et al. 2010), z ~ 1–3 star-forming galaxies (pink filled circles; Genzel et al. 2010; Freundlich et al. 2013; Tacconi et al. 2013), z ~ 2–4 SMGs (blue filled squares; Bothwell et al. 2010; Genzel et al. 2010; Carilli et al. 2010), star-forming dusty galaxies (blue open stars; Villanueva et al. 2017), local blue compact dwarf galaxies (red open triangles; Amorín et al. 2016), and some spatially resolved galaxies – Milky Way clumps (green open diamonds; Heiderman et al. 2010; Evans et al. 2014), star-forming regions of the spiral galaxy NGC 5194 (gray filled circles; Kennicutt et al. 2007), the lensed SMG SDP.81 at z = 3.042 (red pluses; Hatsukade et al. 2015), the young low-metallicity starburst galaxy RCSGA 032727 at z = 1.7 (brown crosses; González-López et al. 2017), massive star-forming galaxies around z = 1.2 (purple open diamonds; Freundlich et al. 2013), the z = 1.5 star-forming galaxy EGS13011166 (dashed red line; Genzel et al. 2013), the lensed SMG SMMJ21352 at z = 2.3 (green filled triangles; Thomson et al. 2015), and the H II regions in the nearby spiral galaxy M33 (black crosses; Miura et al. 2014). The solid gray lines represent gas depletion timescales of 1Myr, 10 Myr, 100 Myr, and 1 Gyr. We applied a 1.36 correction factor to the molecular hydrogen mass for the presence of helium for the comparison samples. For our target, the conversion factor to the gas mass is from the dynamical modeling, which should include the contributions from all other elements (e.g., helium). When deriving the index of the KS power-law with the Σ<sub>Mgas</sub> derived from the [C II] line, we do not use the innermost measurement, as the observed [C II] emission in the central of the galaxy is highly suppressed due to the high dust opacity (see Sect. 5.1.1).</figcaption></figure>

are slightly enhanced on the red side. The enhancements
([S<sub>v, redpeak</sub> − S<sub>v, bluepeak</sub>]/S<sub>v, bluepeak</sub>; from the double-Gaussian fit-
ting) are about 26% ± 10% and 18% ± 14% of the [C II] and
CO (9–8) lines, respectively. The two possible peaks are sepa-
rated by 213 ± 9 and 215 ± 14 km s<sup>−1</sup> for the [C II] and CO (9–8)
lines, respectively. From our kinematic modeling of the [C II]
and CO (9–8) data cube in Sect. 4.3, we found a median gas
dispersion of 54<sup>+6</sup><sub>−6</sub> and 64<sup>+8</sup><sub>−9</sub> km s<sup>−1</sup> for the [C II] and CO (9–8)
lines, respectively (Table 3). The ratios between the separation
of the two peaks in the red and blue parts and the median veloc-

ity dispersion are 3.9<sup>+0.5</sup><sub>-0.5</sub> and 3.3<sup>+0.5</sup><sub>-0.5</sub> for the [C II] and CO (9-
8) lines, respectively. In addition, as shown in Sect. 4.3, the
radial distribution of [C II] and CO (9–8) lines are not steep (e.g.,
the Sérsic index is around 1). Thus, we are able to detect the
double-horned profiles for both the [C II] and CO (9–8) lines
(we consider the two peaks to be resolved if the separation of
the two peaks is twice the velocity dispersion). The redward
enhancement of the spectra may be due to more material on one
side of the galaxy, or temperature differentials from one side to
another.
## 5.4. Ionized and molecular gas kinematics

We studied the gas kinematics with both [C II] and CO (9–8) lines using ³DBAROLO in Sect. 4.3. These two lines present consistent kinematic properties, including the gas disk geometry and the gas motions.

### 5.4.1. [C II] versus CO (9–8)

The [C II] and CO (9–8) emission trace similar gas disk geometries in our quasar J2310+1855 host galaxy. As shown in Table 3, the inferred inclination angle and position angle are identical for the two lines. The modest inclination angle (~40°) makes J2310+1855 an ideal target to investigate the gas kinematics in the early Universe. The rotation speed as well as the velocity dispersion (discussed below) are also consistent between the [C II] and CO (9–8) lines. This suggests that the [C II] and CO (9–8) gas are coplanar and corotating in the host galaxy of quasar J2310+1855.

The rotation speed of the [C II] line presented in the left panel of Fig. 11 is roughly constant (~250 km s⁻¹) with radius. However, the rotation speed of the CO (9–8) line at ~1.2 kpc is somewhat lower than that of the other two rings. The measured rotation speed of the [C II] line is overall consistent with that of the CO (9–8) line, given that the errors in Fig. 11 do not include the covariance with other fitting parameters. The [C II] and CO (9–8) data cubes used in the kinematic modeling in Sect. 4.3 are from the same telescope and follow an identical data reduction process. They have a similar velocity resolution (~65 km s⁻¹), which is sufficient to sample the intrinsic width of both lines, and the line sensitivities are also comparable (i.e., S/N<sub>peak</sub> > 10). The only difference is the angular resolution: that of the [C II] data is roughly two times better than that of the CO (9–8) data. Considering that we measured the gas rotation speed in rings with widths that are half the angular resolution for both lines, and that ³DBAROLO performs well even when the galaxy is resolved with a small number of beams (Di Teodoro & Fraternali 2015), the spatial resolution difference may not affect the inferred rotation curve significantly. Thus, the lower rotational speed of the ~1.2 kpc ring from CO (9–8) kinematical modeling may be due to higher random motions of the CO (9–8) gas in that ring. The maximum rotation velocities (~250 km s⁻¹) of J2310+1855 in this work are smaller than that (~400 km s⁻¹) of quasar ULAS J131911.29+095051.4 at z = 6.13 in our previous work (Shao et al. 2017), but are at the high end of the rotation velocities of a sample of 27 z ≥ 6 quasars (Neeleman et al. 2021). The velocity dispersions of the [C II] and CO (9–8) lines have consistent median values of 54<sup>+6</sup><sub>-6</sub> and 64<sup>+8</sup><sub>-9</sub> km s⁻¹, respectively, as shown in the middle panel of Fig. 11 and Table 3. This is different from what is normally observed in nearby galaxies where the ionized gas (i.e., H α) velocity dispersions are substantially higher than the molecular gas (i.e., CO) velocity dispersions (e.g., Fukui et al. 2009; Epinat et al. 2010). However, because the angular resolution of these two lines are different, the velocity dispersions are measured over different volumes, and in the Milky Way, the velocity dispersion in molecular clouds is proportional to the size and mass of the clouds (e.g., Heyer & Dame 2015). In addition, our velocity dispersions are much smaller than the average velocity dispersion (i.e., 129 ± 10 km s⁻¹) estimated using ALMA [C II] data at a resolution of ~0.′′25 in a sample of 27 z ≥ 6 undisturbed quasar host galaxies (Neeleman et al. 2021).

There is a general trend that the gas velocity dispersion increases with redshift or cosmic time (e.g., Glazebrook 2013).

The predicted velocity dispersion for z = 6 galaxies is about ~80 km s⁻¹ following the empirical relationship for 0 < z < 4 galaxies with measured velocity dispersions from either ionized or molecular gas (Übler et al. 2019). The predicted velocity dispersion is consistent with our measured ones – 36<sup>+11</sup><sub>-13</sub>–80<sup>+8</sup><sub>-6</sub> km s⁻¹ for the [C II] line and 54<sup>+9</sup><sub>-9</sub>–98<sup>+6</sup><sub>-6</sub> km s⁻¹ for the CO (9–8) line. This may indicate that the luminous AGN in the center has little influence on the velocity dispersion of the quasar host galaxy. The physical process to drive and maintain the velocity dispersion over cosmic time is still not clear. But a constant energy input is required to retain the turbulence in the ISM, or it will decay within a few megayears as proposed by theoretical works (e.g., Stone et al. 1998; Mac Low et al. 1998). This energy supply may be associated with either stellar feedback (i.e., winds; expanding H II regions; supernovae) or other instability processes (i.e., clump formation; radial flows; accretion; galaxy interactions). The gas outflow along the line of sight traced by the blueshifted absorption of the OH⁺ line would be very much in line with helping the upkeep of turbulence.

The rotation-to-dispersion ratio $V_{rot}/\sigma$ is a measure of the kinematic support from ordered rotation (i.e., circular motions) versus random motions (i.e., noncircular; turbulence) in a system. We make two kinds of measurements of $V_{rot}/\sigma$ in Table 3 – the ratio between the rotation velocity and velocity dispersion in the flat part of the rotation curve (Col. 10) and the ratio between the maximum rotation velocity and the median velocity dispersion (Col. 11). The ratios are above 4 for both the [C II] and CO (9–8) lines, which are higher than the cutoffs of a rotating system (i.e., a cutoff of 1 from Förster Schreiber et al. 2009; Epinat et al. 2009; a cutoff of 3 from Burkert et al. 2010; Förster Schreiber et al. 2018), indicating that the [C II] and CO (9–8) gas motion is dominated by rotation. These $V_{rot}/\sigma$ values are typical of other z ~ 6 quasar host galaxies traced by resolved [C II] in Neeleman et al. (2021). Similarly, our measured $V_{rot}/\sigma$ ratios are comparable to those of z ~ 2 star-forming galaxies (e.g., Genzel et al. 2008; Förster Schreiber et al. 2009, 2018; Wisnioski et al. 2015), and z ~ 4.5 gravitational lensed dusty star-forming galaxies (e.g., Rizzo et al. 2021).

The inclination angle measured from our kinematic modeling is larger than the one (25°) derived by Tripodi et al. (2022), who also used ³DBAROLO. They constrain the inclination angle <30° considering that $M_{dyn}$ should be larger than $M_{H_2}$, which is derived from the CO (2–1) emission line (Shao et al. 2019) and adopting an $\alpha_{co}$ value of 0.8 $M_{\odot}/(\text{K km s}^{-1} \text{pc}^2)$. As noted in Sect. 4.4, the values of $\alpha_{co}$ may differ from source to source (e.g., Mashian et al. 2015), and based on our dynamical model, we derive an $\alpha_{co}$ value of 0.37 for J2310+1855. The loose constraint on $\alpha_{co}$ (i.e., we adopt a wide value range of 0.2–14; Mashian et al. 2015) allows for a more flexible fitting of the inclination angle in the ³DBAROLO modeling. The initial guess of the inclination angle of 40° in the ³DBAROLO modeling is from the ratio between the photometric minor and major source sizes (which are roughly along the kinematic minor and major axes, respectively) of the CO (9–8) emission assuming an intrinsic round gas disk (before inclination). However, the observed kinematic minor source size is larger than the kinematic major source size (after inclination) for the [C II] line, which indicates that the intrinsic kinematic minor source size of [C II] must be larger than its intrinsic kinematic major source size (before inclination). As we discuss above, the [C II] and CO (9–8) gas are coplanar and corotating in the host galaxy of quasar J2310+1855. Thus, [C II] and CO (9–8) gas are well mixed. Therefore, the intrinsic kinematic minor source size of CO (9–8) might be larger than its intrinsic kinematic major source size. As a result, from
the kinematical aspect, the inclination angle should be larger
than 40°, or conservatively the inclination angle cannot be much
below 40° considering a 5% uncertainty of the initial guess of
the inclination angle.

5.4.2. Outflow traced by OH⁺

The OH⁺ line has been detected in emission, absorption and P-Cygni-shaped profiles in galaxies (e.g., van der Werf et al. 2010; Rangwala et al. 2011; Spinoglio et al. 2012; González-Alfonso et al. 2013; Riechers et al. 2013, 2021a,b; Berta et al. 2021; Butler et al. 2021; Stanley et al. 2021). The production of OH⁺ is mainly driven by H⁺ + O → O⁺ + H followed by O⁺ + H₂ → OH⁺ + H and H⁺ + OH → OH⁺ + H (e.g., van der Werf et al. 2010). The OH⁺ is in a unstable state, and it reacts rapidly with H₂ molecules following the main chemical reactions: OH⁺ + H₂ → H₂O⁺ + H and then H₂O⁺ + H₂ → H₃O⁺ + H (e.g., Gerin et al. 2010). The key initiation for this simple chemical sequence is the production of H⁺, and it is argued that the atomic H ionization is most likely dominated by X-rays and/or cosmic rays (e.g., González-Alfonso et al. 2013). The OH⁺ emission shows double-horned profile with comparable intensities of ~0.30 ± 0.05 mJy of the two peaks. Our measured luminosities of the OH⁺ emission and FIR follow the weak correlation between these two parameters found by Riechers et al. (2021a). At ~384 ± 128 km s⁻¹ (relative to the [C II] redshift from our ALMA Cycle 0 observations; Wang et al. 2013), we detect a ~3σ dip with a peak value of ~0.18 ± 0.05 mJy. This absorption velocity shift is at the high end of the OH⁺ absorption velocity shifts in the range of ~130–360 km s⁻¹ of z = 2–6 starburst galaxies (Riechers et al. 2021a), and is comparable to the widths of the [C II] and CO(9–8) lines of J2310+1855. This P-Cygni profile (i.e., blueshifted absorption and redshifted emission) may indicate the outflow of gas along the line of sight. The outflow signature of broad blue-shifted absorption line of J2310+1855 is recently reported by Bischetti et al. (2022) through VLT/X-shooter observations. The OH⁺ absorption and emission components are co-spatial (left panel of Fig. 3). This is different from, for example, a hyperluminous dusty starbursting major merger ADFS-27 in Riechers et al. (2021b), where the offset is >1 kpc.

In order to further understand the outflow traced by OH⁺
absorption, we estimate the outflow mass with the observed
OH⁺ absorption spectrum and its best-fit Gaussian model. The
optical depth of an unsaturated absorption line can be calcu-
lated with τ<sub>line</sub> = -ln(f<sub>trans</sub>), where f<sub>trans</sub> is the fraction of
the continuum emission that is still transmitted, which can be
expressed as (S<sub>v<sub>cont</sub></sub> + S<sub>v<sub>line</sub></sub>)/S<sub>v<sub>cont</sub></sub>, where S<sub>v<sub>line</sub></sub> is the flux den-
sity of the continuum subtracted OH⁺ absorption, and S<sub>v<sub>cont</sub></sub> is
the flux density of the continuum emission for which we use a
constant value of 1.53 mJy for the whole line. We measure the
peak optical depth τ<sub>OH+</sub> = 0.13 ± 0.04. The velocity-integrated
optical depths τ<sub>OH+</sub> dv are 35 ± 7 and 44 ± 31 km s<sup>-1</sup> by inte-
grating the OH⁺ absorption region and the best-fit Gaussian
region of the OH⁺ absorption, respectively. The OH⁺ column
density N<sub>OH+</sub> is related to the velocity-integrated optical depth
τ<sub>OH+</sub> dv through τ<sub>OH+</sub> dv = λ³ A N<sub>OH+</sub> / 8π, where A is the Einstein
coefficient with a value of 2.11 × 10<sup>-2</sup> s<sup>-1</sup> (e.g., Butler et al.
2021), and λ is the rest frame wavelength of the OH⁺ line (i.e.,
0.029 cm). We thus derive N<sub>OH+</sub> to be (1.7 ± 0.3) × 10<sup>14</sup> and
(2.2 ± 1.5) × 10<sup>14</sup> cm<sup>-2</sup>, respectively. Then assuming a rela-
tive abundance of log<sub>10</sub>(N<sub>OH+</sub>/N<sub>H</sub>) = -7.8 (e.g., Bialy et al.
2019), we measure the neutral hydrogen column density N<sub>H</sub> of

(1.1±0.2)×10²² and (1.4±0.9)×10²² cm⁻², respectively. With a
source radius of 1.3 kpc (the effective radius of CO (9–8) line in
Table 2), we thus derive the neutral hydrogen mass and neutral
gas mass (i.e., with a correction factor of 1.36 for the presence
of helium) to be (4.6±0.9) × 10⁸ and (6.2±1.2) × 10⁸, and
(5.8±4.0)×10⁸ and (7.9±5.5)×10⁸ M☉, respectively. The neu-
tral gas mass of the outflow is 3±1% and 4±3% of the molecular
gas mass in the host galaxy, respectively. Our neutral hydrogen
mass is a few times smaller than the median value of ~2×10⁹ M☉
for a sample of 18 z = 2–6 starburst galaxies that are also traced
by OH⁺ absorption (Riechers et al. 2021a).

Finally, we do not find evidence of gas outflow in the spec-
tra of [C II] or CO (9–8) emission lines in the form of high-
velocity wings (see Figs. 1 and 2). However, Tripodi et al. (2022)
reported the presence of high-velocity components on both the
red and blue sides of the [C II] emission line in J2310+1855.
This discrepancy might be due to differences in the sensitivities
and/or angular resolution of the data used in the analysis. We
have better angular resolution (0.′11 × 0.′09 versus 0.′17 × 0.′15),
although the sensitivities are comparable (0.17 mJy beam⁻¹ per
17 km s⁻¹ versus 0.23 mJy beam⁻¹ per 8.5 km s⁻¹). To under-
stand this discrepancy, we used a nature weighting (robust =
2.0) to TCLEAN our [C II] data, which give more weights on
short baselines. The resulted clean beam size is 0.′16 × 0.′13,
which is similar to the one of Tripodi et al. (2022), and a sen-
sitivity of 0.15 mJy beam⁻¹ per 17 km s⁻¹. The residual spec-
trum between the observed one and the best-fit double-Gaussian
model is only noise-like, with an rms of 0.70 mJy. The reality
of the high-velocity wings reported by Tripodi et al. (2022) is
therefore doubtful. As shown in previous studies, the [C II] emis-
sion line has proven to be a difficult tracer of outflow activity
(e.g., Meyer et al. 2022) and, as shown by Spilker et al. (2020),
strongly lensed dusty star-forming galaxies at z > 4 that show
clear OH outflow do not display high-velocity wings in their
[C II] emission lines. In addition, our CO (9–8) data does not
have enough S/N to look for broad line wings.

5.5. *The dynamics of the quasar-host system*

The dynamical modeling of the rotation curves can yield
detailed knowledge of the composition and distribution of mat-
ter in the quasar-host system. In Sect. 4.4, we decomposed the
J2310+1855 quasar-host system into four dynamical compo-
nents (black hole, stars, gas, and dark matter) and modeled the
potential contribution of each component to the measured circu-
lar velocities of the high-resolution [C II] data cube.

5.5.1. The dynamical structures

The black hole mass is a physical parameter that is critical for
understanding the coevolution of the central SMBHs and their
host galaxies from the local to the high-redshift Universe (e.g.,
Kormendy & Ho 2013; Shen 2013). The mass of the black hole
of the quasar J2310+1855 at z = 6 is on the order of 10⁹ M☉
(Jiang et al. 2016; Feruglio et al. 2018), which is measured using
the width of NIR emission lines and relationships calibrated
with low-redshift AGN. The black hole component is impor-
tant to shape the gravitational potential in the central region of
the host, contributing ~120 km s⁻¹ (for M<sub>BH</sub> = 10⁹ M☉) at a
radius of 0.29 kpc (the innermost position we can reach with
our current resolution), where we measured a circular velocity
of ~250 km s⁻¹. For our target, the value of R<sub>e; star</sub> is strongly
covariant with the black hole mass and the inferred physical
parameters of the stellar component. As shown in Fig. 12 and
Table 4, a factor of 2 change in the black hole mass causes a factor of 4 change in the inferred stellar mass. If we consider an extreme case without a central black hole, the stellar profile needs to become strongly centrally concentrated, with a stellar effective radius of $0.21^{+0.55}_{-0.17}$ kpc to make up for the black hole (as shown in the right panel of Fig. 12 and Col. (c) of Table 4). There is no direct observations to constrain the distribution of the stellar component in this quasar host galaxy. A highly concentrated stellar component might be possible, and thus we cannot definitely remove such an extreme case without a central SMBH. When making the black hole mass a free parameter, we confirm that the black hole mass is on the order of $10^9 M_\odot$; this is the first time that a black hole mass has been dynamically measured at $z \sim 6$. The evidence of the existence of a central SMBH is marginal, and we need more and better data to verify (i.e., James Webb Space Telescope observations to constrain the stellar property, and higher-resolution ALMA [C II] observations). However, we can build a strong case for the existence of a $>10^9 M_\odot$ black hole from our dynamical modeling, unless the stellar distribution in the case without a black hole is strange (but not impossible), as do the other line width metrics (e.g., Jiang et al. 2016; Feruglio et al. 2018). In either scenario (making the black hole mass free or fixed) the gas mass stays stable, as the effective radius of the gas component is fixed based on our ALMA CO (9–8) observations and the gas mass is mainly determined by the circular velocities in the outer regions.

The mass and size of the stellar component can provide essential insights into galaxy evolution (e.g., Hodge et al. 2016). The stellar mass and size correlation of galaxies across cosmic time provides an important view of the assembly history of galaxies (e.g., Kawinwanichakij et al. 2021). However, for quasar host galaxies at $z \sim 6$, it is very difficult to determine these stellar properties through rest frame optical and NIR imaging, as the starlight is overwhelmed by the central quasar. Thus, the comparison with the $z \sim 6$ objects on the stellar mass-size plane is currently not feasible. This also leads to limitations on our dynamical modeling including a stellar component of the quasar-host system in Sect. 4.4. Our modeled stellar mass (around $10^9 M_\odot$) has a large uncertainty. It has been found that not all stellar components are necessarily co-spatial with the gas/dust for SMGs. Thus, the total $M_{star}$ may be under-predicted with the dynamical modeling based on the kinematics of the gas (e.g., Gómez-Guijarro et al. 2018). High-resolution stellar continuum emission observations toward our target and similar luminous quasars in the early Universe by for example the James Webb Space Telescope (JWST) are critical to investigate the $M_{star} - R_{e,star}$ relation and the role of these $z \sim 6$ quasars in galaxy evolution. The stellar effective radius $R_{e,star}$ (in Table 4) is roughly similar to our measured CO and [C II] effective radii, but larger than the dust continuum effective radius, albeit with large uncertainties (in Table 2). This is consistent with the simulated TNG50 star-forming galaxies (Popping et al. 2022). Figure 12 shows dynamical evidence for a massive ($\sim 10^9 M_\odot$) stellar component, especially when we fix the black hole mass to be $1.8 \times 10^9 M_\odot$. In this case, the contribution to $V_c$ from stars rises monotonically to the center of the galaxy, as is typical for nearby massive spiral galaxies (e.g., Simard et al. 2011; Gao et al. 2019, 2020), and such a stellar bulge signature has already been proven to exist using a dynamical method in some dusty star-forming galaxies at $z \sim 4-5$ (e.g., Rizzo et al. 2020, 2021; Tsukui & Iguchi 2021).

The gas component is well constrained by our ALMA CO (9–8) image decomposition and the dynamical modeling (see Sects. 4.1 and 4.4). The gas fraction ($f_{gas} = M_{gas}/M_{bary}$

and $M_{bary} = M_{star} + M_{gas}$) is >70%, and thus it dominates the quasar-host system. The large reservoir of molecular gas ($M_{gas} \sim 10^{10} M_\odot$) is responsible for fueling the bulk of the star formation in the central region of the quasar host galaxy. We should point out that the scale radius of the gas component in the dynamical modeling is from the CO (9–8) intensity map assuming an exponential gas distribution. However, the overall molecular gas (e.g., traced by CO (1–0)) might be more extended than that traced by high-$J$ CO. We adopt a 1.5 times larger gas scale radius to do the dynamical modeling, and get an $\alpha_{CO}$ value of $0.73^{+0.16}_{-0.24}$.

The dark matter component is less constrained. The mass of the dark matter halo is $16^{+26}_{-12}\%$ and $9^{+28}_{-8}\%$ of the total dynamical mass of the quasar-starburst system within a radius of two times the gas half-light radius, when making the black hole mass free and fixed, respectively. The virial velocity and radius of the dark matter halo both have large uncertainties as shown in Table 4. The virial radius of the dark matter halo is much more extended than the gas distribution as shown in Table 2, which is consistent with numerical simulations (e.g., Khandai et al. 2012).

### 5.5.2. The coevolution of the quasar and its host galaxy

To study the coevolution of central black holes and their host galaxies, we typically investigate the evolution of the $M_{BH} - \sigma$ and the $M_{BH} - M_{bulge}$ relations across cosmic time (e.g., Kormendy & Ho 2013), where $\sigma$ and $M_{bulge}$ are the velocity dispersion and total mass of the bulge, respectively.

Considering that the stars and molecular gas trace the same underlying potential, and following Davis et al. (2013), we assume that the second moments of the stars and molecular gas are equal. Thus, we can derive a simple zeroth-order estimate of the velocity dispersion:

$$V_{\text{stars}}^2 + \sigma_{\text{stars}}^2 = V_{\text{gas}}^2 + \sigma_{\text{gas}}^2, \qquad (11)$$

where $V_{gas}$ and $\sigma_{gas}$ are the velocity and velocity dispersion of the molecular gas, for which we use the maximum circular velocity ($\sim 245 \text{ km s}^{-1}$; the maximum rotation velocity corrected by the asymmetric drift) and the median velocity dispersion ($\sim 64 \text{ km s}^{-1}$) of the CO (9–8) line shown in Table 3 as proxies. $V_{stars}$ is the stellar velocity, for which we use the modeled stellar velocity ($40 \text{ and } 99 \text{ km s}^{-1}$ when making the black hole mass free and fixed, respectively) at the half-mass radius of the stellar component (1.76 and 1.39 kpc when making the black hole mass free and fixed, respectively). The predicted stellar velocity dispersions ($\sigma_{stars}$) are 250 and 233 km s$^{-1}$ where we free or fix the black hole mass, respectively. Based on the local relationship for nearby galaxies from Kormendy & Ho (2013), we predict a bulge velocity dispersion of $\sim 335 \pm 20$ and $\sim 300 \pm 16 \text{ km s}^{-1}$ for a central black hole of mass $2.97 \times 10^9$ and $1.8 \times 10^9 M_\odot$ (fixed value in the dynamical modeling), respectively. These bulge velocity dispersions from the local relationship are higher than the predicted stellar velocity dispersions including the scatter in the local relationship, also when we adopt the [C II] kinematic parameters. This may indicate that the black hole evolves faster than the host galaxy of J2310+1855.

We present the baryonic mass $M_{bary}$ from the best-fit dynamical modeling and the ratio of $M_{BH}/M_{bary}$ between the black hole mass and the baryonic mass as a function of radius in Fig. 16. Inside a radius that is two times the half-light radius of the CO emission ($\sim 2.6 \text{ kpc}$), the measured ratio of $M_{BH}/M_{bary}$ is $\sim 28^{+12}_{-9}$.
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Baryonic mass $M_{\text{bary}}$ (red lines) measured from the dynamical modeling in Sect. 4.4 as a function of radius, and the corresponding ratio (green lines) between the black hole mass $M_{\text{BH}}$ and the baryonic mass $M_{\text{bary}}$. The solid and dashed lines represent the scenarios in which the black hole mass is free and fixed during the dynamical modeling, respectively.</figcaption></figure>

and ~$16_{-2}^{+4}$ times higher than the local relationship when making the black hole mass free and fixed in the dynamical modeling, respectively. If we consider a stellar-spheroidal component (i.e., the stellar component in our dynamical modeling in Sect. 4.4), the $M_{\text{BH}}/M_{\text{star}}$ will be >50 times higher than the local relationship (Kormendy & Ho 2013). This may indicate that the central SMBH grows the bulk of its mass before the formation of most of the stellar mass in this quasar host galaxy in the early Universe. Wang et al. (2013) and Tripodi et al. (2022) also found the same result for quasar J2310+1855. This is consistent with what has been found for other high-luminosity $z \sim 6$ quasars (e.g., Wang et al. 2019).

However, as we discussed in Sect. 5.5.1, the black hole mass and the parameters (e.g., the mass, the effective radius, and the Sérsic index) of the stellar component are somewhat degenerate in the dynamical modeling. Thus, the modeled black hole mass, stellar mass and stellar velocity adopted in the discussion above have substantial and coupled uncertainties to the investigation of the quasar-host coevolution. Higher angular resolution and sensitivity observations toward the gas and stellar components are needed to better constrain these parameters with dynamical modeling. In addition, Eq. (11), which we used to predict the stellar velocity dispersion comes with additional uncertainties. Davis et al. (2013) have compared the predicted stellar velocity dispersion using this equation with the measured values. Some data points with stellar velocity dispersion similar to the circular velocity (which is fairly common in elliptical galaxies) appear as obvious outliers, which may indicate that the zeroth order approximations may break down in this situation.

## 6. Summary

We have reported on ALMA high-resolution (sub-kiloparsec-to kiloparsec-scale) observations of the [C II], CO (9–8), and OH$^+$ (1$_{1}$–0$_1$) lines and their underlying dust continuum toward the FIR luminous quasar SDSS J2310+1855 at $z = 6.0031$. We have investigated the ISM distributions, the gas kinematics, and

the quasar-host system dynamics with various models. Below are the main results:

- The [C II] and CO (9–8) lines and their underlying dust continuum are all spatially resolved. The OH$^+$ emission is a point source, and the emission peak is only at a level of ~5σ. The line widths of the [C II] and CO (9–8) emission are consistent, which suggests that these two lines trace the same gravitational potential. We observed asymmetric spectral line profiles for the [C II] and CO (9–8) lines. The enhanced intensity of the lines toward the red part of the spectrum may be due to an asymmetric distribution of gas, or because the temperatures are higher on one side of the galaxy.

- We used 2D elliptical Sérsic models to reproduce the observed intensity maps of the [C II] and CO (9–8) lines and the dust continuum. The extended FIR dust continuum and the CO (9–8) emission can be fitted with a Sérsic distribution with a Sérsic index of ~1 as well as an additional central point component. The point and extended components may represent emission from an AGN dusty molecular torus and the quasar host galaxy, respectively. The derived flux ratio from the AGN dust torus and the star-forming activity for the dust emission are consistent with our rest-frame UV-to-FIR SED decomposition result. The [C II] emission follows a flatter distribution (Sérsic index of 0.59), which may be due to high dust opacity increasing the FIR background of the [C II] line and reducing the [C II] emission at smaller radii. The dust temperature drops with distance from the center. The effective radius of the dust continuum is smaller than that of the line emission and dust mass surface density, but is consistent with that of the SFR surface density. This may indicate that dust emission is a less robust tracer of the dust and gas distribution but is a decent tracer of the obscured star formation activity. The effective radius of the [C II] underlying dust continuum is consistent within ~10% with that of the CO (9–8) underlying dust continuum. The two wavelengths are very close and shortward of the Rayleigh-Jeans tail; therefore, both are dominated by dust at similar temperatures. The half-light radius of the [C II] line is smaller than that of the CO (9–8) line, which may be a reflection of the radial dependence of the gas excitation and the radiation field. A larger [C II]-FIR deficit is observed in the center compared to the outer region, which is likely due to a higher dust temperature and higher dust opacity in the galaxy center. The spatially resolved KS relation ($\Sigma_{\text{SFR}} \sim \Sigma_{M_{\text{gas}}}^{2.50 \pm 0.21}$) with $\Sigma_{M_{\text{gas}}}$ measured from the CO (9–8) line of J2310+1855 is steeper than in the local Universe and $z \sim 2–4$ starburst samples. However, when $\Sigma_{M_{\text{gas}}}$ is measured from the [C II] line, the relation ($\Sigma_{\text{SFR}} \sim \Sigma_{M_{\text{gas}}}^{1.37 \pm 0.07}$) is consistent with the local and low-redshift samples.

- A rotating gas disk is observed in both the [C II] and CO (9–8) lines, as seen from the obvious gradients in the velocity maps and the “S” shaped position-velocity diagrams. We applied 3D tilted ringmodels to both line data cubes with $^3\text{D}\text{BAROLO}$. These two lines show similar kinematical signatures: the gas disk geometry (i.e., consistent inclination angle and position angle of the gas disk) and the gas motions (i.e., rotation dominated). The rotation velocities and velocity dispersions traced by [C II] and CO (9–8) lines are consistent. The circular velocities of [C II] and CO (9–8) lines are in excellent agreement, which is consistent with their identical line width, indicating that they trace the same gravitational potential. We suggest that the [C II] and CO (9–8) gas
are coplanar and corotating in this quasar host galaxy. The velocity dispersions and the $V_{\text{rot}}/\sigma$ ratio of our quasar host galaxy are comparable with other high-$z$ populations, for example $z \ge 6$ quasar hosts and $z \sim 2$ star-forming galaxies. The $\text{OH}^+$ line shows a P-Cygni profile with an absorption at $\sim 400 \text{ km s}^{-1}$, which may indicate an outflow with a neutral gas mass of $(6.2 \pm 1.2) \times 10^8 M_{\odot}$ along the line of sight.

- For the purpose of quantifying the dynamical contributions from different matter components, we decomposed the circular rotation curve measured from the high-resolution [C II] line into four components (black hole, stars, gas, and dark matter). The whole quasar-host system appears baryonic matter dominated. The gas component is well constrained, which reveals a large quantity of gas with a mass on the order of $10^{10} M_{\odot}$. The stellar and dark matter components are not well constrained and depend on whether we allow the black hole mass to be a free parameter. We measured the black hole mass to be $2.97_{-0.77}^{+0.51} \times 10^9 M_{\odot}$. This is the first time that a SMBH mass has been measured dynamically at $z \sim 6$. This result is consistent with that measured from Mg II and C IV lines with the local scaling relations. The stellar component has a dynamical mass measured on the order of $10^9 M_{\odot}$ when the Universe was only $\sim 0.93$ Gyr old. This result is more robust when we fix the black hole mass to be $1.8 \times 10^9 M_{\odot}$ from the quasar spectrum. We note that we may miss the stellar components outside the gas and dust-emitting regions. The relations between the black hole mass and the stellar velocity dispersion and the baryonic mass of this quasar indicate that the central SMBH evolved earlier than its host galaxy.

In the future we plan to improve our understanding of the ISM distribution, gas kinematics, and system dynamics of quasars at the earliest epoch, with higher-resolution observations toward various rest-frame FIR ISM tracers with ALMA and toward the rest-frame optical and NIR stellar emission with the JWST.

**Acknowledgements.** We would like to dedicate this paper to the memory of Dr. Yu Gao. Dr. Gao passed away on 21 May 2022. During his career, Dr. Gao made numerous key contributions to the physics of star formation and the interstellar medium in galaxies; in particular, his work established the star formation law in terms of dense molecular gas content. Dr. Gao left a lasting legacy as a brilliant scientist, as a most respected colleague and caring mentor, and above all, a true gentleman. Dr. Gao will always be in our hearts and thoughts as a truly remarkable human being. We would like to thank the anonymous reviewer for all valuable comments and suggestions, which helped us to improve the quality of the manuscript. This work makes use of the following ALMA data: ADS/JAO.ALMA# 2018.1.00597.S; 2015.1.01265.S; 2015.1.00997.S and 2011.0.00206.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. R.W. acknowledges supports from the National Science Foundation of China (NSFC) grants No. 12173002, 11991052, 11721303, and 11533001, and the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A06. Y.S. acknowledges discussions with Francesca Rizzo, Gang Wu, Aiyuan Yang and Sudeep Neupane.

## References

Akins, H. B., Fujimoto, S., Finlator, K., et al. 2022, ApJ, 934, 64

Alonso-Herrero, A., Pereira-Santaella, M., García-Burillo, S., et al. 2018, ApJ, 859, 144

Amorín, R., Muñoz-Tuñón, C., Aguerri, J. A. L., & Planesas, P. 2016, A&A, 588, A23

Banerji, M., Jones, G. C., Carniani, S., DeGraf, C., & Wagg, J. 2021, MNRAS, 503, 5583

Begeman, K. G., Broeils, A. H., & Sanders, R. H. 1991, MNRAS, 249, 523

Berta, S., Young, A. J., Cox, P., et al. 2021, A&A, 646, A122

Bertoldi, F., Carilli, C. L., Cox, P., et al. 2003a, A&A, 406, L55

Bertoldi, F., Cox, P., Neri, R., et al. 2003b, A&A, 409, L47

Bewketu Belete, A., Andreani, P., Fernández-Ontiveros, J. A., et al. 2021, A&A, 654, A24

Bialy, S., Neufeld, D., Wolfire, M., Sternberg, A., & Burkhart, B. 2019, ApJ, 885, 109

Bianchi, S., & Xilouris, E. M. 2011, A&A, 531, L11

Bigiel, F., Leroy, A., Walter, F., et al. 2010, AJ, 140, 1194

Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, N.J: Princeton University Press)

Bischetti, M., Feruglio, C., D'Odorico, V., et al. 2022, Nature, 605, 244

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Bok, J., Blyth, S. L., Gilbank, D. G., & Elson, E. C. 2019, MNRAS, 484, 582

Bothwell, M. S., Chapman, S. C., Tacconi, L., et al. 2010, MNRAS, 405, 219

Bournaud, F., Combes, F., Jog, C. J., & Puerari, I. 2005, A&A, 438, 507

Burkert, A., Genzel, R., Bouché, N., et al. 2010, ApJ, 725, 2324

Butler, K. M., van der Werf, P. P., Rybak, M., et al. 2021, ApJ, 919, 5

Calistro Rivera, G., Hodge, J. A., Smail, I., et al. 2018, ApJ, 863, 56

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455

Carignan, C., Chemin, L., Huchtmeier, W. K., & Lockman, F. J. 2006, ApJ, 641, L109

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Carilli, C. L., Neri, R., Wang, R., et al. 2007, ApJ, 666, L9

Carilli, C. L., Daddi, E., Riechers, D., et al. 2010, ApJ, 714, 1407

Carniani, S., Gallerani, S., Vallini, L., et al. 2019, MNRAS, 489, 3939

Casasola, V., Cassarà, L. P., Bianchi, S., et al. 2017, A&A, 605, A18

Chen, C.-C., Hodge, J. A., Smail, I., et al. 2017, ApJ, 846, 108

Ciotti, L., & Bertin, G. 1999, A&A, 352, 447

Cochrane, R. K., Hayward, C. C., Anglés-Alcázar, D., et al. 2019, MNRAS, 488, 1779

Combes, F., García-Burillo, S., Audibert, A., et al. 2019, A&A, 623, A79

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Davis, T. A., Alatalo, K., Bureau, M., et al. 2013, MNRAS, 429, 534

de Blok, W. J. G., & Bosma, A. 2002, A&A, 385, 816

de Blok, W. J. G., & Walter, F. 2014, AJ, 147, 96

de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2648

de Blok, W. J. G., Walter, F., Smith, J. D. T., et al. 2016, AJ, 152, 51

Decarli, R., Walter, F., Venemans, B. P., et al. 2017, Nature, 545, 457

Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97

Di Mascia, F., Gallerani, S., Behrens, C., et al. 2021, MNRAS, 503, 2349

Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021

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

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

Draine, B. T. 2003, ARA&A, 41, 241

Epinat, B., Amram, P., Balkowski, C., & Marcelin, M. 2010, MNRAS, 401, 2113

Epinat, B., Contini, T., Le Fèvre, O., et al. 2009, A&A, 504, 789

Esparza-Redondo, D., Gonzalez-Martín, O., Dultzin, D., et al. 2021, A&A, 651, A91

Evans, N. J. I., Heiderman, A., & Vutisalchavakul, N. 2014, ApJ, 782, 114

Fan, X., White, R. L., Davis, M., et al. 2000, AJ, 120, 1167

Ferrarese, L., & Ford, H. 2005, Space Sci. Rev., 116, 523

Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9

Feruglio, C., Fiore, F., Carniani, S., et al. 2018, A&A, 619, A39

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21

Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364

Frank, B. S., de Blok, W. J. G., Walter, F., Leroy, A., & Carignan, C. 2016, AJ, 151, 94

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

Fukui, Y., Kawamura, A., Wong, T., et al. 2009, ApJ, 705, 144

Gallerani, S., Ferrara, A., Neri, R., & Maiolino, R. 2014, MNRAS, 445, 2848

Gao, H., Ho, L. C., Barth, A. J., & Li, Z.-Y. 2018, ApJ, 862, 100

Gao, H., Ho, L. C., Barth, A. J., & Li, Z.-Y. 2019, ApJS, 244, 34

Gao, H., Ho, L. C., Barth, A. J., & Li, Z.-Y. 2020, ApJS, 247, 20

García-Burillo, S., Combes, F., Usero, A., et al. 2014, A&A, 567, A125

García-Burillo, S., Alonso-Herrero, A., Ramos Almeida, C., et al. 2021, A&A, 652, A98

Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

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

Genzel, R., Tacconi, L. J., Kurk, J., et al. 2013, ApJ, 773, 68

Gerin, M., de Luca, M., Black, J., et al. 2010, A&A, 518, L110

Glazebrook, K. 2013, PASA, 30, e056a

Gómez-Guijarro, C., Toft, S., Karim, A., et al. 2018, ApJ, 856, 121

González-Alfonso, E., Fischer, J., Isaak, K., et al. 2010, A&A, 518, L43

González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2013, A&A, 550, A25
González-López, J., Barrientos, L. F., Gladders, M. D., et al. 2017, ApJ, 846, L22

Graham, A. W., Erwin, P., Caon, N., & Trujillo, I. 2001, ApJ, 563, L11

Gullberg, B., De Breuck, C., Vieira, J. D., et al. 2015, MNRAS, 449, 2883

Gullberg, B., Swinbank, A. M., Smail, I., et al. 2018, ApJ, 859, 12

Gullberg, B., Smail, I., Swinbank, A. M., et al. 2019, MNRAS, 490, 4956

Gunn, J. E., Gott, J., & Richard, I. 1972, ApJ, 176, 1

Haas, M., Lemke, D., Stickel, M., et al. 1998, A&A, 338, L33

Hashimoto, T., Inoue, A. K., Tamura, Y., et al. 2019a, PASJ, 71, 109

Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019b, PASJ, 71, 71

Hatsukade, B., Tamura, Y., Iono, D., et al. 2015, PASJ, 67, 93

Heiderman, A., Evans, N. J. I., Allen, L. E., Huard, T., & Heyer, M. 2010, ApJ, 723, 1019

Heyer, M., & Dame, T. M. 2015, ARA&A, 53, 583

Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103

Hönig, S. F. 2019, ApJ, 884, 171

Hönig, S. F., & Kishimoto, M. 2017, ApJ, 838, L20

Ikarashi, S., Ivison, R. J., Caputi, K. I., et al. 2015, ApJ, 810, 133

Iorio, G., Fraternali, F., Nipoti, C., et al. 2017, MNRAS, 466, 4159

Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210

Izumi, T., Onoue, M., Shirakata, H., et al. 2018, PASJ, 70, 36

Izumi, T., Onoue, M., Matsuoka, Y., et al. 2018, ApJ, 908, 235

Jiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222

Jones, G. C., Carilli, C. L., Shao, Y., et al. 2017, ApJ, 850, 180

Jones, G. C., Vergani, D., Romano, M., et al. 2021, MNRAS, 507, 3540

Kawinwanichakij, L., Silverman, J. D., Ding, X., et al. 2021, ApJ, 921, 38

Kenney, J. D. P., Abramson, A., & Bravo-Alfaro, H. 2015, AJ, 150, 59

Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

Kennicutt, R. C., Jr., Calzetti, D., Walter, F., et al. 2007, ApJ, 671, 333

Khandai, N., Feng, Y., DeGraf, C., Di Matteo, T., & Croft, R. A. C. 2012, MNRAS, 423, 2397

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

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

Leipski, C., Meisenheimer, K., Walter, F., et al. 2013, ApJ, 772, 103

Leipski, C., Meisenheimer, K., Walter, F., et al. 2014, ApJ, 785, 154

Lelli, F., Verheijen, M., Fraternali, F., & Sancisi, R. 2012, A&A, 537, A72

Li, Q., Wang, R., Fan, X., et al. 2020a, ApJ, 900, 12

Li, J., Wang, R., Riechers, D., et al. 2020b, ApJ, 889, 162

Li, J., Wang, R., Cox, P., et al. 2020c, ApJ, 900, 131

Liu, L., Weiß, A., Perez-Beaupuits, J. P., et al. 2017, ApJ, 846, 5

Lucero, D. M., Carignan, C., Elson, E. C., et al. 2015, MNRAS, 450, 3935

Luhman, M. L., Satyapal, S., Fischer, J., et al. 1998, ApJ, 504, L11

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

Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754

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

Maiolino, R., Cox, P., Caselli, P., et al. 2005, A&A, 440, L51

Marconi, A., & Hunt, L. K. 2003, ApJ, 589, L21

Mashian, N., Sturm, E., Sternberg, A., et al. 2015, ApJ, 802, 81

Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2022, ApJS, 259, 18

McGreer, I. D., Fan, X., Strauss, M. A., et al. 2014, AJ, 148, 73

Meyer, R. A., Walter, F., Cicone, C., et al. 2022, ApJ, 927, 152

Miller, T. B., Chapman, S. C., Hayward, C. C., et al. 2020, ApJ, 889, 98

Miura, R. E., Kohno, K., Tosaki, T., et al. 2014, ApJ, 788, 167

Molina, J., Wang, R., Shangguan, J., et al. 2021, ApJ, 908, 231

Muñoz, J. A., & Oh, S. P. 2016, MNRAS, 463, 2085

Narayanan, D., & Krumholz, M. R. 2017, MNRAS, 467, 50

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141

Neri, R., Downes, D., Cox, P., & Walter, F. 2014, A&A, 562, A35

Novak, M., Venemans, B. P., Walter, F., et al. 2020, ApJ, 904, 131

Omont, A., Willott, C. J., Beelen, A., et al. 2013, A&A, 552, A43

Pensabene, A., Decarli, R., Bañados, E., et al. 2021, A&A, 652, A66

Planck Collaboration XIII. 2016, A&A, 594, A13

Popping, G., Pillepich, A., Calistro Rivera, G., et al. 2022, MNRAS, 510, 3321

Rangwala, N., Maloney, P. R., Glenn, J., et al. 2011, ApJ, 743, 94

Reynolds, T. N., Westmeier, T., Staveley-Smith, L., Chauhan, G., & Lagos, C. D. P. 2020, MNRAS, 493, 5089

Riechers, D. A., Walter, F., Bertoldi, F., et al. 2009, ApJ, 703, 1338

Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329

Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84

Riechers, D. A., Cooray, A., Pérez-Fournon, I., & Neri, R. 2021a, ApJ, 913, 141

Riechers, D. A., Nayyeri, H., Burgarella, D., et al. 2021b, ApJ, 907, 62

Rizzo, F., Vegetti, S., Powell, D., et al. 2020, Nature, 584, 201

Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS, 507, 3952

Roberts, M. S. 1978, AJ, 83, 1026

Rubin, V. C. 1983, Science, 220, 1339

Rybak, M., Calistro Rivera, G., Hodge, J. A., et al. 2019, ApJ, 876, 112

Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41

Shao, Y., Wang, R., Jones, G. C., et al. 2017, ApJ, 845, 138

Shao, Y., Wang, R., Carilli, C. L., et al. 2019, ApJ, 876, 99

Sharon, C. E., Tagore, A. S., Baker, A. J., et al. 2019, ApJ, 879, 52

Shen, Y. 2013, Bull. Astron. Soc. India, 41, 61

Simard, L., Mendel, J. T., Patton, D. R., Ellison, S. L., & McConnachie, A. W. 2011, ApJS, 196, 11

Simon, J. D., Bolatto, A. D., Leroy, A., Blitz, L., & Gates, E. L. 2005, ApJ, 621, 757

Sofue, Y., Honma, M., & Omodaka, T. 2009, PASJ, 61, 227

Spilker, J. S., Phadke, K. A., Aravena, M., et al. 2020, ApJ, 905, 85

Spinoglio, L., Pereira-Santaella, M., Busquet, G., et al. 2012, ApJ, 758, 108

Stacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373, 423

Stacey, G. J., Hailey-Dunsheath, S., Ferkinhoff, C., et al. 2010, ApJ, 724, 957

Stanley, F., Knudsen, K. K., Aalto, S., et al. 2021, A&A, 646, A178

Stewart, I. M., Blyth, S. L., & de Blok, W. J. G. 2014, A&A, 567, A61

Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, ApJ, 508, L99

Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJ, 842, L15

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

Tadaki, K., Iono, D., Yun, M. S., et al. 2018, Nature, 560, 613

Tadaki, K.-I., Belli, S., Burkert, A., et al. 2020, ApJ, 901, 74

Terzić, B., & Graham, A. W. 2005, MNRAS, 362, 197

Thomson, A. P., Ivison, R. J., Owen, F. N., et al. 2015, MNRAS, 448, 1874

Tripodi, R., Feruglio, C., Fiore, F., et al. 2022, A&A, 665, A107

Tsukui, T., & Iguchi, S. 2021, Science, 372, 1201

Übler, H., Genzel, R., Tacconi, L. J., et al. 2018, ApJ, 854, L24

Übler, H., Genzel, R., Wisnioski, E., et al. 2019, ApJ, 880, 48

Usero, A., García-Burillo, S., Fuente, A., Martín-Pintado, J., & Rodríguez-Fernández, N. J. 2004, A&A, 419, 897

Valluri, M., Merritt, D., & Emsellem, E. 2004, ApJ, 602, 66

van de Hulst, H. C., Raimond, E., & van Werden, H. 1957, Bull. Astron. Inst. Neth., 14, 1

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42

van Eymeren, J., Jütte, E., Jog, C. J., Stein, Y., & Dettmar, R. J. 2011, A&A, 530, A30

Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37

Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJ, 837, 146

Venemans, B. P., Decarli, R., Walter, F., et al. 2018, ApJ, 866, 159

Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130

Villanueva, V., Ibar, E., Hughes, T. M., et al. 2017, MNRAS, 470, 3775

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

Walter, F., Bertoldi, F., Carilli, C., et al. 2003, Nature, 424, 406

Walter, F., Carilli, C., Bertoldi, F., et al. 2004, ApJ, 615, L17

Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563

Walter, F., Riechers, D., Cox, P., et al. 2009, Nature, 457, 699

Walter, F., Neeleman, M., Decarli, R., et al. 2022, ApJ, 927, 21

Wang, R., Carilli, C. L., Beelen, A., et al. 2007, AJ, 134, 617

Wang, R., Carilli, C. L., Wagg, J., et al. 2008, ApJ, 687, 848

Wang, R., Carilli, C. L., Neri, R., et al. 2010, ApJ, 714, 699

Wang, R., Wagg, J., Carilli, C. L., et al. 2011a, AJ, 142, 101

Wang, R., Wagg, J., Carilli, C. L., et al. 2011b, ApJ, 739, L34

Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44

Wang, R., Shao, Y., Carilli, C. L., et al. 2019, ApJ, 887, 40

Warmels, R., Biggs, A., Cortes, P. A., et al. 2018, ALMA Technical Handbook, ALMA Doc. 6.3, ver. 1.0

Watts, A. B., Catinella, B., Cortese, L., & Power, C. 2020, MNRAS, 492, 3672

Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955

Weiß, A., Requena-Torres, M. A., Güsten, R., et al. 2010, A&A, 521, L1

Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209

Wong, T., Blitz, L., & Bosma, A. 2004, ApJ, 605, 183

Yang, C., Gao, Y., Omont, A., et al. 2013, ApJ, 771, L24

Yue, M., Yang, J., Fan, X., et al. 2021, ApJ, 917, 99

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976
# Appendix A: The [C II] channel maps

We present the [C II] and CO (9–8) channel maps of J2310+1855 in Figs. A.1 and A.2. The [C II] emission peak moves in a circular, counterclockwise path from 219 to -281 km/s, a signature of a rotating gas disk.

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

**Fig. A.1.** [C II] channel maps of J2310+1855. The velocity takes the [C II] redshift from Wang et al. (2013) as a reference. The channel width is ~ 17 km/s (corresponding to 15.625 MHz). The contour levels are [-2, 2, 4, 6, 8, 10, 12, 14] × rms, where rms = 0.17 mJy/beam. The white plus sign represents the HST quasar position. The [C II] emission peak moves in a circular, counterclockwise path from 219 to -281 km/s.
<figure><img src="image_19.png" /><figcaption>Fig. A.2. CO (9–8) channel maps of J2310+1855. The velocity takes the [C II] redshift from Wang et al. (2013) as a reference. The channel width is ~ 32 km/s (corresponding to 15.625 MHz). The contour levels are [-2, 2, 4, 6, 8, 10, 12, 14] × rms, where rms = 0.10 mJy/beam. The white plus sign represents the HST quasar position.</figcaption></figure>

## Appendix B: Two-dimensional elliptical Sérsic function

The 1D Sérsic profile as a function of radius, $r$, is given by Sérsic (1963):

$$ \Sigma(r) = \Sigma_0 \exp^{-(r/r_s)^{1/n}}, \qquad (B.1) $$

where $\Sigma_0$ is the central surface brightness, $r_s$ is the scale length, which is the radius at which the surface brightness $\Sigma(r)$ drops by $e^{-1}$, and $n$ is the shape parameter that controls the degree of curvature of the profile, the so-called Sérsic index. When $n$ is large (i.e., > 4), $r_s$ is too small to measure; therefore, one often redefines the profile in terms of the half-light radius, $r_e$, also known as the effective radius:

$$ \Sigma(r) = \Sigma_0 \exp^{-k(r/r_e)^{1/n}}, \qquad (B.2) $$

where the dependent variable $k$ is coupled to $n$, and $r_e = k^n r_s$. $k$ satisfies $\frac{1}{2}\Gamma(2n) = \gamma(2n; k)$, where $\Gamma$ and $\gamma$ are the Gamma function and lower incomplete Gamma function, respectively. $k$ can be approximated as $2n - \frac{1}{3} + \frac{4}{405n} + \frac{46}{25515n^2} + \frac{131}{1148175n^3} + \frac{2194697}{30690717750n^4}$ for $n > 0.36$ (Ciotti & Bertin 1999).

The 2D elliptical Sérsic function can be expressed as

$$ f(x, y) = A \exp^{-[a(x-x_0)^2+2b(x-x_0)(y-y_0)+c(y-y_0)^2]^{1/(2n)}} \qquad (B.3) $$

$$ a = \frac{\cos^2 \theta}{h_x^2} + \frac{\sin^2 \theta}{h_y^2} \qquad (B.4) $$

$$ b = \frac{\sin 2\theta}{2h_x^2} - \frac{\sin 2\theta}{2h_y^2} \qquad (B.5) $$

$$ c = \frac{\sin^2 \theta}{h_x^2} + \frac{\cos^2 \theta}{h_y^2}, \qquad (B.6) $$

where the coefficient $A$ is the height of the peak, $(x_0, y_0)$ is the center, $h_x, h_y$ are the scale heights along the x and y axis before rotating by $\theta$, and $\theta$ is the rotation angle in radians, which increases counterclockwise counting from the positive direction of the x axis, and $n$ is the Sérsic index. We show surface plots, projected filled contour plots and surface brightness distributions in Fig. B.1 for 2D elliptical Sérsic images with different rotation angles and Sérsic indices.

We use a 2D elliptical Sérsic model to fit the CO (9–8) line and the [C II] and CO (9–8) underlying dust continuum. And as motivated by the existence of AGN torus, we also fit these intensity maps with an additional unresolved nuclear component to represent the dusty and molecular torus. The results are shown in Figs. 7, B.2, and B.3.
<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** Surface plots, projected filled contour plots, and surface brightness distributions for 2D elliptical Sérsic images with $(x_0, y_0) = (0, 0)$, $h_x = 4$, and $h_y = 2$. $\theta$ and $n$ are $\pi/3$ and $0.5$ (top left), $\pi/2$ and $1$ (top right), $2\pi/3$ and $2$ (bottom left). The surface brightness as a function of radius shown in the bottom-right panel is measured using elliptical annuli with the same rotation angle of these images. The radius is the one along the major axis. The black, red and blue dots are measured from the top-left, top-right, and bottom-left images, respectively. The black, red, and blue lines are the fits to these surface brightness distributions with 1D Sérsic function in Eq. B.1. The best-fitted Sérsic indices are labeled in the same colors as the fitting lines.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. B.2.** As in Fig. 7 but for the [C II] underlying dust continuum. The peak values and rms within the dashed gray square are 0.101 and 0.028, and 0.099 and 0.028 Jy/beam km/s for the Sérsic and P+Sérsic residual maps, respectively. The shape of the synthesized beam with a FWHM size of 0.113 × 0.092 is plotted as white ellipses. We measured the [C II] underlying dust continuum luminosity surface density using elliptical rings with the ring width along the major axis of half (0.05) that of the major axis of the [C II] clean beam size, the rotation angle equal to PA (=199°) and the ratio of semiminor and semimajor axis – *b*/*a* of cos(î) (î̃ = 42°), where $\overline{PA}$ and $\overline{i}$ come from the [C II] line kinematic modeling (listed in Table 3).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. B.3.** As in Fig. 7 but for the CO (9–8) underlying dust continuum. The peak values and rms within the dashed gray square are 0.043 and 0.014, and 0.041 and 0.014 Jy/beam km/s for the Sérsic and P+Sérsic residual maps, respectively. The shape of the synthesized beam with a FWHM size of 1.92 × 1.56 is plotted as a white ellipse in each panel.</figcaption></figure>
## Appendix C: Tilted ring model

In this section, we illustrate the parameters used in the tilted ring model. We only observe the line-of-sight velocities as a function of position on the sky, which cannot be directly used to investigate the mass distributions in galaxies. We need to measure the rotation velocity $V_{rot}$ as a function of distance $R$ to the center (the rotation curve). If we assume that the gas is rotating symmetrically in a disk, a number of tilted ring model parameters can define the observed velocity field of a galaxy. The tilted ring model (see Fig. C.1) consists of a set of concentric rings characterized by various geometrical and kinematical components. The geometrical components are:

* $x_0, y_0$: sky coordinates of the rotation center of the galaxy.

* $i(R)$: the inclination (i.e., the angle between the normal to the plane of the galaxy and the line-of-sight.

* $\phi(R)$: the position angle of the major axis of a ring projected onto the sky (i.e., an ellipse). This is an angle taken in a counterclockwise direction from the north of the sky to the major axis of the receding half of the galaxy.

* $\zeta(R)$: the azimuthal angle measured between the major axis on the sky plane and the line connecting the galaxy center $(x_0, y_0)$ and the data point $(X, Y)$ in the galaxy plane. It is related to $i(R), \phi(R), (x_0, y_0)$.

The kinematical components are:

* $V_{sys}$: the velocity of the center of the galaxy with respect to the Sun, the so-called systemic velocity.

* $V_{rot}(R)$: the rotation velocity at distance $R$ from the center. This velocity is assumed to be constant in a ring.

The $V_{rot}$ is perpendicular to the ring and has components in the X and Y directions. $V_X (= V_{rot} \sin \zeta)$ is parallel to the xy plane and has no component in the z direction. However, $V_Y (= V_{rot} \cos \zeta)$ has a component in the z direction: $V_z = V_Y \sin i$. This is the radial velocity ($V_{los}$) that we observe at a point $(x, y)$. It must be corrected for the systemic (line-of-sight) velocity $V_{sys}$. Thus,

$$V_{los}(x, y) = V_{sys} + V_z = V_{sys} + V_{rot} \cos \zeta \sin i \quad (C.1)$$

$$\cos \zeta = \frac{X}{R} = \frac{-(x - x_0) \sin \phi + (y - y_0) \cos \phi}{R} \quad (C.2)$$

$$\sin \zeta = \frac{Y}{R} = \frac{-(x - x_0) \cos \phi - (y - y_0) \sin \phi}{R \cos i} \quad (C.3)$$

Here, we also put the $^3D\text{BAROLO}$ modeling result of CO (9–8) line in Fig. C.2.

<figure><img src="image_23.png" /><figcaption>**Fig. C.1.** Simple schematic diagram for the tilted ring model. The pink semiellipse represents the galaxy plane with a coordinate system of XYZ. The red-based color lines are on or perpendicular to the galaxy plane. The light-blue ellipse presents the sky plane with a coordinate system of xyz. The blue color lines are on or perpendicular to the sky plane. We take the x direction as the east and the y direction as the north. The major axis on the sky plane coincides with the X axis on the galaxy plane. The line-of-sight direction is along the z direction from the bottom to the top. The galaxy plane is inclined with an angle of $i$ (the inclination angle, i.e., the angle between the line of sight and the perpendicular line to the galaxy plane XY). The data point $(X, Y)$ with a distance of $R$ to the center rotates counterclockwise with a rotation velocity of $V_{rot}$ on the galaxy plane. What we can see is the line-of-sight component $V_z$ of the rotation velocity $V_{rot}$. The angle between the major axis on the sky plane and the line connecting the galaxy center $(x_0, y_0)$ and the data point $(X, Y)$ is the azimuth angle $\zeta$. The position angle $\phi$ is defined as the angle counting from the north (y axis) to the major axis of the receding half on the galaxy plane. After tilting, the data point $(X, Y)$, the distance $R$ and the azimuth angle $\zeta$ on the galaxy plane project to $(x, y)$, $R'$ and $\beta$ on the sky plane.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. C.2.** Similar to that of Fig. 10, but for the CO (9–8) line kinematic modeling. **Upper panels** – Shape of the synthesized beam with a FWHM size of 0.187 × 0.153 is plotted in the bottom-left corner of each panel. **Lower panels** – Left and middle panels: Contour levels are [−2, 2, 4, 8, 16] × 0.093 mJy/beam for both the data and the model. Right panel: The spectral resolution is 31.25 MHz, corresponding to 64 km/s.</figcaption></figure>