A Spatially Resolved [C ii] Survey of 31 z ∼ 7 Massive Galaxies Hosting Luminous Quasars

The [C ii] 158 μm emission line and the underlying far-infrared (FIR) dust continuum are important tracers for studying star formation and kinematic properties of early galaxies. We present a survey of the [C ii] emission lines and FIR continua of 31 luminous quasars at z > 6.5 using the Atacama Large Millimeter Array (ALMA) and the NOrthern Extended Millimeter Array at sub-arcsec resolution. This survey more than doubles the number of quasars with [C ii] and FIR observations at these redshifts and enables statistical studies of quasar host galaxies deep into the epoch of reionization. We detect [C ii] emission in 27 quasar hosts with a luminosity range of L [C II] = (0.3–5.5) × 109 L ⊙ and detect the FIR continuum of 28 quasar hosts with a luminosity range of L FIR = (0.5–13.0) × 1012 L ⊙. Both L [C II] and L FIR are correlated (ρ ≃ 0.4) with the quasar bolometric luminosity, albeit with substantial scatter. The quasar hosts detected by ALMA are clearly resolved with a median diameter of ∼5 kpc. About 40% of the quasar host galaxies show a velocity gradient in [C ii] emission, while the rest show either dispersion-dominated or disturbed kinematics. Basic estimates of the dynamical masses of the rotation-dominated host galaxies yield M dyn = (0.1–7.5) × 1011 M ⊙. Considering our findings alongside those of literature studies, we found that the ratio between M BH and M dyn is about 10 times higher than that of local M BH–M dyn relation on average but with substantial scatter (the ratio difference ranging from ∼0.6 to 60) and large uncertainties.


INTRODUCTION
Corresponding author: Feige Wang feigewang@email.arizona.edu* Strittmatter Fellow The discovery of the strong correlations between the mass of supermassive black holes (SMBHs) and the physical properties (e.g., stellar velocity dispersion, bulge mass, dynamical mass) of their host galaxies (i.e., the M − σ * relation) in the local Universe indicates that SMBHs co-evolve with their host galaxies (see Ko-rmendy & Ho 2013, and references therein).Such relations are thought to arise from the fact that the energy released by accreting SMBHs is able to affect gas kinematic and star formation activity in the host galaxy, thus regulating both the dynamical properties and the stellar mass of the host (e.g.Di Matteo et al. 2005;Hopkins et al. 2008).Therefore, investigating the SMBHgalaxy relations at high redshift is crucial for understanding galaxy formation and evolution at early cosmic times.However, luminous quasars in the rest-frame UV significantly outshine their hosts, prohibiting the detection of the stellar light or the measurement of the M −σ * relation in luminous z ≳ 6 SMBH-galaxy systems even with the Hubble Space Telescope (Mechtley et al. 2012;Marshall et al. 2020).
Thanks to the superb sensitivity of the Atacama Large Millimeter/submillimeter Array (ALMA) and the NOrthern Extended Millimeter Array (NOEMA), the host galaxies of more than 60 z ≳ 6 quasars have been detected in rest-frame far-infrared dust continuum and the fine structure [C II] 158 µm line (e.g., Walter et al. 2009;Wang et al. 2013Wang et al. , 2016;;Willott et al. 2015;Venemans et al. 2016;Decarli et al. 2018;Feruglio et al. 2018;Venemans et al. 2018;Izumi et al. 2019Izumi et al. , 2021;;Wang et al. 2019b;Yang et al. 2019b;Eilers et al. 2020;Pensabene et al. 2020).Several key findings have been established based on these observations.Firstly, copious amounts of dust and gas are present in the majority of these z > 6 quasar host galaxies (see Carilli & Walter 2013;Fan et al. 2023, for reviews).These galaxies, significantly enriched with metals and powered by prodigious star formation with star formation rate (SFR) of ≳ 100 M ⊙ yr −1 , are significantly more massive than typical star-forming galaxies found in deep fields (e.g., Bouwens et al. 2022).Secondly, kinematic modeling of several high spatial resolution (sub-kpc scale) observations of the [C II] line suggests that the most luminous systems seem to have over-massive SMBHs compared with expectations from the local M − σ * relation while fainter systems still follow the local M −σ * relation (e.g., Bañados et al. 2019;Wang et al. 2019b;Venemans et al. 2019;Pensabene et al. 2020;Neeleman et al. 2021;Farina et al. 2022), suggesting that we are witnessing the establishment of the M − σ * relation in the early Universe.Thirdly, the existence of dust-rich companion galaxies adjacent to some fraction of z ≳ 6 quasars indicates that the earliest SMBHs grew in galaxy-rich environment, and that major mergers may be important drivers for rapid SMBH and host galaxy growth (e.g., Decarli et al. 2017;Mazzucchelli et al. 2019;Venemans et al. 2020;Meyer et al. 2022).
The majority of the discoveries mentioned above were based on observations of z ∼ 6 quasars, while observations of higher redshift (e.g., z > 6.5) systems are still limited to only a dozen (e.g., Venemans et al. 2016;Neeleman et al. 2021;Izumi et al. 2021).On the other hand, recent studies show that the spatial density of quasars rises extraordinarily rapidly (by a factor of six) from z ∼ 7 to z ∼ 6, significantly faster than that from z ∼ 6 to z ∼ 3 (Wang et al. 2019a), indicating a rapid buildup of the earliest SMBHs in this short timescale.In the last few years, the number of known z > 6.5 quasars has dramatically increased (e.g., Wang et al. 2017Wang et al. , 2018Wang et al. , 2019aWang et al. , 2021a;;Bañados et al. 2018;Matsuoka et al. 2018Matsuoka et al. , 2019b;;Reed et al. 2019;Yang et al. 2019aYang et al. , 2020aYang et al. ,b, 2021)).This motivated us to perform a comprehensive [C II] survey of quasar host galaxies at z ∼ 7. The main science goals of this survey are to characterize star formation activities in the host galaxies of the earliest SMBHs, to constrain the dynamical masses of quasar host galaxies and thus investigate the M BH − M dynamical relation at z ∼ 7, to perform a census of [C II] emitters in quasar vicinities, and to probe cosmic reionization history and black hole growth by combining the [C II] redshift measurements with optical and infrared spectroscopy.
In this paper, we present an overview of our spatially resolved [C II] survey of 31 luminous z ∼ 7 quasars and discuss the properties of quasar host galaxies based on [C II] and dust emission.In §2, we will describe the quasar sample construction, ALMA and NOEMA observations, and data reduction.In §3, we will present the [C II] line properties, continuum luminosity, SFR, as well as the size measurements of quasar host galaxies.In §4, we will highlight the diverse kinematic properties of quasar host galaxies, present first-order constraints on their dynamical masses, and discuss the M − σ * relation of these systems.Finally, we will conclude in §5.Throughout the paper, we adapt a flat cosmology model with H 0 = 70.0km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.

