© The Authors 2022

1. Introduction

Almost 400 quasars at z ≥ 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 ∼ 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¹¹ to 10¹³ L_⊙ and star formation at rates of hundreds to thousands M_⊙ yr⁻¹ in these young quasar host galaxies. A few of these z ∼ 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¹⁰ M_⊙ 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 ∼ 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 ∼ 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 ∼ 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 ≤ 3 (z ≥ 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[[INLINE38]] are in good agreement in a z = 1.4 star-forming galaxy. However, Simon et al. (2005) find that the H[[INLINE40]] 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 ≤ z ≤ 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₁–0₁) 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 Å = 19.30 mag in the Sloan Digital Sky Survey (SDSS) imaging data. The black hole mass has been measured to be (4.17 ± 1.02)×10⁹ M_⊙ and (3.92 ± 0.48)×10⁹ M_⊙ 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 ± 0.5)×10⁹ M_⊙ 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_⊙ yr⁻¹ 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₂O (2₀₂–1₁₁) and para-H₂O (3₂₂-3₁₃), OH⁺ (1₁–0₁) 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 ³P₂-³P₁), 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 (; 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 × 10¹⁰ M_⊙ 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 × 10¹⁰ M_⊙ 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₁–0₁) 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 ΛCDM cosmology with H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_M = 0.3089, and Ω_Λ = 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 (ν_rest = 1900.5369 GHz; redshifted to ν_obs = 271.3851 GHz), along with band-4 observations of the CO (9–8) line (ν_rest = 1036.9124 GHz; redshifted to ν_obs = 148.0648 GHz) and the OH⁺ (1₁–0₁) line (ν_rest = 1033.0582 GHz; redshifted to ν_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¹). The original channel width is 15.625 MHz, corresponding to ∼17 and 32 km s⁻¹ 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₁–0₁) 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σ 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σ. 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 × and × , 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⁻¹ for the [C II] and CO (9–8) lines, respectively. The root mean 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.

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 × 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−1 km s−1. 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−1. 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.

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 resolution (8.83 ± 0.44 Jy km s⁻¹; 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⁺ lines and the possible detection of para-H₂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⁻¹ is seen as in the [C II] line.

Fig. 2. Similar to Fig. 1, but for the CO (9–8) line. The FWHM size of the restoring beam is × . The contour levels are [−3, 3, 6, 9, 12, 15, 18, 21] × rmsco, where rmsco = 0.020 Jy beam−1 km s−1. The spectral resolution is 15.625 MHz, corresponding to 32 km s−1. The rms of the spectrum is 0.31 mJy. The CO (9–8) spectral profile is also asymmetric, with enhancement on the red side.

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⁻¹ 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 resolution (1.31 ± 0.06 Jy km s⁻¹; 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⁺ velocity-integrated intensity map and spectrum are shown in Fig. 3. The peak flux density is 0.066 ± 0.012 Jy beam⁻¹ km s⁻¹ (> 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⁻¹. 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⁻¹ from the Gaussian fit) of the OH⁺ spectrum may indicate an outflow.

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 × is plotted in the bottom-left corner. The contour levels in black are [3, 4, 5] × rmsOH+_em, where rmsOH+_em = 0.012 Jy beam−1 km s−1. The peak flux density is 0.066 ± 0.012 Jy beam−1 km s−1. The overplotted dashed white contours are from the OH+ absorption line, with contour levels of [−3, −2] × rmsOH+_ab, where rmsOH+_ab = 0.008 Jy beam−1 km s−1. 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−1. 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.

Our new ALMA band-4 observations serendipitously partially cover the para-H₂O (1₁₁–0₀₀) line at ν_rest = 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₂O (2₀₂–1₁₁) line from the ALMA archive. The para-H₂O (1₁₁–0₀₀) 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_dust ∼ 20–30 K) and widespread diffuse ISM component gives rise to the ground-state para-H₂O line (Weiß et al. 2010; Liu et al. 2017). Liu et al. (2017) find the ratio of para-H₂O (1₁₁–0₀₀) and para-H₂O (2₀₂–1₁₁) to be < 1 in 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.

Fig. 4. Spectra of the para-H2O (111–000) line (black histogram) and the scaled para-H2O (202–111) 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−1. The observations may only cover the red tail of the para-H2O (111–000) emission line.

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.

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 × and × , respectively. The black contours are [−3, 3, 9, 27, 81] × rmscon for both maps. The rmscon at 262 and 147 GHz are 0.02 and 0.01 mJy beam−1, 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.

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 data pixel-by-pixel using Markov chain Monte Carlo (MCMC) method with the emcee² 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).

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 is plotted in the bottom-left corner of each panel. We measured the residual rms to be 0.037 Jy beam−1 km s−1 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 () that of the [C II] clean beam size, the rotation angle equal to (=199°), and the ratio of semiminor and semimajor axis – b/a of (), where 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 107 L⊙ kpc−2.

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−1 km s−1 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 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 () that of the CO (9–8) clean beam size, the rotation angle equal to (=198°), and the ratio of semiminor and semimajor axis – b/a of (), where and 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 106 L⊙ kpc−2.

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):