OBSERVATIONS AND DATA REDUCTION
The parent sample of our [C II] survey was built to include all 56 z > 6.5 quasars with M 1450 < −25.0 known at the time of designing this program.The redshift and M 1450 distributions of these quasars and all other known quasars at z > 6.5 are shown in Figure 1.In the literature, 13 quasars had either ALMA or NOEMA observations at the time of designing this program and thus were not targeted by this survey.An additional two quasars, J0244−5008 and J0020−3653 (Reed et al. 2019), were proposed to ALMA in Cy- a There is no infrared spectroscopic observation for this object; the redshift was measured from the Lyα line by Matsuoka et al. (2018).
b References for the quasar discovery (first number) and for zMgII and M1450 measurements (second number).References: 1: cle 5 by the quasar discovery team and were also excluded from this survey.We observed the remaining 31 z > 6.5 quasars with M 1450 < −25.0 with ALMA and NOEMA: 26 quasars with Decl.< 30 • were targeted with ALMA and the other 5 quasars with Decl.> 30 • were observed with NOEMA.One quasar in our sample, J0252-0503, was observed by both ALMA and NOEMA.
The NOEMA observation of the gravitationally lensed quasar J0439+1634 at z = 6.52 was published in Yang et al. (2019b) and the high-resolution (C43-5) ALMA observation of this object was published in Yue et al. (2021).In this paper, we present the ALMA C43-3 observation of J0439+1634 obtained from this program.
Although the ALMA observations of J1007+2115 and J0313-1806 were reported in Yang et al. (2020a) and Wang et al. (2021a), respectively, we include the observations of these two quasars in this paper for completeness.The basic information of these 31 quasars as well as the observation details are listed in Table 1.In the following sections, we only reports the measurements for objects targeted by this survey but includes objects from literature for statistically analyses.

ALMA observation and data reduction
Our ALMA observations span two cycles with one Cycle 6 program (program ID: 2018.1.01188.S, PI:  These observations were designed to reach a flux sensitivity of ∼0.3 mJy/beam per 100 km s −1 and a resolution of 0. ′′ 4-0.′′ 7, in order to detect the [C II] emission of nearly all quasars and marginally resolve the most extended quasar host galaxies.We tuned two spectral windows (SPW, 0 and 1) centered at the expected frequency of [C II] and the other two SPWs (2 and 3) centered at about 15 GHz away from the expected [C II] line for observing the continuum emission.Observations were carried out between 2018 December and 2020 March with 42-50 12m antennas.The on-source exposure time for each target was designed to be ∼ 15 minutes and targets were observed with either the C43-3 or C43-4 configuration, with the final on-source time depends on the observing conditions.The detailed ALMA observational information is listed in Table 1.
All the ALMA data were processed using the CASA (McMullin et al. 2007) pipeline for ALMA using the default calibration procedure.The Cycle 6 data were calibrated with CASA version 5.4.0 and the Cycle 7 data were calibrated with CASA version 5.6.1.All the data were then analyzed using CASA version 5.6.1.We imaged the data cubes using Briggs cleaning via the CASA task tclean with robustness parameter r = 2.0, corresponding to natural weighting, to maximize the signal-to-noise ratio (S/N) of our observations.Since we will compare our results with those of lower redshift (z ∼ 6) quasars from Decarli et al. (2018), the imaging process was designed to closely follow that used by Decarli et al. (2018).The imaging process includes the following steps: (1) We generate a dirty continuum + line map for each target by collapsing the data in SPWs 0&1 with tclean.This map is used for identifying the positions of detected sources, creating box masks for detected sources and estimating the sensitivity of the observations.
(2) We then create data cubes for SPWs 0&1 and 2&3 with tclean.During this step, we use the object mask regions generated from last step and re-bin the data cube to 23.4 MHz (∼ 30 km s −1 at the quasar redshifts).
(3) We extract a 1D spectrum of each object with an aperture diameter of 1. ′′ 5. Note that Decarli et al. (2018) extracted 1D spectra on a single-pixel basis because their data have a lower resolution.
(4) The spectra extracted from the SPWs 0&1 data cubes are then fitted with a flat continuum and a Gaussian function.The continuum flux, line flux, line peak and line width are listed in Table 2.
(5) We then subtract the continuum for the data cube from SPWs 0&1 using the uvcontsub task in the UV plane.
(6) The [C II] line map is collapsed from the data cube in the frequency range set by the line peak ±1.4σ line from the 1D spectral fitting generated in Step 4. For a Gaussian line profile, this integration corresponds to recovering ∼83% of the total line flux.For the object (J1129+1846) without a successful Gaussian fitting from the 1D spectra, we assume the line has the same redshift as that derived from the Mg II line and a width of FWHM = 300 km s −1 .
(7) The continuum maps are generated by collapsing the entire line-free channels in SPWs 0&1&2&3.
The [C II] line maps and the continuum contours are shown in Figure 2 and the extracted 1D spectra are shown in Figure 3.

NOEMA observation and data reduction
Our NOEMA observations were obtained in the W18 cycle (program ID: W18EJ, PIs: B. Venemans and J. Yang).All quasars were observed using 10 antennas with the C array configuration during winter 2018.We used the wideband correlator PolyFiX, which delivers two ∼7.7 GHz wide sidebands.The expected frequency of the [C II] line of each quasar was tuned at the center of the inner baseband of the upper sideband.The total telescope time was designed to be ∼ 4 hours (onsource ∼ 2.5 hours) per target in order to reach ∼0.7     mJy/beam per 100 km s −1 .The detailed NOEMA observational information is also listed in Table 1.
The NOEMA data were analyzed using the mapping software from the GILDAS 1 suite.We use natural weighting and resampled the spectral axis in 50 km s −1 wide channels.We extract the [C II] 1D spectra from the peak positions for those [C II] detected targets and extract the spectra from the optical positions for the [C II] undetected targets.The extracted 1D spectra are shown in Figure 3. Similar to the ALMA data reduction, the [C II] line intensity map is collapsed from the data cube in the frequency range of ±1.4σ line from the line peak frequency (using values derived from Mg II redshifts for [C II] undetected sources).The continuum maps are generated by collapsing all the line-free channels.The [C II] line maps and the continuum contours of the NOEMA targets are also shown in Figure 2.

[C II] and Infrared Luminosity Measurements
As described in §2, the extracted spectra are shown in Figure 3.We fit the spectra by assuming a flat dust continuum emission and a Gaussian profile for the [C II] emission line.The derived line peak frequency, full width at half maximum (FWHM), and the line flux as well as their associated 1σ uncertainties are listed in 1 http://www.iram.fr/IRAMFR/GILDASTable 2.We consider a line detection to be significant if its integrated line flux is > 3σ and marginal if its integrated line flux is between 2σ and 3σ.Out of 31 targeted objects, only four quasars (J0829+4117, J0839+3900, J1129+1846, and J1216+4519) are undetected in [C II] and two quasars have marginal line detections at a 2.6σ level (J0213−0626) and a 2.1σ level (J0837+4929), respectively.Three (J0829+4117, J0839+3900, and J1216+4519) out of the four [C II]-undetected quasars are from the NOEMA observations, which have a depth that is two times shallower than that of the ALMA observations.The FIR continuum flux and flux uncertainty measured from the 1D spectral fitting are also listed in Table 2.All objects observed with ALMA, except for J1129+1846, are detected in the dust continuum at high significance.Similar to the [C II] results, the FIR continua of J0829+4117, J0839+3900, and J1216+4519 are not detected by NOEMA.J1129+1846, the only known radio-loud quasar in our sample, is intrinsically faint in both [C II] and FIR (see also Khusanova et al. 2022).
We also derived the line flux and the size of the [C II] emitting region for the ALMA-detected objects by fitting the [C II] line maps shown in Figure 2. The fit is performed with the imfit task within CASA by assuming a 2D Gaussian profile with the centroid, integrated flux, major and minor axes, and the position angle as free parameters.We successfully fit all the detected [C II] emission except for J0213−0626 because of the low significance of its detection.Since the [C II] maps are created by integrating over ±1.4σ line , which encloses ∼ 83% of the total line flux for a Gaussian profile, we take this factor into account when computing the total line flux (see also Decarli et al. 2018).The line flux and the uncertainty measured from the 2D fitting are also listed in Table 2.For J0213-0626, we use the line flux calculated from the intensity map with a 1. ′′ 5 diameter aperture after correcting for the missing flux.We note that the line fluxes measured from the 2D map could be different from those measured in 1D spectra in some cases, especially when the galaxy is very extended (e.g., J0411-0907).This is mainly because the 1D spectral extraction could not enclose all flux from the quasar host galaxies.Therefore, we use the 2D line fluxes in the following analyses (see also Decarli et al. 2018).For the NOEMA-detected objects, we calculate the [C II] flux from the intensity map with a 3. ′′ 0 diameter aperture after correcting for the missing flux.The 3σ upper limits are reported in Table 2 for the [C II]-undetected targets.We then measure the line luminosity following   a Flux was estimated from a circle centered at the quasar position with a diameter of 1. ′′ 5.
b 3σ upper limit.
c Flux was estimated from the taper image with a beam size of 1. ′′ 0.
d Flux was estimated from a circle centered at the quasar position with a diameter of 3. ′′ 0.
where  a Dynamical mass measured from Equation 7.
b Dynamical mass measured from Equation 8.
significantly increases the number of [C II] detections in quasar host galaxies at z > 5.7 and more than doubles the number of [C II] detections in quasar host galaxies at z > 6.5.This program provides a complete coverage in [C II] of the most luminous quasars (i.e., M 1450 < −25.0) known at z > 6.5.In the left panel of Figure 5, we show the scatter plot between L [CII] and L bolometric for both our sample and the literature sample.In order to examine whether these two qualities correlate with each other, we did a Spearman test on measurements from both this work and literature.The Spearman test gives a correlation coefficient of ρ = 0.43 and a chance probability of p = 1.7 × 10 −4 , indicating a moderate correlation between these two quantities though a large scatter is present.This is consistent with what has been found with a smaller quasar sample at z ∼ 6 (e.g.Decarli et al. 2018;Izumi et al. 2019).This suggests that more luminous quasars in general hosted by galaxies with brighter [C II] emission, however, the correlation could be biased since the [C II] observation were only performed for UV luminous quasars (i.e., UV faint quasars with bright [C II] emission is missed by current surveys).
Since most of the ALMA-observed targets are marginally resolved (i.e., emission extent is ≳ 1.5× the beam size), we fit their continuum maps with imfit to derive the continuum flux and size.The measured dust continuum fluxes are listed in Table 2. J1129+1846 is only marginally detected in the 2D continuum map and we estimate the continuum flux from a taper image (with a 1. ′′ 0 beam) without measuring the size.The continua of the quasar host galaxies observed by NOEMA are not spatially resolved, therefore we use the peak values from the 2D continuum maps as their continuum fluxes.The FIR rest-frame 42.5-122.5µm luminosities, L FIR , and the total infrared (TIR, rest-frame 8-1000 µm) luminosities, L TIR , are then estimated from the 2D continuum flux measurements.In order to infer the FIR and TIR luminosities, we used two different models: a galaxy template from Haro 11 which has a relatively low metallicity and is thought to be a good analogue of high-redshift galaxies (e.g.Lyu et al. 2016), and a modified blackbody in the optically thin regime, i.e., assuming the dust optical depth τ dust ≪ 1 (e.g.Beelen et al. 2006;Venemans et al. 2018).For the case of the Haro 11 template, we simply scale the template to the measured continuum flux of each target and then derive L FIR and L TIR by integrating the galaxy template from 42.5 to 122.5 µm and from 8 to 1000 µm, respectively.For the case of the graybody model, we assume a dust temperature of T dust = 47 K and an emissivity index of β = 1.6 (Beelen et al. 2006).Since the cosmic microwave background reduces the flux densities that we can measure from sources at high redshift, we also take the CMB effect into account, following da Cunha et al. ( 2013) when fitting the graybody.The measured infrared luminosities for both cases are listed in Table 3.The L TIR measured using the galaxy template from Haro 11 is in general 1.5 times higher than that was measured from the graybody fitting.Since the modified graybody has been more broadly used in the literature, we only use the luminosities derived from the graybody fitting in the following analyses.There is no precise magnification measurement in FIR of the gravitationally lensed quasar, J0439+1634, thus we only report the flux measurements of this source and discard the lumi-nosity and size measurements in the following analyses.We refer the readers to Yang et al. (2019b) and Yue et al. ( 2021) for a detailed characterization of the host galaxy of this quasar.
The L FIR (measured from the graybody fitting) of our targets has a broad range between ∼ 5 × 10 11 L ⊙ and ∼ 1.3 × 10 13 L ⊙ , similar to what has been found for slightly lower redshift quasar host galaxies (e.g.Venemans et al. 2020) but significantly brighter than z ≳ 6 optically-selected star-forming galaxies (e.g.Smit et al. 2018;Bouwens et al. 2022).In the right panel of Figure 5, we show the scatter plot between L FIR and L bolometric .The Spearman test gives a correlation coefficient of ρ = 0.48 and a chance probability of p = 7.4 × 10 −6 .This indicates a mild positive correlation exists between L FIR and L bolometric although with significant scatter.This is also consistent with similar studies at slightly lower redshifts by Izumi et al. (2018).We note that Venemans et al. ( 2020) claimed a weaker or even no relationships between L FIR and L bolometric , but that could be caused by the lack of low-luminosity objects in their sample.The L bolometric is a proxy for ongoing black hole growth while the L FIR and L [CII] are good tracers of star formation in the quasar host galaxies.The mild positive L FIR − L bolometric and L [CII] − L bolometric correlations of these quasars indicate that the more active SMBHs in the early Universe generally resides in galaxies with more active star forming activities.However, the large scatters in both relations could be caused by a variety of effects, including different timescales for black hole accretion and star formation.However, we noted that these correlations could be biased since current observations were targeted for UV luminous quasars and a population of UV faint quasars with strong activate star formation activities could be missed by our survey.

[C II] and FIR Continuum Size Measurements
As described above, we fit the 2D maps of both [C II] and continuum emission with imfit for the sources detected with ALMA.The derived sizes are listed in Table 4.The ratio between the [C II] major axis size and the beam major axis size ranges from ∼1.3 to ∼3.0 (see the column S/B [CII] of Table 4), suggesting that the [C II] emission regions of all quasars detected by ALMA are marginally resolved.The deconvolved sizes of the [C II] range from 0. ′′ 44×0.′′ 34 (or 2.20 kpc × 1.70 kpc) to 1. ′′ 85×1.′′ 41 (or 9.82 kpc × 7.48 kpc) with a median major axis FWHM of 0. ′′ 7 or 3.5 kpc.The ratio between the continuum emitting region size and the beam size ranges from ∼1.1 to ∼2.4 (see the column S/B cont of Table 4).The deconvolved major axis size of the continuum emitting region ranging from 0. ′′ 20×0.′′ 16 (or 1.06 kpc × 0.85 kpc) to 1. ′′ 34×0.′′ 79 (or 7.23 kpc × 4.26 kpc) with a median major axis FWHM of 0. ′′ 4 or 2.1 kpc.In Figure 6, we show the comparison of the [C II] effective sizes and continuum effective sizes, where the effective sizes were derived with FWHM maj FWHM min .In this plot, we also include all effective size measurements collected from the literature.In general, the continuum emission is more concentrated than the [C II] emission, consistent with previous studies of both quasar host galaxies (e.g.Venemans et al. 2020;Novak et al. 2020;Walter et al. 2022) and high-redshift star-forming galaxies (e.g.Capak et al. 2015;Fujimoto et al. 2020).This could be caused by the fact that the radiation in the central region closer to the quasar is more intense which yields a lower L [CII] /L FIR (i.e., the so-called [C II] deficit).In §3.5, we will discuss the [C II] deficit in more detail.
Studies of the dynamical mass of quasar host galaxies with low-resolution [C II] observations usually assume an inclination angle (e.g.i ∼ 55 • ; Willott et al. 2015;Decarli et al. 2018) and a disk size for the hosts (e.g.D ∼ 4.5 kpc; Wang et al. 2016).These values are derived from the early ALMA observations of five quasar host galaxies at z ∼ 6 (Wang et al. 2013).Using the spatially resolved observations from ALMA, we can update the median size and inclination angle of quasar host galaxies based on a much larger quasar sample (23 quasars from this work and 49 quasars from the literature).Following previous work (e.g., Wang et al. 2013;Willott et al. 2015), we estimate the galaxy diameter D to be 1.5× the FWHM major axis from the [C II] intensity map (i.e., full width at 20% of the peak intensity for a Gaussian profile), and we estimate the inclination angle i from the ratio of minor (a min ) and major (a maj ) axes according to i = cos −1 (a min /a maj ).The median values are D = 5.0 kpc and i = 46 • , slightly different from but still consistent with the typical values assumed in previous works (e.g., Wang et al. 2013;Willott et al. 2015).Note that when creating the moment zero maps of the [C II] emission (as described in §2), we only collapsed the data cube within the frequency range of ±1.4σ line .To assess whether this choice might miss any extended gas emission with higher velocities in the quasar host galaxies, we also measured the effective sizes of the [C II] emitting regions by collapsing the data cube within the frequency range of ±3.0σ line .From this experiment, we discovered that the ratio of the effective size derived from the wider cube width to that of the size derived from ±1.4σ line falls within the range of 0.86 to 1.27, with a median of 1.04 ± 0.10.This indicates the robustness of our size measurement.However, we acknowledge the possibility that we may have missed fainter emission at the outskirts of the galaxies due to limited data quality.

Star Formation Rates
One of the main goals of this program is to measure the SFR of the targeted quasar host galaxies, which is important to our understanding of how the growth of the central SMBHs is related to the assembly of their host galaxies.In this work, we infer the SFR of quasar host galaxies using both [C II] and dust continuum emissions.
The [C II] line is one of the primary fine-structure cooling lines and has been suggested as a good SFR tracer (e.g., De Looze et al. 2014;Herrera-Camus et al. 2015).We use the empirical relation calibrated by De   2014) for high-redshift galaxies to estimate the SFR of these [C II] detected quasar host galaxies: (2) For the dust continuum-based SFR, we use two different templates as discussed above when measuring the L FIR : a modified graybody and the Haro11 galaxy template.In the case of the graybody in the optically regime, we assume a dust temperature of T dust = 47 K and an emissivity index of β = 1.6 (Beelen et al. 2006) and use the relation calibrated by Murphy et al. (2011): (3) We use the star formation law calibrated by Lyu et al. (2016) for the Haro 11 template: (4) The SFRs derived from all three estimators are listed in Table 3.In Figure 7, we show the relation between the SFR derived from the [C II] lines and the dust continuum-based SFR inferred from the graybody model.The SFRs derived from these two methods are broadly consistent with each other, with a range of SFR ∼ 100 − 3000 M ⊙ yr −1 and both are about two times lower than the SFR derived from the Haro 11 template.However, the [C II] based SFR at the highluminosity end is systematically lower than that derived from the FIR-based SFR.This is likely due to the fact that the empirical relation from De Looze et al. (2014) was calibrated to star-forming galaxies, while the FIRluminous galaxies have a significant [C II] deficit effect at the high-luminosity end.We will discuss this in §3.5 in more detail.Since the SFR graybody is more broadly used in the literature and more consistent with the SFR [CII] , we will only use the SFR graybody and SFR [CII] in the following analyses.We note that the SFRs derived from all three methods have large uncertainties given that we have no knowledge about the metallicity and dust temperature of the quasar host galaxies.In particular, the scaling relations were derived from nominal star-forming galaxies (L [CII] and L TIR,graybody ) or local galaxies (Haro 11) which have not been validated as accurate templates for high-redshift quasar host galaxies.Such uncertainties can be better understood with future ALMA band 8/9/10 and JCMT/SCUBA-2 high frequency observations and JWST mid-infrared observations by capturing the peak of the dust continuum emission and improving our understanding of the AGN contribution to dust heating, respectively.The inferred SFR surface densities (i.e., the graybody-based SFR divided by our estimates of the continuum emitting size) are in the range of ∼1-800 M ⊙ yr −1 kpc −1 , below the star formation Eddington limit (∼ 1000 M ⊙ yr −1 kpc −1 ; Walter et al. 2009Walter et al. , 2022)).