(1)

and

(2)

where S_ν is the observed flux density for a target at redshift z. B_ν(T_dust) and B_ν(T_CMB) are the black-body functions with dust temperature T_dust and the cosmic microwave background (CMB) temperature T_CMB, respectively. Ω_app is the apparent solid angle. In principle, Ω_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 ∼ 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. τ_ν is the optical depth, which is a function of frequency ν (rest-frame) and dust mass surface density . M_{dust,  app} is the apparent dust mass. D_A is the angular distance. β is the emissivity index. We here adopted the absorption coefficient per unit dust mass κ₀ of 13.9 cm² g⁻¹ at the reference frequency ν_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_dust of K, a total dust mass M_{dust,  app} of and an overall emissivity index β of . 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 and , 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 yr⁻¹ calculated from the above-mentioned IR luminosity and the formula in Kennicutt (1998).

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(−m2/2) (m = 0.5, 1, 1.5, 2) of the volume. The vertical dashed lines show 16th, 50th (median), and 84th percentiles.

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 . 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 β 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_dust ∝ r^{−0.41 ± 0.20}. The radial dependence of Σ_{dust_mass} and Σ_SFR are consistent with an exponential distribution Σ(r) = Σ₀exp^−r/r_s, where r_s is the exponential scale length and Σ₀ is the central surface density. The best-fit relation for Σ_{dust_mass} is shown in the top-right panel of Fig. 9 with r_{s,  dust_mass} = 0.77 ± 0.27 kpc, corresponding to an effective radius of 1.29 ± 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 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 ± 0.69)×10⁹ M_⊙, which is consistent with the value measured from the SED fitting. An exponential fit for Σ_SFR is shown in the bottom-right panel of Fig. 9 with r_s,  SFR = 0.38 ± 0.09 kpc, corresponding to an effective radius of 0.64 ± 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 yr⁻¹, 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., Σ_{dust_mass} > 10⁸ M_⊙ kpc⁻²), and the optical depth for the dust continuum emission at both frequencies is above 0.1 inside ∼1 kpc.

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.

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 (^3DBAROLO³; Di Teodoro & Fraternali 2015).

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 3DBAROLO, 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 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 3DBAROLO 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 3DBAROLO 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 3DBAROLO 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.