Dust and Gas Mass
The dust masses of quasar host galaxies can be estimated from the observed FIR flux density with an assumption of the dust temperature, emissivity index with an optically thin case at the observed wavelength using the following equation: where F ν is the continuum flux density measured from the 2D continuum map as listed in Table 2, D L is the luminosity distance, κ ν (β) is the dust mass opacity coefficient, and B ν is the Planck function.Following Venemans et al. ( 2018), the opacity coefficient is set by κ ν (β) = 0.77 (ν/352 GHz) β cm 2 g −1 .To be consistent with our measurements of L FIR and the reports from most of the literature, we assume T = 47 K and β = 1.6 in Eq. 5, and take the CMB effect, B ν (ν, T CMB ), into account when estimating the dust mass.The estimated dust masses of these continuum-detected quasar host galaxies range from ∼ 3 × 10 7 M ⊙ to ∼ 8 × 10 8 M ⊙ and are listed in Table 3.The dust contents of these quasar host galaxies are similar to those found in z ∼ 6 quasar host galaxies (e.g.Izumi et al. 2018;Venemans et al. 2018).By assuming a dust-to-gas ratio of 100 and a molecular-to-total gas mass fraction of 0.75 (e.g.Neeleman et al. 2021), we can roughly estimate the gas mass  in these galaxies to be in the range of 4 × 10 9 M ⊙ to 1 × 10 11 M ⊙ , suggesting that these quasar host galaxies are among the most massive and gas-rich galaxies in the early universe.However, we note that the dust mass inferred from a single continuum measurement has large uncertainties and could increase by a factor of two if we assume a lower dust temperature (i.e.T = 35 K).Future high-frequency FIR continuum observations and CO line observations are needed to improve the constraints on the dust and gas masses.
The mass of singly ionized carbon can be derived from the strength of the [C II] emission line.Following Venemans et al. (2017b), we calculate this quantity using the following equation by assuming that the emission is optically thin Weiß et al. (2005).
where T ex is the excitation temperature, and Q(T ex ) = 2+4e −91.2/Tex is the [C II] partition function.By assuming T ex = 100 K (e.g.Meijerink et al. 2007), we derive the mass of singly ionized carbon from the measured [C II] luminosity, L ′ [CII] (in units of K km s −1 pc 2 ).The measured M C + has a broad range between 8 × 10 5 M ⊙ and 1.6 × 10 7 M ⊙ .The M C + values for [C II]-detected galaxies are also given in Table 3.

[C II] Deficit
The ratio between L [CII] and L FIR is an indicator of the contribution of the [C II] line to the cooling of the interstellar medium.The L [CII] /L FIR was found to be ∼ 10 −3 for Milky Way-like star-forming galaxies, while it is about an order of magnitude smaller in FIR-luminous systems like ULIRGs (see Hodge & da Cunha 2020, for a recent review).Although the underlying physical causes of the observed '[C II] deficit' are still debated, subsequent studies have expanded the investigation of the [C II] deficit to high-redshift galaxies and to sub-kpc scales, indicating that the physical scale of [C II] deficits must be sub-kpc (e.g.Lamarche et al. 2018) and could be caused by a locally intense radiation field (e.g.Rybak et al. 2019).
In Figure 8, we show L [CII] /L FIR as a function of both L FIR and the FIR luminosity surface density, Σ FIR .We compute Σ FIR = L FIR /(2πR 2 eff,cont ) based on the 2D Gaussian fit of the continuum map as detailed in §3.2, where R eff,cont = R maj,cont R min,cont .Following Decarli et al. (2018), we include a factor 2 in the denominator when computing Σ FIR to account for the fact that R eff,cont roughly encompasses half of the total light.The L [CII] /L FIR ratios of these high-redshift quasar host galaxies span nearly two orders of magnitude.In the right panel of Figure 8, we also show the [C II] deficit relations found from local infrared-luminous galaxies.The relation between L [CII] /L FIR and Σ FIR in quasar host galaxies is similar to that of local infraredluminous galaxies, indicating that the [C II] deficit is not strongly affected by the presence of luminous AGN (see also Venemans et al. 2020).However, we need to acknowledge that the current dataset does not allow us to understand the physical origins of the [C II] deficit.In this work, we assumed T = 47 K for all quasar host galaxies which could introduce significant uncertainties on the L FIR measurements.In addition, the current observations only marginally resolve the quasar host galaxies and do not allow us to study Σ FIR as a function of location.Furthermore, we do not have metallicity constraints on these quasar host galaxies.Future multiplefrequency continuum observations, high-resolution observations of both dust continuum and [C II] and multiline diagnostics are necessary to constrain the physical origins of the [C II] deficit.
As discussed in §3.North is up and east to the left.In the moment 0 map, the solid black contours mark the +3, 4, 5, 6 isophotes for those objects with continuum peak at < 7σ, or the 50%, 60%, 70%, 80%, 90% of the continuum peak emission for those objects with ≥ 7σ peak detections.The synthesized beam of the observations is shown in the bottom-left corner of each panel.The solid contours and dashed contours in the moment 1 maps mark the positive velocities and negative velocities, respectively.The contour levels for both positive velocities and negative velocities are in steps of 20 km s −1 .tio of L [CII] /SFR is different for galaxies with different Σ SFR , consistent with the [C II] deficit observed in these galaxies: log 10 (L [CII] /SFR) = 7.03 ± 0.21 for galaxies with Σ SFR in the range of 10 8 − 10 11.2 L ⊙ kpc −2 and log 10 (L [CII] /SFR) = 6.53 ± 0.30 for galaxies with Σ SFR in the range of 10 11.2 − 10 12 L ⊙ kpc −2 .The comparison of these measurements with the L FIR based SFR is shown in Figure 9.This figure indicates that the [C II]based SFR after considering the Σ SFR is more consistent with the L FIR -based SFR when compared with Figure 7.This test suggests that the [C II] deficit in quasar host galaxies is broadly consistent with star-forming galaxies at lower redshifts (see also Figure 8) and it is necessary to consider the FIR radiation field intensities when estimating the SFR of quasar host galaxies using the [C II] line.
4. KINEMATICS AND DYNAMICAL MASS 4.1.Kinematics and Morphology Neeleman et al. (2021) presented kinematic analyses of 27 quasar host galaxies at z ∼ 6 with ALMA at ∼0. ′′ 1-0.′′ 25 resolution.They found that about one-third of the z ∼ 6 quasar host galaxies (10 out of 27) show a smooth velocity gradient, one-third (9 out of 27) of the quasar host galaxies have disturbed [C II] emission profiles, and the remaining one-third (8 out of 27) show a velocity profile without a clear velocity gradient, which is consistent with the emission arising from a dispersiondominated system.They also found that all galaxies have intrinsically large velocity dispersions despite different kinematics.
To investigate the kinematics of the quasar host galaxies from our z ∼ 7 quasar sample, we re-image the ALMA [C II] data cube using the tclean task with robust parameter r = 0.5 to obtain images with slightly higher spatial resolution (∼0.′′ 2-0.′′ 5).Three objects We then create a mask image by selecting all pixels with S/N> 3 in the intensity map.We finally create a velocity map (moment 1) and a velocity dispersion map (moment 2) of each target using the mask image created from previous step and selecting frequencies between −2.0σ line and +2.0σ line with the immoments task within CASA package.The moment maps are shown in Figure 10 and Figure 11.Similar to what has been found by Neeleman et al. (2021), all galaxies show large velocity dispersions as indicated by the moment 2 maps, although the marginally resolved maps are affected by beam smearing and by the uncertainty of the inclination angles.Förster Schreiber et al. (2009) proposed a general method for distinguishing rotationdominated galaxies from dispersion-dominated systems by using the ratio of the projected velocity gradient (v proj , observed difference between the maximum and minimum relative velocities without correcting for inclination) and the integrated line width (σ int , measured Figure 12.Early universe SMBH-host galaxy co-evolution.In all panels, the orange circles and orange squares denote our measurements of z > 6.5 quasars with rotation-dominated kinematics and dispersion-dominated kinematics or disturbed galaxies, respectively.The open blue points represent data of z ∼ 6 quasars collected from literature, and the gray points are the measurements of local galaxies from Kormendy & Ho (2013).The best fits for the local MBH − σ and MBH − M bulge relations derived by Kormendy & Ho (2013) are shown as a red dashed line with the shaded regions denoting the 1σ scatter of MBH and σ/M bulge .The error bars of MBH in this plot include the uncertainty from the Mg II black hole mass scaling relation, which is ∼ 0.5 dex (Vestergaard & Osmer 2009).Left: Black hole mass versus σ derived from FWHM [CII] .Middle: The black hole mass versus σ derived using the method proposed by Ho (2007) which corrects for the effects of inclination and gas turbulence.The uncertainties of σHo includes both measurement uncertainty and the relation uncertainty in Ho (2007).Right: The black hole mass versus dynamical mass estimated from Equation 7.
from 1D line profile) for marginally resolved observations.Based on simulations of disk galaxies, they found that galaxies with v proj /2σ int > 0.4 can be classified as rotation-dominated systems.Such method has been successfully used in ALMA observations of high-redshift galaxies (Smit et al. 2018).We measure v proj from the velocity map shown in Figure 10 and Figure 11 and σ int from the measured line widths in Table 2.The measured v proj /2σ int for the 23 quasar host galaxies are listed in Table 2. Nine of the 23 galaxies have v proj /2σ int > 0.4, however, J0923+0753 shows obvious distorted morphology.Further more, two galaxies have v proj /2σ int < 0.4 but with clear velocity gradient (J0246-5219 and J0910-0414).We therefore conclude to classify J0038-1527, J0224-4711, J0229-0808, J0246-5219, J0439+1634, J0525-2406, J0706+2921, J0910-0414, J1058+2930, and J2002-3013 (43%) as rotationdominated galaxies (Fig. 10).The remaining galaxies (57%) are classified are likely either dispersiondominated galaxies or disturbed by galaxy mergers.Such fraction is similar to that found in the z ∼ 6 quasar host galaxies in Neeleman et al. (2021) if we also separate their galaxies into dispersion-dominated kinematics (13 out of 27) and galaxies showing clear velocity gradients (14 out of 27) solely based on velocity maps.We note that such method could cause some misclassification limited by the low-resolution of the data used in this work and future high-resolution observations are needed to improve our classifications.