The ^3DBAROLO 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 ϕ; the rotation velocity V_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⁴). The ^3DBAROLO 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 ring, ^3DBAROLO 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 ^3DBAROLO 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_ring as half of the clean beam size (i.e., Nyquist sampling): and for the modeling of the [C II] and CO (9–8) lines, respectively. The central radius of each ring is W_ring × 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 kinematic 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_size): i = cos⁻¹(r_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.

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.

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 (Φ) and can be measured by correcting the rotation curve for the asymmetric drift:

(3)

where V_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 ^3DBAROLO, and Φ 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

(4)

where V_BH, V_star, V_gas and V_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_BH is

(5)

where G and M_BH are the gravitational constant and the black hole mass, respectively. We measured a circular velocity of km s⁻¹ 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 . The black hole mass for our target J2310+1855 measured from the Mg II and C IV lines is (3.15–5.19) ×10⁹ M_⊙ (Jiang et al. 2016), which corresponds to a V_BH range of 216–277 km s⁻¹ at R = 0.29 kpc. The black hole contribution to the total potential is significant in the innermost region, but some values of V_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 where σ_star is the velocity dispersion of the surrounding stellar population. When M_BH = 2.97 × 10⁹ M_⊙ (from this dynamical modeling) and σ_star = 250 km s⁻¹ (see Sect. 5.5.2), r_h = 0.21 kpc. When M_BH = 10⁹ M_⊙ (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 σ_star = 100 km s⁻¹ (assuming it is equal to the maximum velocity dispersion measured from our kinematic modeling of the gas), r_h = 0.43 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_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_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_res range of 0.4–1700. Our r_h/r_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.

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 × 109 M⊙ (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.

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

(6)

where M_star, R_e;  star and n_star are the total stellar mass, the effective radius and the Sérsic index. γ and Γ are the incomplete and complete gamma functions. p and b are related to the Sérsic index by (when 0.6 < n_star < 10, and 10⁻² ≤ R/R_e;  star ≤ 10³) and b = 2n_star − 1/3 + 0.009876/n_star (when 0.5 < n_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 degeneracies 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,

(7)

where M_gas and R_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_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 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 (νS_ν = 0.18 ± 0.02 Jy km s⁻¹; Shao et al. 2019), and assuming (Carilli & Walter 2013), we are able to measure the gas mass with an assumed CO-to-gas mass conversion factor (α_CO), namely: . The only free parameter for the gas component is α_CO.

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

(8)

where M_DM and R_DM are the virial mass and radius (M₂₀₀ and R₂₀₀ in Navarro et al. 1996), respectively, c is the concentration parameter, and x = R/R₂₀₀ is the radius in units of the virial radius. The virial mass and radius are correlated with the critical density (ρ_crit): , where and (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_DM.

In summary, for the dynamical modeling, we have six free parameters (M_BH, M_star, n_star, R_e;  star, α_CO, M_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⁸, 10¹⁰] M_⊙ is adopted for M_BH. The quantities of M_star, n_star and R_e;  star are tightly coupled together, and we use prior ranges of [10⁷, 10¹¹] M_⊙, [0.5, 10], and [0, 2.5] kpc, respectively. As for α_CO, we used a range of [0.2, 14] M_⊙/(K km s⁻¹ pc²), which is measured from a sample of nearby AGN, ultra-luminous infrared galaxies, and starburst galaxies (Mashian et al. 2015). For M_DM, we consider a uniform prior of 10¹⁰–10¹³ M_⊙ 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 is consistent with the measurements – (4.17 ± 1.02)×10⁹ and (3.92 ± 0.48)×10⁹ M_⊙ from the Mg II and C IV lines (Jiang et al. 2016), and (1.8 ± 0.5)×10⁹ M_⊙ 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 (α_[C II]) by M_gas = α_[C II]L_[C II]. The derived gas mass corresponds to a α_[C II] of . Given that there are only two degrees of freedom, and the strong degeneracies among some parameters (i.e., M_star, n_star and R_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 × 10⁹ M_⊙ (Feruglio et al. 2018), and present the dynamical results in the middle panel of Fig. 12 and Col. (b) of Table 4. In both situations, the derived stellar mass is on the order of 10⁹ M_⊙ (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 ∼ 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 (Σ_{dust_mass} > 10⁸ M_⊙ kpc⁻²) and the [C II] underlying dust optical depth (τ ∼ 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_ν(T_dust)(1 − exp^−τ_ν)). For both situations, we also consider the CMB as a background (i.e., a black-body radiation of B_ν(T_CMB); T_CMB = 2.73 × (1 + z) K). It should be noted that we only model collisional excitation with molecular hydrogen H₂, 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 . The number density of H₂ is fixed to be 10⁵ cm⁻³. The column density of [C II] for each ring can be calculated with a fixed gas-to-dust mass ratio of ∼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 × 10⁻⁶ 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 Σ_{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 Σ_{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 ∼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.

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 Σ[C II] and Σ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.

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 ∼30° 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 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 ∼ 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 ∼ 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 suggest 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 ∼ 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 different 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) (β 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_⊙. 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_⊙.

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 – α_CO, assuming (Carilli & Walter 2013). We derive an average α_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 α_CO that is 10 times smaller than the overall value, the molecular gas mass is derived to be (2.44 ± 0.28)×10⁸ M_⊙, 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 α_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 and 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. Σ_[C II]/Σ_FIR

Local spiral galaxies have a typical [C II]-FIR luminosity ratio of ∼5 × 10⁻³ when Σ_FIR < 10¹¹ L_⊙ kpc⁻². 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 Σ_[C II] and Σ_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 Σ_FIR than do the comparison samples, but do not probe deeper into the faint end (i.e., Σ_FIR < 10¹⁰ L_⊙ kpc⁻²) 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).

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 ∼ 6 quasars (green triangles connected by dashed lines; Wang et al. 2019) and z ∼ 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 Σ[C II] and ΣFIR as a function of Tdust, respectively.

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 Σ_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_gas > T_[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 Σ_FIR and Σ_[CII] on T_dust. As shown in the right panel of Fig. 14, Σ_FIR is an integral over the gray-body relating to the dust temperature with a power-law index (; black line), larger than that (∼T⁴) of the Stefan-Boltzmann law. However, Σ_[CII] scales as (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. Σ_Mgas and Σ_SFR

The Kennicutt-Schmidt (KS) power-law relation between the gas mass and the SFR surface density, , 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 ∼ 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.00 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 Σ_Mgas 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 (Carilli & Walter 2013) and a conversion factor α_CO ∼ 0.37 M_⊙/(K km s⁻¹ pc²) from the dynamical modeling in Sect. 4.4. In addition, we also measure Σ_Mgas from the [C II] intensity map with α_[C II] ∼ 3.2 M_⊙/L_⊙ from the dynamical modeling in Sect. 4.4. We similarly measure the SFR surface density Σ_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).

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 ΣMgas 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 ΣMgas 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 ΣMgas 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).

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 Σ_Mgas measured from the CO (9–8) line for J2310+1855 is

(9)

and the best-fit KS relation with Σ_Mgas measured from the [C II] line is

(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 ± 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₂ conversion factor α_CO can vary for different galaxy types, position scales, and metallicities. We adopted the α_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 ± 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 α_[C II] might be variable within the galaxy. The depletion times defined by Σ_Mgas/Σ_SFR for the ISM of our quasar host galaxy are 1–100 Myr, increasing from the inner region (corresponding to the highest Σ_Mgas and Σ_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 are slightly enhanced on the red side. The enhancements ([S_{ν, redpeak} − S_{ν, bluepeak}]/S_{ν, bluepeak}; from the double-Gaussian fitting) are about 26%±10% and 18%±14% of the [C II] and CO (9–8) lines, respectively. The two possible peaks are separated by 213 ± 9 and 215 ± 14 km s⁻¹ 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 and km s⁻¹ 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 velocity dispersion are and 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 ^3DBAROLO 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_peak > 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 ^3DBAROLO 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 and 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[[INLINE914]]) 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 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 – – km s⁻¹ for the [C II] line and – 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/σ 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/σ 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/σ 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/σ 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 ^3DBAROLO. They constrain the inclination angle < 30° considering that M_dyn should be larger than M_H2, which is derived from the CO (2–1) emission line (Shao et al. 2019) and adopting an α_co value of 0.8 M_⊙/(K km s⁻¹ pc²). As noted in Sect. 4.4, the values of α_co may differ from source to source (e.g., Mashian et al. 2015), and based on our dynamical model, we derive an α_co value of 0.37 for J2310+1855. The loose constraint on α_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 ^3DBAROLO modeling. The initial guess of the inclination angle of 40° in the ^3DBAROLO 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 calculated with τ_line = −ln(f_trans), where f_trans is the fraction of the continuum emission that is still transmitted, which can be expressed as (S_νcont + S_νline)/S_νcont, where S_νline is the flux density of the continuum subtracted OH⁺ absorption, and S_νcont 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 τ_OH⁺ = 0.13 ± 0.04. The velocity-integrated optical depths τ_OH⁺dv are 35 ± 7 and 44 ± 31 km s⁻¹ by integrating the OH⁺ absorption region and the best-fit Gaussian region of the OH⁺ absorption, respectively. The OH⁺ column density N_OH⁺ is related to the velocity-integrated optical depth τ_OH⁺dv through τ_OH⁺dv = λ³AN_OH⁺/8/π, where A is the Einstein coefficient with a value of 2.11 × 10⁻² s⁻¹ (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_OH⁺ to be (1.7 ± 0.3)×10¹⁴ and (2.2 ± 1.5)×10¹⁴ cm⁻², respectively. Then assuming a relative abundance of log₁₀(N_OH⁺/N_H) = − 7.8 (e.g., Bialy et al. 2019), we measure the neutral hydrogen column density N_H 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 neutral 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 spectra 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 ( versus ), although the sensitivities are comparable (0.17 mJy beam⁻¹ per 17 km s⁻¹ versus 0.23 mJy beam⁻¹ per 8.5 km s⁻¹). To understand 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 , which is similar to the one of Tripodi et al. (2022), and a sensitivity of 0.15 mJy beam⁻¹ per 17 km s⁻¹. The residual spectrum 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] emission 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 matter in the quasar-host system. In Sect. 4.4, we decomposed the J2310+1855 quasar-host system into four dynamical components (black hole, stars, gas, and dark matter) and modeled the potential contribution of each component to the measured circular 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 important to shape the gravitational potential in the central region of the host, contributing ∼120 km s⁻¹ (for M_BH = 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_e;  star 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 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⁹ M_⊙; this is the first time that a black hole mass has been dynamically measured at z ∼ 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 built a strong case for the existence of a > ​10⁹ M_⊙ 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 ∼ 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 ∼ 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⁹ M_⊙) 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 ∼ 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 (∼10⁹ M_⊙) stellar component, especially when we fix the black hole mass to be 1.8 × 10⁹ M_⊙. 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 ∼ 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 ∼ 10¹⁰ M_⊙) 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 α_CO value of .

The dark matter component is less constrained. The mass of the dark matter halo is 16 and 9 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 − σ and the M_BH − M_bulge relations across cosmic time (e.g., Kormendy & Ho 2013), where σ 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:

(11)

where V_gas and σ_gas are the velocity and velocity dispersion of the molecular gas, for which we use the maximum circular velocity (∼245 km s⁻¹; the maximum rotation velocity corrected by the asymmetric drift) and the median velocity dispersion (∼64 km s⁻¹) 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 and 99 km s⁻¹ 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 (σ_stars) are 250 and 233 km s⁻¹ 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 ∼335 ± 20 and ∼300 ± 16 km s⁻¹ for a central black hole of mass 2.97 × 10⁹ and 1.8 × 10⁹ M_⊙ (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 (∼2.6 kpc), the measured ratio of M_BH/M_bary is and 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_BH/M_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 ∼ 6 quasars (e.g., Wang et al. 2019).

Fig. 16. Baryonic mass Mbary (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 MBH and the baryonic mass Mbary. The solid and dashed lines represent the scenarios in which the black hole mass is free and fixed during the dynamical modeling, respectively.

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₁–0₁) 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 () with Σ_Mgas measured from the CO (9–8) line of J2310+1855 is steeper than in the local Universe and z ∼ 2 − 4 starburst samples. However, when Σ_Mgas is measured from the [C II] line, the relation () 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 ^3DBAROLO. 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_rot/σ ratio of our quasar host galaxy are comparable with other high-z populations, for example z ≥ 6 quasar hosts and z ∼ 2 star-forming galaxies. The OH⁺ line shows a P-Cygni profile with an absorption at ∼–400 km s⁻¹, which may indicate an outflow with a neutral gas mass of (6.2 ± 1.2)×10⁸ M_⊙ 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¹⁰ M_⊙. 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 . This is the first time that a SMBH mass has been measured dynamically at z ∼ 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⁹ M_⊙ when the Universe was only ∼0.93 Gyr old. This result is more robust when we fix the black hole mass to be 1.8 × 10⁹ M_⊙ 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.

Acknowledgments

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

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.

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.

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.

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):

(B.1)

where Σ₀ is the central surface brightness, r_s is the scale length, which is the radius at which the surface brightness Σ(r) drops by e⁻¹, 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:

(B.2)

where the dependent variable k is coupled to n, and r_e = kⁿr_s. k satisfies , where Γ and γ are the Gamma function and lower incomplete Gamma function, respectively. k can be approximated as for n > 0.36 (Ciotti & Bertin 1999).

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

(B.3)

(B.4)

(B.5)

(B.6)

where the coefficient A is the height of the peak, (x₀,  y₀) is the center, h_x, h_y are the scale heights along the x and y axis before rotating by θ, and θ 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.

Fig. B.1. Surface plots, projected filled contour plots, and surface brightness distributions for 2D elliptical Sérsic images with (x0,  y0) = (0, 0), hx = 4, and hy = 2. θ and n are π/3 and 0.5 (top left), π/2 and 1 (top right), 2π/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.

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.

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 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 () that of the major axis of the [C II] clean beam size, the rotation angle equal to (=199°) and the ratio of semiminor and semimajor axis – b/a of (), where and come from the [C II] line kinematic modeling (listed in Table 3).

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 is plotted as a white ellipse in each panel.

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:

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 takes 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 Vrot on the galaxy plane. What we can see is the line-of-sight component Vz of the rotation velocity Vrot. The angle between the major axis on the sky plane and the line connecting the galaxy center (x0,  y0) and the data point (X,  Y) is the azimuth angle ζ. The position angle ϕ 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 ζ on the galaxy plane project to (x,  y), R′ and β on the sky plane.

• x₀,  y₀: 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.

• ϕ (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.

• ζ (R): the azimuthal angle measured between the major axis on the sky plane and the line connecting the galaxy center (x₀,  y₀) and the data point (X,  Y) in the galaxy plane. It is related to i (R), ϕ (R), (x₀,  y₀).

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 ζ) is parallel to the xy plane and has no component in the z direction. However, V_Y (=V_rot cos ζ) 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,

(C.1)

(C.2)

(C.3)

Here, we also put the ^3DBAROLO modeling result of CO (9–8) line in Fig. C.2.

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 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.

All Tables

All Figures

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 × 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−1 km s−1. 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−1. 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.

In the text

	Fig. 2. Similar to Fig. 1, but for the CO (9–8) line. The FWHM size of the restoring beam is × . The contour levels are [−3, 3, 6, 9, 12, 15, 18, 21] × rmsco, where rmsco = 0.020 Jy beam−1 km s−1. The spectral resolution is 15.625 MHz, corresponding to 32 km s−1. The rms of the spectrum is 0.31 mJy. The CO (9–8) spectral profile is also asymmetric, with enhancement on the red side.
In the text

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 × is plotted in the bottom-left corner. The contour levels in black are [3, 4, 5] × rmsOH+_em, where rmsOH+_em = 0.012 Jy beam−1 km s−1. The peak flux density is 0.066 ± 0.012 Jy beam−1 km s−1. The overplotted dashed white contours are from the OH+ absorption line, with contour levels of [−3, −2] × rmsOH+_ab, where rmsOH+_ab = 0.008 Jy beam−1 km s−1. 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−1. 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.

In the text

Fig. 4. Spectra of the para-H2O (111–000) line (black histogram) and the scaled para-H2O (202–111) 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−1. The observations may only cover the red tail of the para-H2O (111–000) emission line.

In the text

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 × and × , respectively. The black contours are [−3, 3, 9, 27, 81] × rmscon for both maps. The rmscon at 262 and 147 GHz are 0.02 and 0.01 mJy beam−1, 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.

In the text

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 is plotted in the bottom-left corner of each panel. We measured the residual rms to be 0.037 Jy beam−1 km s−1 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 () that of the [C II] clean beam size, the rotation angle equal to (=199°), and the ratio of semiminor and semimajor axis – b/a of (), where 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 107 L⊙ kpc−2.

In the text

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−1 km s−1 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 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 () that of the CO (9–8) clean beam size, the rotation angle equal to (=198°), and the ratio of semiminor and semimajor axis – b/a of (), where and 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 106 L⊙ kpc−2.

In the text

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(−m2/2) (m = 0.5, 1, 1.5, 2) of the volume. The vertical dashed lines show 16th, 50th (median), and 84th percentiles.

In the text

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.

In the text

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 3DBAROLO, 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 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 3DBAROLO 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 3DBAROLO 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 3DBAROLO 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.

In the text

	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.
In the text

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 × 109 M⊙ (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.

In the text

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 Σ[C II] and Σ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.

In the text

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 ∼ 6 quasars (green triangles connected by dashed lines; Wang et al. 2019) and z ∼ 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 Σ[C II] and ΣFIR as a function of Tdust, respectively.

In the text

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 ΣMgas 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 ΣMgas 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 ΣMgas 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).

In the text

	Fig. 16. Baryonic mass Mbary (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 MBH and the baryonic mass Mbary. The solid and dashed lines represent the scenarios in which the black hole mass is free and fixed during the dynamical modeling, respectively.
In the text

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.

In the text

	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.
In the text

Fig. B.1. Surface plots, projected filled contour plots, and surface brightness distributions for 2D elliptical Sérsic images with (x0,  y0) = (0, 0), hx = 4, and hy = 2. θ and n are π/3 and 0.5 (top left), π/2 and 1 (top right), 2π/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.

In the text

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 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 () that of the major axis of the [C II] clean beam size, the rotation angle equal to (=199°) and the ratio of semiminor and semimajor axis – b/a of (), where and come from the [C II] line kinematic modeling (listed in Table 3).

In the text

	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 is plotted as a white ellipse in each panel.
In the text

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 takes 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 Vrot on the galaxy plane. What we can see is the line-of-sight component Vz of the rotation velocity Vrot. The angle between the major axis on the sky plane and the line connecting the galaxy center (x0,  y0) and the data point (X,  Y) is the azimuth angle ζ. The position angle ϕ 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 ζ on the galaxy plane project to (x,  y), R′ and β on the sky plane.

In the text

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 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.

In the text