Dynamical Mass
In the following, we estimate the dynamical mass for those rotation-dominated galaxies (except for the lensed galaxy J0439+1634) by assuming that the [C II] emission is originating from a thin rotating disk.As mentioned in §3.2, the disk size, D, is given by 1.5× the FWHM major axis from the [C II] intensity map, and the factor 1.5 is used to account for spatially extended low-level emission (e.g.Wang et al. 2016).The inclination angle, i, is given by the ratio of minor (a min ) and major (a maj ) axes, i = cos −1 (a min /a maj ).The circular velocity is expressed as v circ = 0.75 × FWHM [CII] /sin i.The dynamical mass within D can then be estimated by: We use the size and the inclination angle measured in §3.2.The dynamical masses of these galaxies estimated using Equation 7 are listed in Table 3. Neeleman et al. (2021) note that the dynamical mass measured from Equation 7 using marginally resolved observations could be overestimated because of the presence of high velocity dispersion.They suggested an empirical relation based on their high-resolution [C II] observations: The dynamical masses measured using Equation 8, which are ∼ 2.3 times lower than those estimated from Equation 7, are also listed in Table 3.The dynamical mass spans a broad range of ∼ (0.1−7.5)×10 11 M ⊙ .We emphasize that the M dyn estimated from both equations have large uncertainties, we will improve the M dyn measurements based on kinematic modeling of future highresolution ALMA observations (Yang et al. in prep).

SMBH-host Co-evolution
The recent measurement of the spatial density of the brightest quasars suggests an extraordinary rapid rise of the bright quasar population from z ∼ 6.7 to z ∼ 6, indicating a rapid buildup of SMBHs in a short time (∼120 Myr) from z > 6.5 to z ∼ 6 (Wang et al. 2019a).Therefore, investigating the co-evolution between these earliest SMBHs and their host galaxies is crucial for understanding the assembly of early massive galaxies.Such investigation has been extensively explored in the past few years at z ∼ 6 (e.g.Wang et al. 2013;Willott et al. 2015;Izumi et al. 2019;Neeleman et al. 2021), but limited to a small sample of quasars at z > 6.5 (e.g.Venemans et al. 2017a;Neeleman et al. 2021).In this section, we investigate the early co-evolution of SMBHs and their host galaxies with both the M BH − σ and M BH − M dyn relations by combining the sample of z > 6.5 quasars studied in this work with data collected from the literature for z ∼ 6 quasars.
The BH masses of all z > 6.5 quasars studied in this work, except for J0213-0626, J0229-0808, and J0430-1445, are measured from the Mg II emission line and are collected from the literature (Bañados et al. 2021;Wang et al. 2021b;Yang et al. 2021).We also collected the BH masses for other [C II]-detected quasars where available (e.g.Onoue et al. 2019;Shen et al. 2019;Schindler et al. 2020;Wang et al. 2021b;Farina et al. 2022).To be consistent, all M BH values are calculated using the single-epoch black hole mass estimator (Vestergaard & Osmer 2009) and under the same cosmological model used in this work.
Since detecting the stellar light of these quasar host galaxies is challenging (e.g.Marshall et al. 2020), σ and/or M dyn are usually estimated from the [C II] observations.First, we simply assume σ = FWHM/(2 2 ln(2)) following most of the literature (e.g.Venemans et al. 2017a).Secondly, Ho (2007) suggests that one can relate molecular or atomic gas in a disk to stellar bulges using the line width measured at 20% of the peak intensity after correcting for the inclination and other effects.Following Willott et al. (2015), we set the [C II] line full-width at 20% equal to 1.5 times the FWHM since the lines are approximately Gaussian.We then estimate σ following Ho (2007) by correcting the turbulent broadening and the inclination factor (assuming a median inclination angle of 45 • ).
In Figure 12, we show the M BH − σ and M BH − M dyn relations for all z > 5.7 quasars with both [C II] detections and Mg II-based black hole mass measurements as well as the measurements and derived relations from local galaxies (Kormendy & Ho 2013).Overall the newly available measurements of z > 6.5 quasars from this work occupy the same regions as those of z ∼ 6 quasars collected from the literature in all the plots.In the left panel of Figure 12 we show the σ derived from FWHM/(2 2 ln(2)) versus M BH .Similar to previous findings (e.g.Wang et al. 2010Wang et al. , 2016;;Venemans et al. 2017a;Wang et al. 2019b), most of the SMBHs are significantly over-massive compared to the local M BH − σ relation.On the other hand, the deviation is smaller if we use the σ derived following Ho (2007) as shown in the middle panel of Figure 12.The reason is that the method proposed by Ho (2007) corrects for the inclination effect and thus has a larger σ than the simple estimate from FWHM [CII] (see also Wang et al. 2010;Willott et al. 2015).Nevertheless, we caution that both methods have large uncertainties given the high turbulent velocities and unclear kinematics of the [C II] emitting gas, and the large uncertainty on the inclination angle estimates.In the right panel of Figure 12, we show the M BH − M dyn relation, where M dyn of rotation-dominated quasar host galaxies are estimated using Equation 7. The majority of z ≳ 6 SMBHs are slightly over-massive compared with expectations from the local M BH − M bulge relation, although some of them still follow the local trend.
Again as noted by Neeleman et al. (2021), the dynamical mass estimated from Equation 7 could be overestimated by a factor of ∼ 2.3.We show the M BH − M dyn relation after correcting for this factor using Equation 8 in Figure 13.In Figure 13, we present the corrected M dyn values for galaxies classified as rotation-dominated systems, encompassing both literature sources and our observations.Most of the quasars studied in this work, with the exception of J0910-0414, exhibit overmassive SMBHs when compared to local relations.The ratio M BH /M dyn is approximately three times higher than the local relation on average, spanning a broad range of 1.2-14.5.The literature sample, comprising ten rotationdominated quasar host galaxies, shows an even more pronounced BH mass excess on average (around 20 times higher).We observe that the disparities between these two samples arise from differences in the [C II] emission sizes, as depicted in Figure 6, and the slightly more rounded morphology of the quasar host galaxies studied in our work compared to those in the literature (e.g., Venemans et al. 2020;Neeleman et al. 2021).These differences can be attributed to the slightly higher resolution of literature observations, although with similar depth, leading to smaller [C II] sizes due to missing extended emission and more complicated and resolved structures.Considering both samples jointly, the deviation of the M BH − M dyn relationship for these highredshift quasar host galaxies from local relations could be a factor of approximately 10 (indicating black holes being overmassive by a factor of 10), with a broad range spanning from approximately 0.6 to 60.We acknowledge that the current M dyn measurements for both the literature sample and our work are subject to considerable uncertainties.Thus, we conclude that distant luminous quasars may indeed host overmassive black holes, but with substantial scatter and uncertainty in M dyn measurements.We further note that the M BH − M dyn correlation could be affected by the Lauer bias (Lauer et al. 2007), since overmassive objects are tend to be brighter and thus would be overrepresented in the selected sample (see Zhang et al. 2023, for a detailed discussion).

SUMMARY
In this work, we present a spatially resolved [C II] survey of 31 z ∼ 7 galaxies hosting luminous quasars with ALMA and NOEMA.The purpose of this paper is to describe the program design and to present the data and analyses of [C II] properties, continuum luminosity, star formation rate (SFR), the size and morphology of quasar host galaxies, as well as the first-order constraints on dynamical masses of the quasar host galaxies.We summarize our main findings below.
• We detect the [C II] emission of 26 quasar host galaxies with a [C II] luminosity range of L [CII] = (0.3−5.5)×10 9 L ⊙ , and detect the dust continuum of 27 quasar host galaxies with a far-infrared luminosity range of L FIR = (0.5−13.0)×10 12 L ⊙ .The SFR of these quasar host galaxies are in the range of 100-3000 M ⊙ yr −1 with a fiducial assumption of a graybody model with T = 47 K and β = −1.6.
• Both L [CII] and L FIR are correlated (ρ ≃ 0.4) with the quasar bolometric luminosity but with substantial scatter.
• The [C II]/FIR luminosity ratio ranges from ∼ 10 −4 to ∼ 10 −2 and anti-correlates with the FIR luminosity (L FIR ) or the FIR luminosity surface density (Σ FIR ).This so-called '[C II] deficit' is consistent with that found in lower redshift quasar host galaxies and infrared luminous galaxies.
• The majority of our quasar host galaxies are clearly resolved with a median diameter of ∼5 kpc, and the [C II] emission is usually less concentrated than the dust emission.
• The quasar host galaxies show diverse kinematics and morphologies, with ∼ 40% of them showing a clear velocity gradient, consistent with a rotating gas disk, while the rest of them show either dispersion-dominated compact morphology or disturbed kinematics.
• Our first-order estimates of the dynamical masses of the rotation-dominated quasar host galaxies give M dyn = (0.1 − 7.5) × 10 11 M ⊙ .Considering our findings alongside those of literature studies, we conclude that the majority of black holes at the centers of distant luminous quasars are overmassive compared to the local M BH − M dyn relation but with significant scatter, with the ratio between M BH and M dyn ranging from approximately 0.6 to 60 times that of local relations, and large uncertainties.
In the next steps, we will perform kinematic modeling of [C II] emissions to better constrain the dynamical mass of host galaxies.We will also stack the [C II] emissions to detect the low surface brightness emissions and search for outflow signatures which has been suggested to be strong in these luminous quasars by the strong blueshift of high-ionization UV broad emission lines and the high broad absorption line quasar fraction (Yang et al. 2021), as well as the presence of blueshifted [O III] lines (Yang et al. 2023).Furthermore, most of these quasars will be observed by ALMA mosaic observations from the ALMA Cycle 9 Large Program (2022.1.01077.L) which will enable us to characterize the Mpc-scale environment of these quasars.The majority of these quasars will also be observed by JWST which will allow us to detect the stellar light of these quasar host galaxies.Finally, combining the [C II] systemic redshifts with the quasar optical-to-infrared spectroscopy (Yang et al. 2020b(Yang et al. , 2021) ) will allow us to constrain the cosmic reionization and quasar life time by measuring the quasar proximity zone sizes and modeling damping wing profiles.

Figure 1 .
Figure 1.The redshift and M1450 distribution of quasars known at z ≥ 6.5.This [C II] survey aims to observe all z ≥ 6.5 quasars with M1450 < −25.0 (black dotted line).The quasars J0244−5008 (z = 6.724) and J0020−3653 (z = 6.834) were proposed in ALMA Cycle 5 by the discovery team (Reed et al. 2019) and were excluded from this survey.Another quasar J2356+0017 at z = 7.01 was discovered (Matsuoka et al. 2019b) after our survey begins.The ALMA observation of the gravitationally lensed quasar, J0439+1634, is also reported in this work.

Figure 2 .
Figure 2. ALMA (top 26 panels) and NOEMA (bottom 6 panels) thumbnails of z > 6.5 quasars showing the continuumsubtracted [C II] line maps in color and the continuum in contours.Each panel is 5 ′′ × 5 ′′ wide.North is up and east to the left.The solid black crosses on the [C II] undetected objects mark the expected emission position.The solid black contours mark the +3, 4, 5, 6, 7, 8, 9 isophotes for those objects with continuum peak detected at < 10σ, or 30%, 40%, 50%, 60%, 70%, 80%, 90% of the continuum peak emission for those objects with ≥ 10σ peak detections.The synthesized beam of the observations is shown in the bottom-left corner of each panel.

Figure 3 .
Figure 3. ALMA spectra of [C II] and the underlying continuum of the quasars in our sample (black lines).The blue lines are 1σ error vectors.The best-fit Gaussian line+flat continuum models are shown as solid red lines.Dotted lines indicate the best fits in those cases where no significant line emission is detected.

Figure 3 .
Figure 3. Continued.The bottom six spectra are NOEMA observations while the rest are ALMA observations.

Figure 4 .
Figure 4.The [C II] luminosity as a function of redshift for the quasars from our survey (orange points) and all other quasars with [C II] detections at z > 5.7 (blue circles).Our ALMA and NOEMA surveys detect the [C II] emission of 26 z > 6.5 quasars, representing 67% of all quasars at z > 6.5 with [C II] detections.Note that the brightest object is the lensed quasar J0439+1634 without correcting the lensing magnification.
F [CII] is the total line flux, ν obs is the observed frequency of the [C II] line, and D L is the luminosity distance.The [C II] line luminosities, L [CII], of all detected quasar host galaxies in our sample are listed in Table3and they span a range of L [CII] = (0.3 − 5.5) × 10 9 L ⊙ with the majority of them having L [CII] > 10 9 L ⊙ .To compare with observations available in the literature, we collected the properties of both [C II] and FIR continua of z > 5.7 quasars.There are 65 z > 5.7 quasars with [C II] and FIR continuum observations in the literature at the time of doing our analyses.Among them, 58 systems are detected in the [C II] line, with 13 at z > 6.5.Since most of the non-detections are caused by shallower observations, we only include the [C II]-detected sources in the following analyses.Figure 4 shows the L [CII] distribution of quasar hosts detected from this survey and all other known z > 5.7 quasars with published L [CII] measurements.Our [C II] survey

Figure 5 .
Figure 5.The relations between quasar bolometric luminosity and the [C II] luminosity (lef t) as well as the far-infrared luminosity (right) at z > 5.7, respectively.The Spearman test gives a correlation coefficient of ρ = 0.43 and a chance probability of p = 1.7 × 10 −4 for the L [CII] − L bolometric relation, and ρ = 0.48, p = 7.4 × 10 −6 for the LFIR − L bolometric relation, respectively.A moderate correlation exits for both relations, though both relations have significant scatter.

Figure 6 .
Figure 6.Comparison of the sizes of the continuum and [C II] emissions.The sizes are measured after deconvolution with the beam.The black dashed line denotes the one-toone relation.In general, the continuum emission is more compact than the [C II] emission.

Figure 7 .
Figure 7.The relation between the [C II] luminosity and the far-infrared luminosity.The star formation rates estimated from the [C II] luminosity (De Looze et al. 2014) and from the far-infrared luminosity (Murphy et al. 2011) are also shown.The black dashed line is the one-to-one case for the SFR [CII] and the SFRFIR.

Figure 8 .
Figure 8.The [C II]-to-FIR luminosity ratio as a function of FIR luminosity (left) and the FIR luminosity surface density, ΣFIR (right).The [C II]-to-FIR luminosity ratios of these high redshift quasars span a wide range (nearly two orders of magnitude).In the right panel, we also plot the best fits found by Lutz et al. (2016) and Díaz-Santos et al. (2017) from local-infrared luminous galaxies.The [C II] deficit is clearly seen in these plots but the relations have substantial scatter.

Figure 9 .
Figure 9.The [C II] based SFR is consistent with FIR based SFR after considering the FIR surface brightness of the quasar host galaxies.The x-axis in this plot is the SFR estimated from the LFIR, same as Figure 4.The y-axis is the SFR estimated from L [CII] using the relations from Herrera-Camus et al. (2018) after considering the FIR surface brightness of the quasar host galaxies.

Figure 10 .
Figure10.High-resolution ALMA moment 0 (left), 1 (middle), and 2 (right) maps of 12 objects with clear [C II] velocity gradients.Each panel is 15 × 15 kpc 2 wide.North is up and east to the left.In the moment 0 map, the solid black contours mark the +3, 4, 5, 6 isophotes for those objects with continuum peak at < 7σ, or the 50%, 60%, 70%, 80%, 90% of the continuum peak emission for those objects with ≥ 7σ peak detections.The synthesized beam of the observations is shown in the bottom-left corner of each panel.The solid contours and dashed contours in the moment 1 maps mark the positive velocities and negative velocities, respectively.The contour levels for both positive velocities and negative velocities are in steps of 20 km s −1 .

Figure 11 .
Figure 11.Same as Figure 10, but for 11 galaxies without clear velocity gradients.observed with ALMA (J0213-0626, J0319-1008, and J1129+1846) have relatively low significance detections, hence we only focus on the 23 objects observed with ALMA.The [C II] intensity map is collapsed from −1.4σ line to +1.4σ line of the emission-line data cube.We then create a mask image by selecting all pixels with S/N> 3 in the intensity map.We finally create a velocity map (moment 1) and a velocity dispersion map (moment 2) of each target using the mask image created from previous step and selecting frequencies between −2.0σ line and +2.0σ line with the immoments task within CASA package.The moment maps are

Figure 13 .
Figure13.Same as the right panel of Figure12, but for the dynamical mass estimated from Equation8.Most of the galaxies fall above the locally derived MBH-M bulge correlation (e.g.Kormendy & Ho 2013) with the black holes being over-massive by a factor of ∼ 10, although with substantial scatter.

Table 1 .
The basic information and observing log of 31 z > 6.5 quasars studied here.

Table 2 .
The measured frequency, flux and width of the [C II] line and the flux of underlying continuum.