© ESO 2020

1. Introduction

Supermassive black holes (BHs; M_BH ∼ 10⁶ − 10¹⁰ M_⊙) are believed to reside at the centre of all nearby galaxies and are likely the relics of a past quasar (QSO) activity (e.g. Soltan 1982; Hopkins et al. 2008). Such BHs have likely played a key role in shaping galaxies during their assembly at early epochs, with the implication that BH growth and galaxy formation are closely linked (Heckman & Best 2014).

The discovery of the strong correlations (in the local Universe) between the mass of the central black hole (M_BH) and the physical properties of host galaxies (e.g. stellar velocity dispersion of the bulge stars, mass of the bulge, etc.; Tremaine et al. 2002; Häring & Rix 2004; see also Kormendy & Ho 2013 for an extensive review and references therein) has been one of the most significant breakthroughs of the past decades and represents a key building block for our understanding of galaxy formation and evolution across the cosmic time. In the framework of co-evolution between BHs and their host galaxies, the observed local relations are believed to arise from the balance between the energy released by the active galactic nucleus (AGN), which generates galactic-scale outflows expelling gas from the galaxy, and the gravitational potential that keeps the galactic system bound. According to current galaxy evolution models (see e.g. Lamastra et al. 2010; Sijacki et al. 2015), AGN are able to regulate the star formation activity in the host and constrain both the final stellar mass and dynamical properties of the galaxy (e.g. Di Matteo et al. 2005; Menci et al. 2008; Hopkins et al. 2008; Kormendy & Ho 2013). Therefore, investigating the onset of BH-galaxy relations at high redshift is fundamental to exploring the interplay between BH accretion and star formation activity in the host galaxies, and to constrain, accordingly, galaxy formation and evolution models.

In this work, we focus on the relation between BH mass and that of the host galaxy (M_BH−M_gal relation). The latter has been widely sampled for active and quiescent galaxies in the local Universe (z <  1), indicating that BH mass is a defined fraction of the bulge stellar mass (M_BH ∼ 10⁻³M_gal, e.g. Marconi & Hunt 2003; Häring & Rix 2004). More recently, several groups (Treu et al. 2004, 2007; Walter et al. 2004; Peng et al. 2006a,b; Shields et al. 2006; Woo et al. 2006, 2008; Ho 2007; Decarli et al. 2010; Merloni et al. 2010; Wang et al. 2010; Bennert et al. 2011; Canalizo et al. 2012; Targett et al. 2012; Bongiorno et al. 2014) attempted to sample this relation beyond the local Universe, showing that there are indications for a possible evolution with redshift. In particular, these works suggest a parameterisation of the ratio Γ = M_BH/M_gal as a function of redshift, Γ ∝ (1 + z)^β. The published values of β span the range 0.7−2 (McLure et al. 2006; Bennert et al. 2010, 2011; Decarli et al. 2010; Merloni et al. 2010) with the implication that, at higher redshifts, galaxies host black holes that are more massive than the local counterparts (e.g. a factor of ∼7 at z ∼ 3; Decarli et al. 2010). Therefore, during the competitive accretion of matter from the galactic halo that occurred at early epochs, black hole growth possibly must have preceded that of the host galaxy (e.g. Decarli et al. 2010; Bongiorno et al. 2014; Wang et al. 2016).

However, the aforementioned results are affected by observational biases and instrumental limits. The selection of host galaxies revealed at high redshift (z >  3), is driven by AGN luminosity, so more massive black holes are preferably selected (Lauer et al. 2007; Vestergaard et al. 2008; Volonteri & Stark 2011; Portinari et al. 2012; Schulze & Wisotzki 2014; Volonteri & Reines 2016). Then, in these sources, the luminosity of the central region overwhelms the emission from the host galaxy, and the disentangling of the two components is challenging even with high-resolution observations. Since the galaxy stellar mass estimates used to derive the M_BH−M_gal relation are based upon photometric methods, they are significantly contaminated by light from the central non-stellar source, and are thus very uncertain. Finally, since the average gas fraction of galaxies increases with the redshift (e.g. Magdis et al. 2017; Tacconi et al. 2018), primordial galaxies may not yet have converted a large fraction of their gas into stars, therefore their stellar mass content may not be a reliable tracer of the total mass (but there is also evidence of luminous QSOs with low gas fractions possibly related to the effect of an AGN-driven feedback mechanism, see e.g. Carniani et al. 2017; Kakkad et al. 2017; Brusa et al. 2018; Cresci & Maiolino 2018; Perna et al. 2018). The galaxy’s capability of retaining its gas under the influence of AGN activity, is indeed determined by the gravitational potential of the whole galaxy traced by the total (dynamical) mass.

The recent advent of ALMA (Atacama Large Millimeter and Sub-Millimeter Array) opened a new era of cold gas observations. Thanks to its unparalleled capability in terms of sensitivity, signal-to-noise ratio (S/N), and angular resolution, it is now possible to spatially resolve the gas kinematics in quasar host galaxies up to the higher redshifts targeting the brightest line emission of the cold gas, such as [CII]_158 μm or CO rotational line transitions with sub-mm spectroscopic observations (see Carilli & Walter 2013; Gallerani et al. 2017, for a comprehensive review). In fact, the emission of radio-quiet AGN in the sub-mm band is dominated by the cold gas mass and the dust continuum in their hosts, thus allowing observations that are not affected by the non-stellar emission of the central source. Therefore, thanks to the efforts of many groups, ALMA has made it possible to trace the BH-galaxy relation at very high redshift using dynamical mass estimations of host galaxies (e.g. Wang et al. 2013, 2016; Willott et al. 2013, 2015a; Venemans et al. 2012, 2016, 2017a; Decarli et al. 2017; Trakhtenbrot et al. 2017; Feruglio et al. 2018). The dynamical masses provided in these works are estimated assuming rotating disc geometry and by simply combining the full width half maximum (FWHM) of the observed line emission, the observed size of the emitting region, and the inclination angle of the galaxy disc with respect to the sky plane. However, it is hard to test the basic assumption that the cold atomic/molecular gas of the galaxy is a rotating disc. Furthermore, the disc inclination is calculated from the observed morphology by using the axial ratio of the flux map, and is thus affected by significant uncertainties.

In this work, we studied a large sample (∼70) of 2 <  z <  7 quasars observed by ALMA targeting the [CII]_158 μm atomic fine-structure line or the CO rotational line emission, which we exploited in order to trace the morphology and kinematics of quasar host galaxies. Overall, we identified ordered rotational motion in a sample of ten quasars (for which high S/N allowed a spatially resolved analysis). By carefully modelling the kinematics with rotating discs, we were able to measure their host galaxy dynamical mass (M_dyn), at variance with previous work where rough estimates are usually adopted. Our dynamical mass measurements, combined with M_BH estimates obtained from the literature allowed us to trace the evolution of the M_BH−M_dyn relation and to study the trend of Γ = M_BH/M_dyn across the cosmic time.

The paper is organised as follows: in Sect. 2, we outline our starting sample and the data reduction performed on the raw data. In Sect. 3, we illustrate the methods of data analysis to retrieve the information on the morphology and kinematics of the host galaxies. In Sect. 4, we present the kinematical model used to measure the galaxy dynamical mass. In Sect. 5, we obtain the BH masses from the literature and from LBT data. In Sect. 6, we recap the different sub-selections of the starting sample that occurred during this work. In Sect. 7, we compare our dynamical mass estimates with previous similar studies and discuss the uncertainties on our measurements. Then, we investigate limits of validity of the assumptions. In Sect. 8, the M_BH−M_dyn relation and the trend of M_BH/M_dyn ratio across cosmic time are presented. Then, in Sect. 9, we discuss our results and compare them with previous works. We also examine how possible additional uncertainties and biases could affect the results both from observational and theoretical points of view. Then, we compute the virial masses of our final sample, and we compare them with our dynamical mass estimates. Finally, in Sect. 10, we draw our conclusions.

Throughout the paper, we assume a standard ΛCDM cosmology with H₀ = 69.3 km s⁻¹ Mpc⁻¹, Ω_m = 0.287, Ω_Λ = 1 − Ω_m from Hinshaw et al. (2013).

2. ALMA data selection and reduction

We started by collecting all [CII]_158 μm and CO(J → J − 1) (rotational quantum number J = 4, 3) observations of z >  1.5 QSOs on the ALMA data archive public as of June 2017 for a total of 72 QSOs in the redshift range 1.5 <  z <  7.1. Different ALMA bands were involved according to the targeted atomic/molecular transition and the redshift of the sources. The collected data were calibrated using the ALMA pipeline in the Common Astronomy Software Applications, CASA (McMullin et al. 2007), by executing the appropriate ALMA calibration scripts corresponding to each specific observation. Continuum images were produced for each quasar from the calibrated visibilities, by combining the line-free channels from all spectral windows in multi-frequency synthesis mode using the CASA task tclean and briggs weighting scheme (with robustness parameter R = 0.5) to maximise both the signal-to-noise ratio and angular resolution. The line-free channels were determined by inspecting the visibilities in all the frequency sidebands. For those quasars in which the FIR line was not detected, we selected the line-free channel by adopting a line width of 300 km s⁻¹ and the redshift from literature.

These same channels were also used to produce a UV plane model by fitting the continuum emission with a zeroth order polynomial¹ that was then subtracted from the spectral windows containing the line using the CASA task uvcontsub. The continuum-subtracted line visibilities were then imaged using tclean. In order to recover all the information within the resolution element, the pixel size was commonly set to ∼B_min/7, where B_min is the minor FWHM of ALMA’s synthesised beam. Therefore, we obtained cubes with a typical pixel size of 0.025″ − 0.05″ and with a spectral bin width set to 40−70 km s⁻¹. Self-calibration was attempted but showed no additional improvement for almost all observations and was not used for the final cubes. Finally, both continuum images and the line cubes were corrected for the primary beam response.

Among these observations, we selected the cubes in which the line detection was significant (≳3σ). This first selection reduced the sample to 32 QSOs at 2.2 <  z <  7.1 on which we performed all of the analyses described in the following sections. Different sub-selections occurred at each step of the analysis (see Sect. 6) and the final sample is composed of only ten sources, for which we obtained constraints on the host galaxy’s dynamical mass. We thus picked deeper observations from the archive for this final sub-sample of sources that became public while the work was in progress (by the end of February 2018). In Table 1, we list the starting sample of 32 objects and the characteristics of the observations, including the aforementioned deeper observations for the final sub-sample. The distribution of redshifts of our quasars is illustrated in Fig. 1.

Fig. 1. Redshift distribution of 32 QSOs listed in Table 1. The histogram colours indicate the emission line detected in the ALMA data. The figure also shows the sample cumulative redshift distribution (dotted line and right axis). Individual sources are marked with blue circles (sizes are proportional to the ALMA synthesised beam FWHM).

3. Methods of data analysis

Our goal is to measure the dynamical mass of our sample of host galaxies (listed in Table 1) by modelling the gas kinematics as traced by [CII]_158 μm or CO line emission with rotating discs. Therefore, in order to obtain the kinematical maps, we performed a spaxel-by-spaxel fit of the emission line profile by adopting a single Gaussian model with three free parameters: the amplitude A, the central frequency ν_obs, and the standard deviation σ. For this purpose, we designed an algorithm to achieve a robust residual minimisation in each pixel. Since the beam smearing affects the observed emission, we expect the spatial shape of the line to change smoothly from one pixel to the adjacent one, and the signal-to-noise ratio to decrease as a function of the distance from the centre of the galaxy. The underlying idea of the procedure is the subsequent performance of the line fit in all spaxels starting from the central pixel and moving away following a spiral-like path. The basic operations are: (1) performing a 2D Gaussian fit on the continuum image and defining the central spaxel; (2) extracting the spectrum from the central spaxel and computing a 1D Gaussian line fit in which the starting points are properly chosen by inspecting the line shape; (3) following a spiral-like path to select the next spaxel, extracting the spectrum, and performing a 1D Gaussian line fit by using the best-fitting results from the neighbour spaxels as starting points for the spectral fit; (4) continuously repeating step3 for the consecutive spaxel until the end of the spiral-like path. We used the minimum chi-square method to estimate the best-fitting parameters.

The result of the fit in each pixel is accepted or rejected on the basis of criteria illustrated in Sect. 3.2, while the stopping criteria to break the entire fitting procedure can be fixed by setting the dimension d of the spiral path, that is the distance from the central pixel. In the case of our datacubes, a typical value of d ∼ 20−25 pixels (∼0.5″ − 1.25″ depending on the pixel size, see Sect. 2) turned out to be adequate to fit the line throughout the emitting region with a total of ∼1200−2000 pixels analysed for each source. This fitting strategy enables a more robust minimisation compared with using a unique set of initial guess parameters for all the pixels, thus avoiding numerical problems arising from incorrectly chosen starting-points. Finally, we retrieved the information regarding the line together with the uncertainties on each spectral channel of the cube by measuring the rms of the noise (r_ν) over a wide spatial region where no emission is detected.

3.1. Integrated spectra and derived quantities

We obtained the integrated spectra of all the sources by adding all the fitted spectra in spaxels selected based on criteria illustrated in Sect. 3.2 (e.g. Fig. 2, see also Appendix E). Then, the resulting integrated spectrum was fitted using a Monte Carlo method in order to estimate the redshift uncertainty. Firstly, we collected a large number (e.g. 2000) of different integrated spectra obtained by adding a random value extracted from a normal distribution defined by a zero mean, and a standard deviation equal to the corresponding rms in that channel (r_ν) to each channel of the original spectrum. Then, we performed the fit of each spectrum with a single Gaussian, and we estimated the redshift of the line as ν_obs = ν_rest/(1 + z), where ν_obs is the mean of the Gaussian model and ν_rest is the line rest-frame frequency. Finally, all the estimates of z obtained with this method were histogrammed and its distribution was fitted with a Gaussian model. We finally assumed the mean and the standard deviation of the best-fit model as the best value of redshift and its uncertainty, respectively. In addition, the fit of integrated spectra allowed us to determine the line FWHM and flux. In Appendix A, we use these quantities to derive the line luminosity, the [CII] mass (M_[CII]), the total gas mass (M_gas) and the star formation rate (SFR) of the quasar host galaxies.

Fig. 2. Integrated spectrum of source SDSS J0923+0247. Top panel: the observed data is shown in light yellow with error bars in grey (rms in each channel). The red solid curve is the best-fit Gaussian model. The velocity scale (top axis) has as its referecence the central frequency of the best fit. Bottom panel: the fit residuals (model-data), the grey filled area shows the rms along the spectral axis.

3.2. Flux, velocity, and velocity-dispersion maps

The cube fitting procedure provides the best-fit values of the Gaussian parameters (A, ν_obs, σ) in each pixel. We used these values to obtain the line-integrated velocity and the velocity-dispersion maps along the line of sight (LOS).

In order to produce the maps, among all the spaxels in which we performed the line fit, we selected those satisfying the following conditions: (1) the peak of the best-fitting Gaussian is ≥1.5 × rms in the corresponding channel; (2) the percentage relative error on the flux value is ≤50%.

Condition 1 represents the signal-to-noise cut-off we used to reject pixels in which the line emission is not clearly detected. However, in case of poor signal-to-noise ratio, the fit process possibly fails, resulting in a bad Gaussian model for which condition 1 could be still satisfied. Therefore, we also imposed condition 2 in order to avoid this kind of effect and to consequently reject the corresponding pixels when producing the maps.

We also manually masked bad pixels far away from the galaxy centre, which are clearly associated with spikes of noise. Finally, we obtained error maps using the uncertainties on the best-fit Gaussian model parameters of each pixel. As an example, in Fig. 3, we report the maps obtained from the continuum-subtracted cube of SDSS J0923+0247 (see also Appendix E).

Fig. 3. Observed maps of SDSS J0923+0247. From left panel to right: flux map, velocity map, and velocity-dispersion map along the line of sight, respectively. At the bottom-left corner of each panel, we report the ALMA synthesised beam FWHM. The coordinates indicate offsets with respect to the map centre.

3.3. Red and blue residuals maps

The angular resolution may not be high enough to spatially resolve the rotation of the emitting gas in host galaxies. In order to assess if the kinematics are spatially resolved or not, we performed the same analysis computed by Carniani et al. (2013) for ALMA [CII] observations of a QSO at z = 4.7. We replicated the spaxel-by-spaxel fit of continuum-subtracted cubes with a single Gaussian component, using the amplitude (A) as the only free parameter, and by fixing the values of ν_obs and σ to the best values obtained from the fit of the integrated spectrum (see Sect. 3.1). Then, we computed the residuals of fits in each channel, which is the data-model, and we obtained two maps by collapsing all channels in the blue-shifted and red-shifted (with respect to the central frequency) side of the residual spectrum. If the kinematics are consistent with a spatially resolved rotating disc, we expect the blue and red residual maps to show two symmetric lobes: a positive and negative lobe at the opposite side with respect to the map centre (e.g. as shown in Fig. 4). Otherwise, if rotation is not resolved then we expect a random distribution of negative and positive residuals on both the maps. After performing this test on the 32 objects listed in Table 1, we conclude that 14 of them ( ∼ 45%) show no evidence of spatially resolved kinematics. We therefore excluded them from the final sample (see Sect. 6 for a summary of sample sub-selections). In Sect. 7.4, we investigate biases possibly occurring while excluding these objects.

Fig. 4. Red and blue residual maps of SDSS J0923+0247. This test reveals spatially resolved kinematics consistent with what we would expect from a rotating disc.

4. Kinematical modelling

In Sect. 3, we obtained all the necessary information about galaxy morphology (line-integrated maps) and kinematics (velocity and velocity-dispersion maps). In order to measure the dynamical masses of the host galaxies, we designed a kinematical model to perform a 2D fit of the maps. Therefore, we assumed that:

1. The observed line emission (i.e. [CII]_158 μm or CO transition) traces cold gas distributed in a rotating thin disc.

2. The gas mass surface density Σ(r), is an exponential distribution that also tracks the distribution of surface brightness I(r), that is:

   (1)

   where I₀ is a normalisation constant and R_D is the scale radius.

3. The galaxy stellar mass is distributed as the gas mass component (Eq. (1)).

4. The contribution of the dark matter is negligible.

Under these assumptions, Freeman (1970) showed that the corresponding circular velocity is given as:

(2)

where ℐ and 𝒦 are the modified Bessel functions evaluated at y = r/2R_D, and Σ₀ is the normalisation constant of the mass distribution that accounts for both gas and stars contribution.

The total mass of the disc, which is the dynamical mass of the galaxy (M_dyn), is thus obtained by integrating the mass surface density over all the radii; . By inserting this expression in Eq. (2), we can relate the total dynamical mass with the velocity curve: V²(r) = 2(M_dyn/R_D)Gy²[ℐ₀(y)𝒦₀(y)−ℐ₁(y)𝒦₁(y)]. Therefore, by estimating R_D from the flux map, we can infer the galaxy dynamical mass by performing a 2D fit of the velocity field.

4.1. Details on the kinematical model and the strategy of the analysis

The kinematical model is calculated using Monte Carlo methods. At first, the 3D space is randomly filled with N ≫ 1 point-like sources uniformly distributed in a thin disc. Each source represents a “cloud” that contributes with a unit of flux in the computation of the total observed flux. Then, the 3D disc model is projected on the sky plane and convolved with the appropriate instrumental point spread function (PSF) and line spread function (LSF) of the observation. Thus, the flux map, the flux-weighted velocity map and the velocity-dispersion map are obtained through 2D-weighted histograms by using the pixel size of the corresponding observed map as the bin width. By properly choosing the weights, we can set the flux contribution of each cloud forming the model in order to reproduce any brightness (density) and velocity profiles. We set the weights in order to create an exponential thin disc defined by Eqs. (1) and (2). As an example, in Fig. 5, we show the simulated flux, velocity and velocity-dispersion maps obtained with our kinematical model in three different geometrical configurations.

Fig. 5. Simulated maps obtained with the kinematical model for a galaxy thin disc with an exponential brightness profile defined by a scale radius RD = 0.125″ and dynamical mass Mdyn = 5.0 × 1010 M⊙. The flux, velocity, and velocity-dispersion maps along the LOS are indicated by F(xp, yp′), V(xp, yp′), and σ(xp, yp′), respectively. From top to bottom panel: disc inclination and position angle (β, γ) are, respectively, equal to (40, −45) deg, (40, 0) deg, and (60, 45) deg. The other parameters defining the model are: the FWHM of the synthesised beam (0.4″ × 0.275″), the position angle of the beam (BPA = −50 deg), the FWHM of the LSF (σLSF = 20 km s−1), the bin size (0.05″ × 0.05″), and the angular radius of the disc (R = 1.25″). The number of clouds used in this model is Np = 7.5 × 106.

In order to recover the galaxies dynamical masses, we basically adopted the same method used in Cresci et al. (2009) and Carniani et al. (2013). We first performed a 2D fit of the observed flux map using a thin disc model with an exponential brightness profile, and we recovered the best value of scale radius R_D (see Eq. (1)); then, by using the resulting R_D value, we computed the velocity field of our disc model accordingly with Eq. (2), and we performed a 2D fit of the LOS velocity map, thus recovering the best estimate of M_dyn.

4.2. Estimation of R_D: 2D fit of the flux maps

Following the method illustrated in the previous section, we first estimated the scale radius R_D on the 18 flux maps of the sources with spatially-resolved kinematics. The typical angular extension of the observed maps is ∼1″, with a pixel size depending on the beam size of the observation (see Sect. 2) resulting in a typical map size of ∼15−20 pixels in linear diameter. We thus generated simulated maps using a 3D disc model with radius of R = 20 pixels filled by N_p = 5 × 10⁶ clouds. These values turned out to be the best compromise to smooth the stochastic oscillations of the cloud’s numerical density and to avoid spurious numerical effects at the boundary of the model, while simultaneously keeping the computational time relatively short.

The 3D disc model is then projected on the sky plane, and, according to Eq. (1), the observed image of the simulated flux map depends on the normalisation constant I₀, the scale radius R_D, and on the geometrical parameters: the coordinates of the map centre (x₀, y₀), the inclination with respect to the sky plane (β) and the position angle of the line of nodes (γ) measured clockwise with respect to the east. We note that, for the purposes of this work, we were not interested in the physical value of I₀.

The aforementioned parameters were variable during the fit procedure. Thus, to retrieve their best estimations, we carried out the 2D map fit using the cap-mpfit PYTHON procedure part of the pPXF package by Cappellari & Emsellem (2004) based on minpack-1 (Moré et al. 1980), performing a Levenberg-Marquardt least-squares minimisation between the data and the model. For each map, the best model minimised the following function:

(3)

where , F and σ_F are, respectively, the observed flux map, the flux model map, and the flux error map. We note that, in addition to the standard χ² function, we inserted a penalty term 𝒫 = 5 × N_out in Eq. (3), where N_out is the number of pixels defined in the data but not in the model. Indeed, the sum in Eq. (3) is computed only taking into account the pixels (x_p, y_p′) on the observed map in which the model is defined. Therefore, unless the model is defined in all the pixels in which the observed data are present, during the minimisation process, the penalty term 𝒫 ensures the adequate penalisation of configurations for which the model cannot reproduce the data in all the points (e.g. since the disc model is thin, completely edge-on disc configurations are highly unlikely unless the PSF is large enough).

In order to obtain a robust χ² minimisation avoiding the convergence to a possible local minimum, the 2D flux map fit is performed multiple times by fixing the disc inclination angle (β) to 5 deg and up to 90 deg with a step size Δβ = 5 deg. For each value of β, the χ² function in Eq. (3) is minimised with respect to the free parameters [I₀, R_D, γ, x₀, y₀]. At the end of each step, we retrieve the minimum of χ²(β) function and the corresponding set of best values of free parameters. We then used them as starting points for the next step. Once the minimisation is performed for all the inclination values in the range [5, 90] deg, we sampled the curve of the minima as a function of the disc inclination angle, which is . Finally, by finding the absolute minimum of , we retrieved the best set of [I₀, R_D, β, γ, x₀, y₀]. Following the method illustrated above, we estimated the best value of R_D measured in arcseconds for all 18 objects indicated in Table 1 with flag [rot] = “r”. Finally, we computed R_D values in physical length, using the redshift estimates obtained in Sect. 3.1.

In the next section, we use the R_D estimates to compute the kinematical model in order to perform the 2D fit of the velocity fields. For this purpose, we use the curve as a prior knowledge on disc inclination angle. As an example, in Fig. 6, we show flux map modelling results and the correspondent curve for SDSS J0923+0247. The R_D values are listed in Table 2 for those objects with constraints on dynamical mass (see also Appendix E).

Fig. 6. Kinematical modelling performed on SDSS J0923+0247. Upper panels [a.1] and [a.2]: fit result of the flux map, bottom panels [b.1] and [b.2]: result of the 2D velocity field fit. [a.1]: curve of χ2 minima as a function of the disc inclination (see Sect. 4.2 for details). The magenta hexagon indicates the absolute minimum (=35 deg) of χm(β). The green area (see Sect. 4.3 for the definition) indicates the allowed inclination values used as a prior in performing a 2D velocity field fit. [a.2]: 2D best-fit model of the flux map. From left to right, we report the observed map, the model, and the residuals ([data-model]/error). Residuals values are color-coded according to the horizontal colorbar displayed at the bottom-left corner of the panel. [b.1]: the posterior probability distributions of the free parameters in 2D velocity field fits retrieved with the MCMC algorithm with the best values and their uncertainties. [b.2]: 2D best-fit model of the kinematical map. From left to right, we report the observed velocity map, the model, and the velocity curves extracted from a long-slit of two pixels in width aligned with the line of nodes. The slit is superimposed on maps; red circles and solid black lines in the right panel indicate the observed flux-weighted velocity values in each bin of the slit and the model, respectively.

4.3. Estimation of M_dyn: 2D fit of the velocity maps

In order to estimate the dynamical masses of the quasars sample, we performed the 2D fit of LOS velocity maps. At high-z, uncertainties on the dynamical mass estimates are mainly driven by the poor angular resolution of observations. As the integrated flux map of line emission, even the velocity maps are affected by beam smearing, thus introducing additional uncertainties in the fitting parameters (as also pointed out by other authors, see e.g. Bosma 1978; Begeman 1987; de Blok & McGaugh 1997; O’Brien et al. 2010; Swaters et al. 2000, 2009; Epinat et al. 2009, 2010; Carniani et al. 2013; Kamphuis et al. 2015). This effect leads to the disc inclination angle and the dynamical mass becoming almost degenerate parameters, meaning the observed velocity field can be similarly reproduced by different couples (β, M_dyn) with similar M_dynsin²β, thus providing very near values of χ² function (see also e.g. Epinat et al. 2010). In Fig. 7, we show the effect of the beam smearing on the iso-velocity curves of simulated velocity fields.

Fig. 7. Simulated velocity field corresponding to four different couples of inclination angle β and dynamical mass Mdyn = 5.0 × 1010/​​sin2β. Increasing the dimension of the beam FWHM (shadowed area), the iso-velocity curves are increasingly smoothed, and different configurations of the disc appear almost indistinguishable. The simulated maps are obtained with the kinematical model for a galaxy thin disc with RD = 0.125″, γ = −90 deg. The other parameters defining the model are set as equal to the model in Fig. 5.

Consistently with the method illustrated in Sect. 4.1, we performed fits of velocity maps using thin rotating disc models defined by exponential mass distributions with the R_D values estimated in Sect. 4.2. In order to retrieve the best-fitting model and to estimate the parameter uncertainties, we used the PYTHON affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013). For this purpose, we defined the likelihood function of a model, given the data, as:

(4)

where , V, and σ_V are the observed velocity map, simulated velocity map, and velocity error map, respectively.

The simulated field depends on geometrical parameters, the scale radius, and the dynamical mass (see Eq. (2)). However, instead of M_dyn, we used μ = log(M_dynsin²β) as free parameter, since it is decoupled from disc inclination, thus it uniquely determines the intrinsic velocity field. Here, we also took into account the systemic velocity of the galaxy, which is V_sys, as an additional free parameter. Furthermore, defining the likelihood function in Eq. (4), we assumed Gaussian priors for the coordinates of the galaxy centre with maximum probability corresponding to (x₀, y₀) obtained from the flux map modelling, as illustrated in Sect. 4.2. Here, Δx₀ and Δy₀ are the standard deviations assumed equal to 0.1 pixel. In addition, we assumed a box-like prior on the position angle of the disc (−180 ≤ γ(deg)≤180) and on the dynamical mass (M_dyn >  0). Finally, we accounted for prior knowledge on disc inclination from morphology by imposing that, during the fitting process, where with the absolute minimum of curve resulting from the fit of the flux map (see Sect. 4.2).

We thus maximised Eq. (4) using [V_sys, μ, γ, sin β, x₀, y₀] as free parameters, and we recovered their posterior probability distributions². Finally, the best values of parameters and related uncertainties were estimated by computing the 16th, 50th, and 84th percentile of the distributions. As an example, Fig. 6 shows the kinematical modelling of SDSS J0923+0247 (see also Appendix E).

Due to the poor spatial resolution of some observations, we successfully constrained the disc inclination and the dynamical mass for only 13 objects of the sample (see Table 2). We note that for two of them we find bimodal distributions for sin β and μ, thus not permitting us to define unique values.

5. Determination of the black hole masses

Currently, the only possible technique to carry out black hole mass estimates at high-z is the use of single epoch (SE) virial relation, which combines the FWHM or the line emission that originated in the broad line region (BLR) of quasar, with the continuum luminosity emitted from the BH accretion disc (e.g. McLure & Dunlop 2004; Vestergaard & Peterson 2006; Vestergaard & Osmer 2009; but see also e.g. Trevese et al. 2014; Grier et al. 2019). This approach assumes that the BLR is virialised and that there is a tight relation between the BLR radius (R_BLR) and the continuum luminosity of the AGN (L_AGN) (e.g. Kaspi et al. 2005; Bentz et al. 2006, 2009). Under these assumptions, L_AGN and the FWHM of the broad emission lines are used as proxies for R_BLR and virial velocity, respectively.

To date, thanks to the effort of various groups (e.g. McLure & Jarvis 2002; McLure & Dunlop 2004; McGill et al. 2008; Wang et al. 2009; Shen et al. 2011; Denney 2012; Park et al. 2013; Coatman et al. 2017), many relations have been calibrated by employing different broad lines in order to infer the M_BH and by assuming that such BH-mass estimates are in agreement with reverberation mapping BH masses (Vestergaard & Peterson 2006), which, in turn, are in agreement with local BH-galaxy scaling relations for normal galaxies (Onken et al. 2004). This is motived by the unknown geometry and kinematics of the BLR (e.g. McLure & Dunlop 2001; Onken et al. 2004). In high-redshift quasars, the atomic transitions of MgII and CIV are the most common and brightest BLR lines that are observed in optical (rest-frame UV) spectra, and thus they are extensively used as virial mass estimators with the corresponding continuum luminosity measured by convention at 3000 Å and 1550 Å for MgII and CIV, respectively. However, the reliability of CIV line is still strongly debated. Firstly, the CIV scaling relation is based on very few measurements (Kaspi et al. 2007; Saturni et al. 2016; Park et al. 2017); secondly the CIV line is often associated with broad and blueshifted wings likely resulting from outflows (e.g. Richards et al. 2011; Denney 2012) that may affect the measurement of the line width biasing the BH-mass estimates. We note that the aforementioned calibrations have intrinsic uncertainties of ∼0.3−0.4 dex (see, e.g. Vestergaard & Peterson 2006; Denney 2012; Park et al. 2017) that are usually larger than the errors associated with line width and flux measurements.

5.1. Black hole masses from the literature

In this work, we adopted a unique SE virial relation to estimate the BH masses of our sample homogeneously. In detail, we used the relation by Bongiorno et al. (2014), which was calibrated by assuming the BH-galaxy scaling relations by Sani et al. (2011). The latter is consistent with the relation used as a z = 0 reference for studying the redshift evolution (e.g. Kormendy & Ho 2013; de Nicola et al. 2019):

(5)

Thus, where available, we retrieved MgII FWHM and the continuum luminosity estimates from the literature and, if they were unavailable, we assumed M_BH estimates as provided in the literature.

In summary, we recovered BH masses from the literature for everything except PSO J308−21 and PSO J183+05, using Hβ, MgII, and CIV BLR lines. For SDSS J0129−0035 and SDSS J2054−0005, the black hole masses were estimated from the bolometric luminosity (L_bol) assuming Eddington accretion (L_bol/L_Edd = 1). We note that this assumption is supported by some evidences showing that black holes at z ≳ 6 accrete matter at a rate comparable to the Eddington limit (De Rosa et al. 2014; Mazzucchelli et al. 2017). The data are listed in Table 3. The BH masses computed in this work using Bongiorno et al. (2014) calibrations are a factor ∼2 smaller than those reported in literature using different calibrations (see references in Table 3 for full details). However, both estimates are consistent within the typical uncertainties (∼0.4 dex).

5.2. Black hole mass from LBT data

The study of the redshift evolution of M_BH−M_dyn relation can be severely affected by reliability of M_BH estimates; in particular, the available spectroscopic information for our sample did not allow us to derive M_BH measurements with a unique method. We therefore started an observational campaign with the LBT (Large Binocular Telescope), of those sources with estimated M_BH assuming L_bol/L_Edd = 1 (SDSS J0129−0035, SDSS J2054−0005, and SDSS J2310+1855), and of two additional targets without previous M_BH estimates from the literature (PSO J308−21 and PSO J138+05; even though dynamical masses of these two sources are tentative). We thus obtained NIR spectra of the quasars with LUCI (LBT Utility Camera in the Infrared) targeting the CIV line and the adjacent continuum, which are redshifted into the zJ filter (central λ = 1.17 μm), for everything except PSO J183+05. In fact, for all these sources, the MgII line falls in a spectral region with a very low atmospheric transmission. For PSO J183+05, instead, we targeted the BLR MgII line, which is redshifted at ≈2.0817 μm and can be observed with the HK filter (central λ = 1.93 μm).

Unfortunately, due to poor weather conditions, we did not achieve the requested sensitivities. No BLR emission lines have been detected in any of the quasars except J2310+1855, from which we derived a CIV-based M_BH estimate of 6 × 10⁹ M_⊙ (see Appendix D for full details of the observations). Our estimation is consistent, within the error, with that reported by Feruglio et al. (2018) and Shen et al. (2019). We note that Shen et al. (2019), who published NIR spectra of a large sample of z ∼ 5.7 QSOs, also provide a M_BH measurement for J2310+1855 through virial relation based on MgII as well. In the following, we refer to the MgII estimate from Shen et al. (2019), because of the aforementioned issues related to CIV-based measurements.

6. Summary of sample sub-selections

The data analysis described in Sects. 3 and 4 was performed on the 32 continuum-subtracted cubes of the sources listed in Table 1. Each step of the analysis has led to the rejection of a number of objects that turned out not to be suitable for the method adopted in this work. Here, we briefly summarise the different sub-selections used throughout this work:

1. By inspecting the velocity maps and red/blue residuals maps (see Sect. 3.3 for details), we found that 14 out of 32 objects (∼45%, flagged with [rot] = “u” in Table 1) do not show spatially resolved kinematics or rotating disc kinematics. This is possibly due to the presence of outflows or merging events, or of a companion located in proximity (projected on the sky plane) of the quasar. For this purpose, velocity-dispersion maps provide additional information on the kinematics. However, a comprehensive interpretation of the complex velocity fields observed in these sources is beyond the scope of this work. As a result of this analysis, the sample has been reduced from 32 to 18 objects. On the other hand, excluding these objects from the final sample may introduce bias in the final results (see Sect. 7.4).

2. We then performed the fits of the flux and velocity maps (see Sects. 4.2 and 4.3 for details) on the remaining 18objects selected in the previous step. As a result, for five objects (∼30%), the kinematical modelling has not enabled us to constrain the disc inclination, and consequently the dynamical mass. It is possible that incorrect assumptions on the mass distribution (see Eq. (1)) and/or the poor angular resolution of the observations making inclination and dynamical mass almost degenerate parameters (see Sect. 4.3) have prevented the determination of the mass in these host galaxies. In addition, the iso-velocity curves of the kinematical fields are typically distorted, suggesting the presence of non-circular motion. In particular, for two objects (PSO J308−21, VIKING J0305−3150), the posterior probability distributions of β and μ are bimodal, preventing us from constraining these parameters. For three objects (PSO J183+05, SDSS J0129−0035, ULAS J1319+0950), we derived an upper limit on the disc inclination, which is a lower limit on mass M_dyn. In summary, we obtained 8 constrained, 2 bimodal, and 3 lower limit measurements of the dynamical mass (see Table 2).

3. The final step is to relate our dynamical mass estimates with M_BH retrieved from the literature. We illustrate this step in Sects. 5 and 8. Despite several studies performed in this field aiming to estimate M_BH even for high-redshift quasars, for two objects (∼15% of the remaining 13 resulted from the previous step), black hole mass estimates were not available at the time this paper was written. Therefore, we rejected these objects from the final sample (see Table 3). These quasars are flagged with [M_BH] = “u” in Table 1.

These selection steps are shown in the scheme drawn in Fig. 8. Overall, we were able to obtain a measurement of galaxy dynamical mass and retrieve black hole mass only for eight sources, that is ∼10% of the initial sample of 72 QSOs.

Fig. 8. Summary of the different sub-selections used throughout this work. Starting from an initial sample of 72 quasar host galaxy observations extracted from the ALMA data archive, the final sample is composed of ten high-z objects for which we study the BH-galaxy relation. The characters indicating the selection type are in accordance with Table 1. The number of objects rejected from the sample is also indicated at each step.

7. Comparison of results and uncertainties on dynamical mass estimates

In Sect. 7.1, we compare our results with those obtained by other authors who attempted to perform a full kinematical modelling of individual sources that belong to our sample. Other works highlighted the presence of companion sources in the close environment of a few QSOs analysed in this work. Such satellite galaxies can disturb the gas kinematics of the host through tidal interaction, thus introducing additional uncertainties in measuring the galaxy dynamical mass. We discuss this point in Sect. 7.2. In Sect. 7.3, we discuss the limit of validity of the assumption of a rotating disc model and the possible mass contribution arising from random motions throughout the galaxy (see also Appendix B). Finally, in Sect. 7.4, we investigate observational biases possibly arising from the sub-selection of the sample.

7.1. Comparison of results from other kinematical modelling in the literature

In the following sections, we discuss the results obtained from the kinematical modelling of J1319+0950, J0305−3150, J1044−0125, and J0129−0035, making a comparison between our results and those obtained in previous works.

7.1.1. ULAS J1319+0950

Jones et al. (2017) and Shao et al. (2017) performed a kinematical characterisation of [CII] emission of J1319+0950 by using a tilted rings warpless model and assuming purely circular rotation. They inferred, respectively, a dynamical mass of 15.8 × 10¹⁰ M_⊙ and 13.4 × 10¹⁰ M_⊙ and an inclination angle of 29 deg and 33 deg (which are roughly consistent with the value estimated by Wang et al. 2013) by using the axial ratio of the [CII] flux map to estimate the disc inclination and the [CII] FWHM as estimate of maximum circular velocity. Furthermore, by fitting the dust continuum emission in UV plane, Carniani et al. (2019) inferred an inclination angle of ∼15 deg. In contrast, we are not able to constrain the disc inclination angle from our kinematical modelling, resulting in a lower limit on dynamical mass (> 4 × 10¹² M_⊙). The disagreement between our result and the previous ones may arise from the beam smearing effect that is not taken into account in the model used by Jones et al. (2017) and Shao et al. (2017). As discussed in Sect. 4.3, beam smearing strongly affects the observed velocity gradients and introduces additional uncertainties in the fitting parameters. In addition, Shao et al. (2017) mentioned that current data cannot fully rule out the presence of a bidirectional outflow, which introduces additional uncertainty regarding the dynamical mass. In such a case, the strong deviation of the ratio M_BH/M_dyn could be associated with an incorrect kinematical modelling of the observed data, for which we also assume rotating disc kinematics, like in Jones et al. (2017) and Shao et al. (2017).

7.1.2. VIKING J0305−3150

High angular resolution (0.076″ × 0.071″) ALMA [CII] observations of J0305−3150 were recently presented and analysed by Venemans et al. (2019). The resulting analysis highlights that the distribution and kinematics, as traced by the [CII] emission, are highly complex and include the presence of cavities and blobs.

Venemans et al. (2019) attempted to model the kinematics adopting different 3D models (thin rotating disc with constant velocity, Keplerian disc, truncated disc, and a simple AGN model embedded in a uniform rotating spherical gas) taking into account beam smearing effects and pixel correlation within the beam with a Bayesian approach. The results show that the gas kinematics in J0305−3150 appear to be dispersion-dominated, with some overall rotation in the central kiloparsecs, and cannot easily reproduced by a simple rotating disc model with the implication that most of the gas has not settled in a disc yet. In addition, authors point out that energy injection into the ISM produced by AGN feedback processes, and the presence of a companion in the close environment, may play an important role in producing the observed [CII] cavities and in perturbing the gas kinematics. In conclusion, a simple model of a rotating disc is not sufficient to match the [CII] observations of J0305−3150 also derived by our analysis, where M_dyn is unconstrained by the simple model assumed.

7.1.3. SDSS J1044−0125 and SDSS J0129−0035

In the work by Wang et al. (2019), authors carried out observations of J1044−0125 and J0129−0035 through the ALMA program 2012.1.00240.S (the same dataset used in this work for the latter source) with angular resolution of ∼0.2″. The authors show that gas in J1044−0125, as traced by [CII] emission, does not show a clear sign of rotation, suggesting a very turbulent gas velocity field. Furthermore, the [CII] spectrum reveals offset components that could be associated with a node of outflowing gas or the dense core of a satellite galaxy, which contribute to increasing the velocity-dispersion component of the gas in the host galaxy. On the other hand, the lower angular resolution data used in our work (∼0.6″ × 0.5″, ALMA programme 2011.0.00206.S, see Wang et al. 2013) reveal the presence of a velocity gradient. This could be the result of beam smearing effects producing a smoothing of the rapidly changing velocity gradients. In the case of J1044−0125, we find that the observed velocity field is roughly consistent with a nearly edge-on rotating disc model. Therefore, we conclude that our dynamical mass estimate is tentative. We also note that Wang et al. (2019) show that [CII] and dust emissions in the nuclear region of J1044−0125 and J0129−0035 follow an exponential light profile, in accordance with the hypothesis at the base of our model.

In the case of J0129−0035, the observations analysed in Wang et al. (2019) reveal that [CII]-emitting gas shows clear velocity gradients likely associated with a rotating disc with additional gas clumps, thus suggesting complex kinematics in the nuclear region. They attempted to constrain the host galaxy dynamics adopting the same method as in the works of Jones et al. (2017) and Shao et al. (2017). The results show that the kinematics are consistent with a nearly face-on rotating disc with an inclination angle of β = (16 ± 20) deg and a lower limit on the dynamical mass equal to M_dyn = 2.6 × 10¹⁰ M_⊙. The results are consistent with what we found in this work. The BH mass of J0129−0035 is estimated as in Wang et al. (2013), assuming Eddington accretion, and is the same one that we used in this work. Hence, Wang et al. (2019) estimated an SMBH to host a galaxy dynamical mass ratio of M_BH/M_dyn = 0.0066 to be compared with M_BH/M_dyn = 0.0022, which is roughly consistent with the local ratio predicted in Decarli et al. (2010), unlike the most luminous quasars with massive BHs (10⁹ − 10¹⁰ M_⊙) at this redshift that show ratios from a few to ≳10 times higher (Venemans et al. 2016; Decarli et al. 2018). Therefore, as pointed out by Wang et al. (2019), this result may suggest that the BH-galaxy coevolution of a less massive system (M_BH ∼ 10⁷ − 10⁸ M_⊙) in the early Universe is closer to the trend of local galaxies (see also, Willott et al. 2010, 2015a, 2017; Izumi et al. 2018, 2019).

7.2. Possible contamination due to the presence of companion sources in the quasar’s local environment

Decarli et al. (2017) serendipitously discovered companion galaxies in the ALMA field of SDSS J0842+1218, CFHQS J2100−1715, PSO J231−20, and PSO J308−21. Such companions appear similar to the host galaxies of quasars in terms of [CII] brightness and implied dynamical mass, but do not show evidence of AGN activity. In our work, we analysed the same dataset as Decarli et al. (2017) (ALMA programme 2015.1.01115.S), concluding that the kinematics are unresolved (flag [rot] = “u”; see Table 1) in the case of J0842+1218 and J2100−1715 (beam size of ∼1.0″ × 0.9″ and ∼0.7″ × 0.6″, respectively); and marginally resolved (disc inclination angle is unconstrained; flag [β] = “u”; see Table 1) in the case of J231−20 (beam size of ∼1.0″ × 0.8″). This last source together with J308−21 has a [CII]-bright companion at small projected separation of ∼10 kpc, suggesting a strong gravitational interaction between quasar and companion able to alter the disc kinematics increasing the velocity-dispersion component of the gas. In particular, Decarli et al. (2017) show that the [CII] emission of J308−21 stretches over about 4″ (≈25 kpc) and more than 1000 km s⁻¹ connecting the companion source suggesting that is undergoing a tidal disruption due to the interaction or merger with the quasar host. This scenario is successively supported by higher angular resolution (∼0.3″; ALMA programme 2016.A.00018.S) follow-up observation of J308−21 presented in Decarli et al. (2019); the same dataset analysed in this work. However, the bulk of [CII] emission of the quasar host galaxy shows a spatially resolved velocity gradient, which, in our work we attempt to model with a rotating disc by excluding pixels that are clearly not associated with the quasar host. Nevertheless, our analysis leads to a bimodal posterior probability distribution of disc inclination angle and dynamical mass parameters of J308−21. We can thus conclude that the complex kinematics of this system highlighted in the previous analysis presented in Decarli et al. (2017, 2019), cannot be easily interpreted with a simple rotating disc, likely due to the perturbed kinematics caused by the strong interaction with the satellite galaxy.

Willott et al. (2017) analysed the source PSO J167−13 observed in ALMA Cycle 3 project 2015.1.00606.S; the same dataset analysed in this work. This source shows an asymmetric continuum emission that is more extended to the south–east than north–west of the peak. This excess is located at ≈0.9″ (projected distance ≈5.0 kpc), and it is associated with a companion galaxy whose [CII] blueshifted (270 km s⁻¹) emission corresponds to about 20% of the QSO [CII] luminosity. The P–V diagram of the source shows a positive velocity gradient, suggesting a rotating disc geometry. With this assumption, Willott et al. (2017) infer the dynamical mass of M_dyn = 2.3 × 10¹¹ M_⊙ using the axial ratio of the quasar (excluding the companion source) [CII] flux map as an estimate of the disc inclination angle. The black hole mass of J167−13, M_BH = (4.0 ± 2.0)×10⁸ M_⊙, is estimated in Venemans et al. (2015) using calibration based on MgII broad emission line (Vestergaard & Osmer 2009). By comparison, we measure a dynamical mass of and a black hole mass of M_BH = (2.7 ± 1.4)×10⁸ M_⊙ (see, Sect. 5 for details) resulting in a ratio of M_BH/M_dyn = 0.0016, which is completely consistent with the value found by Willott et al. (2017) (M_BH/M_dyn = 0.0017) and with the prediction of the local relation (Decarli et al. 2010; Kormendy & Ho 2013).

Neeleman et al. (2019) further investigated the aforementioned four quasar host-companion galaxy pairs of J0842+1218, J2100−1715, J231−20, J167−13 by analysing high angular resolution (∼0.4″ × 0.3″) ALMA observations of [CII] emission. They observe tidal interactions disturbing the gas in these high-z galaxies making the ISM turbulent and thus confirming previous results of Decarli et al. (2018, 2019) and Willott et al. (2017). Furthermore, these high angular resolution observations reveal that [CII] emission of SDSS J1306+0356 arises from two spatially and spectrally distinct sources with a physical separation of 5.4 kpc that are interacting tidally. Neeleman et al. (2019) modelled the [CII] kinematics of the galaxy pairs with a rotating thin disc model, taking into account the beam smearing and the correlation between pixels. They obtained upper limits on dynamical masses for all the sources except J167−13 and J2100−1715. In particular, they measured a dynamical mass of (3.5 ± 0.4)×10¹⁰ M_⊙ for the J167−13 quasar. This value is about one order of magnitude lower than the result of Willott et al. (2017) and our work. In our estimate, we also take into account the emission from the companion galaxy, thus possibly overestimating the quantities derived from the total integrated spectrum (FWHM_[CII], [CII] flux, luminosity, etc., see Table A.1), the scale radius of the mass profile and the total mass content of the quasar host galaxy. This could explain the inconsistency in our dynamical mass measurements of J167−13 compared to the works of Willott et al. (2017) and Neeleman et al. (2019).

7.3. Limit on the assumption of thin rotating disc

The comparison of our results discussed in Sects. 7.1 and 7.2 highlights that, at least in some cases, the assumption of a thin rotating disc is too simplified to properly describe the observed complex kinematic field. Furthermore, from the analysis of the velocity maps, we find extreme cases in which the disc inclination with respect to the sky plane is very low (e.g. ULAS J1319+0950), compatible with a face-on configuration. However, for these objects, the observed velocity dispersions are still high (∼100−200 km s⁻¹, as is clear from the figures in Appendix E), which is not expected for thin face-on discs.

The observed velocity dispersion can be produced both by instrumental effect and random motions throughout the galaxy (see e.g. Flores et al. 2006; Weiner et al. 2006; Epinat et al. 2010). Different authors (e.g. Cresci et al. 2009; Epinat et al. 2009; Taylor et al. 2010) pointed out that random motions can support part of the mass. In this case, modelling the kinematics with a rotating disc could result in underestimating the galaxy’s dynamical mass. In Appendix B, we investigate the contribution of random motions to the dynamical mass and find that, in our sample, the mass supported by non-rotational motions is negligible, meaning it is included within the dynamical mass uncertainties. Therefore, we conclude that the rotating gas disc model provides an overall good description for the gas kinematics of our QSO host galaxies.

7.4. Potential observational biases in excluding the unresolved objects

In Sect. 3.3, we investigated whether the observed kinematics are spatially resolved. Out to 32 objects with a line detection, 14 (≳40%) were rejected from the final sample (see also Sect. 6). Excluding those objects that are spatially unresolved could result in an observational bias. In fact, if these sources were significantly less massive than the others, the final results might therefore be biased towards more massive host galaxies.

The observed size of the line emitting region may depend on both the achieved sensitivity and the angular resolution. Therefore, in the case of spatially unresolved emission, it is difficult to asses if this is due to the intrinsic compact size of the galaxy or to the low sensitivity level of the observations. For this purpose, deeper observations of these objects with similar observational setups could help us to make a fair comparison of the observed size. However, we do not observe a strong correlation between the spatial size of the FIR line emission and the dynamical mass of the galaxy (see Table 2). Therefore, we conclude that we cannot safely argue that a possible bias is introduced in rejecting the spatially unresolved objects.

8. The M_BH-M_dyn relation at high redshift

In order to trace the relation between black hole mass and dynamical mass for the final high-z QSOs sample, we compared the M_dyn measurements obtained through kinematical modelling illustrated in Sect. 4 (see Table 2), with the black hole masses obtained from literature as we explained in Sect. 5 (see Table 3). The relation is shown in the plot of Fig. 9. We also report two reference relations obtained with samples of local quiescent galaxies (Kormendy & Ho 2013; de Nicola et al. 2019) and AGN (Reines & Volonteri 2015). In order to infer the average redshift evolution of the M_BH−M_dyn, we adopted the relation log M_BH = α + β(log M_dyn − 10.8), and we performed the fit assuming fixed slope β = 1.01 ± 0.07 as found by de Nicola et al. (2019), and the normalisation α as the only free parameter. Furthermore, to reduce the impact of any possible outliers, we executed the fit adopting the bootstrap method on the standard χ² minimisation.

Fig. 9. High-redshift relation between the black hole mass (MBH) and the dynamical mass of the host galaxy (Mdyn). The dashed black line and the dotted red line represent the reference local relation inferred using samples of local galaxies (E = ellipticals, S/S0/Sb = spirals) shown as black dots (Kormendy & Ho 2013, also indicated as KH+13) and red triangles (de Nicola et al. 2019, or DN+19). The green line is the relation found by Reines & Volonteri (2015, or RV+15) by measuring the total stellar mass in a sample of the local AGN (green crosses). The solid blue line is the best fit to our data. The shadowed areas show the 1σ uncertainty. In the case of SDSS J0129−0035 and ULAS J1319+0950, we inferred a lower limit on the dynamical mass. We do not take these data into account in the fit. The circles of our data points indicate the sources for which the BH masses are estimated from bolometric luminosity assuming Eddington accretion.

Using 10⁴ bootstrap iterations, we obtained the best value of α and its uncertainties by computing the 16th, 50th, and 84th percentiles, respectively:

(6)

Our result is in agreement with those reported by other high-z works (e.g. Decarli et al. 2010, 2018; Trakhtenbrot et al. 2015, 2017; Venemans et al. 2016, 2017a) suggesting that the M_BH−M_dyn relation evolves with redshift. It should be noted that the local reference relation (e.g. Kormendy & Ho 2013; de Nicola et al. 2019) is obtained using bulge stellar mass in spiral and elliptical galaxies (where, in the latter case, bulge stellar mass corresponds to the total stellar mass). As a result, the galaxy dynamical masses estimated in this work should be treated as an upper limit of the total stellar mass. By comparing our results with the relation by Reines & Volonteri (2015), who adopted the total stellar mass of the AGN host galaxy (green line in Fig. 9), we find an even stronger evolution with redshift.

The evolution of the ratio Γ = M_BH/M_dyn as a function of redshift provides key information about the relative time scale between black hole growth and galaxy mass assembly. For this purpose, we show the M_BH/M_dyn ratios as a function of z estimates obtained from the integrated spectra of the lines provided in Table 2. The final result is shown in Fig. 10, where we also overplot the relation found by Decarli et al. (2010) using galaxy stellar masses of a sample of quasars at z ≲ 3, extrapolated up to z = 7:

(7)

Fig. 10. Evolution of Γ = MBH/Mdyn as a function of redshift z. The black dotted and dash-dotted lines represent, respectively, the relation found by Decarli et al. (2010) at z ≲ 3, and the corresponding ratio at z = 0. The inset panel shows the same plot at 2 <  z <  3. In the case of SDSS J0129−0035 and ULAS J1319+0950 we inferred lower limits on dynamical masses (i.e. an upper limit on the ratio Γ). The circles indicate those sources which BH masses are estimated from bolometric luminosity assuming Eddington accretion.

We conclude that the trend of Γ that we inferred at high redshift is roughly consistent with that of Eq. (7), and therefore this result confirms the evidence that M_BH−M_dyn appears to evolve with the redshift, as has been highlighted in previous works (Walter et al. 2004; Decarli et al. 2010, 2018; Merloni et al. 2010; Venemans et al. 2012, 2017a). The ratio Γ appears to be ∼10× the local value logΓ(z = 0) = 0.28 at z ∼ 4−6. However, from Fig. 10, we can infer that for SDSS J0923+0247, SDSS J0129−0035, ULAS J1319+0950, and PSO J167−13 at redshift 4.6 ≲ z ≲ 6.6, the M_BH/M_dyn is consistent with the value observed in galaxies in the local Universe. Although a very preliminary result, this possibly suggests a decreasing of Γ at z ≳ 6.

The discussions reported in Sects. 7.1 and 7.2 point out that at least some of the galaxy masses estimated in this work could suffer from large uncertainties associated with the simple assumptions that are the basis of the fitting method. Therefore, although M_BH estimates are strongly affected by the large (systematic) uncertainties associated with M_BH measurements at high-z (up to ∼0.4 dex, see Sect. 5), we conclude that the observed M_BH/M_dyn values are also likely affected by uncertainties in M_dyn.

Our method did, however, allow us to obtain accurate galaxy dynamical mass estimates at such high redshift. We find that the spread in M_BH/M_dyn values at z ∼ 4−7 is much greater compared with that of local galaxies. This suggests that the observed spread could not arise from the large uncertainties associated with rough galaxy (virial) mass estimates at high-z usually adopted, but it could have a physical reason.

9. Discussion

In order to extend the context of our work, we compare our results with both observational and theoretical predictions of the BH-galaxy relation obtained in other works. Then, we investigate observational biases possibly affecting the results, and we test the reliability of galaxy virial mass estimates.

In Sect. 9.1, we compare the M_BH−M_dyn and Γ − z relation presented in this work with other results on the high-z BH-galaxy relation and we discuss the effect of observational biases. In Sect. 9.2, we compare the BH-galaxy relation prediction from recent simulations of galaxy evolution, and we discuss potential issues on dynamical mass estimates at high redshift outlined from theoretical models. Finally, in Sect. 9.3, we test the reliability of virial mass estimates by comparing them with the dynamical mass measurements.

9.1. Observational biases: comparison with other results on the early BH-galaxy relation

Previous works that aimed to study the BH-host galaxy co-evolution in the early epochs show that z  ≳  6 luminous quasars have BH-to-host-galaxy-mass ratios ∼10 times larger than the typical value observed in the local Universe, implying that these SMBHs formed significantly earlier than their hosts (e.g. Walter et al. 2004; Maiolino 2009; Merloni et al. 2010; Decarli et al. 2010, 2018; Venemans et al. 2012, 2016; Wang et al. 2013, 2016; Trakhtenbrot et al. 2015). However, these results may be affected by observational biases.

Because most luminous quasars are powered by the most massive BHs at high redshift, if there is a scatter in BH-host galaxy mass relation, for a given M_BH, the selection of objects with low galaxy mass is favoured due to the steepness of the galaxy mass function at its high-mass end (see e.g. Grazian et al. 2015; Song et al. 2016), thus producing an artificially high average M_BH/M_dyn (Lauer et al. 2007; Schulze & Wisotzki 2014). In order to investigate this selection bias effect, in Fig. 11, we compare the distribution of our dynamical mass estimates, with the galaxy stellar mass function at different redshift. For this purpose, our dynamical mass measurements represent upper limits on galaxy stellar masses (M_⋆). Most of the quasars are at the knee of the quasar luminosity function (Song et al. 2016), indicating that they represent the bulk of the quasar population at such high redshifts. On the other hand, J183+05, J167−13, and the most extreme J1319+0950 are at the massive end of the M_⋆-function, and they may be more affected by the “Lauer” bias. Notwithstanding all considerations of the reliability of dynamical mass estimates of the aforementioned objects, we can conclude that these three quasars may be considered the most evolved system known at z ∼ 6 in terms of galaxy mass.

Fig. 11. Comparison between our dynamical mass measurement distribution and galaxy stellar mass function. Upper panel: Galaxy stellar mass function at different redshifts (Song et al. 2016). The dashed vertical lines show the position of the M⋆ values that represent the “knee” of the mass function. Lower panel: stacked histogram of our dynamical mass estimates.

Interestingly, J167−13 together with J0923+0247 and J0129−0035 in our sample have a black hole mass M_BH <  5.0 × 10⁸ M_⊙ (see Table 3), and all of them are consistent with the local BH-galaxy relation at z ∼ 0. As discussed in Venemans et al. (2016), Wang et al. (2016), and Willott et al. (2015a, 2017), this may suggest that, while the z ∼ 6 quasars with BH masses to the order of 10⁸ M_⊙ are close to the relation valid for their local counterpart, the most massive BHs (M_BH >  10⁹ M_⊙) at the earliest epochs grow faster than the quasar host galaxy and tend to be above the trend of local galaxies. This may imply that actually, there is no strong correlation between the two properties in high-redshift quasars, but that the scatter was much larger in the early Universe than today. However, to confirm this conclusion, a wide range of BH masses and larger sample are required to overcome the observational bias due to the intrinsic scatter of the BH-galaxy relation. In this context, Izumi et al. (2018, 2019) studied the M_BH−M_dyn using a sample of seven z ≳ 6 low-luminosity quasars (absolute magnitude at 1450 Å, M₁₄₅₀ >  −25) targeted in [CII]_158 μm by ALMA. They derived the quasar host galaxy dynamical (virial) masses assuming rotating disc geometry, and the axial ratio of [CII] flux map as a proxy of the disc inclination angle. Furthermore, they estimated M_BH through a SE virial relation or assuming Eddington-limited accretion. Izumi et al. (2018, 2019)’s results show that while the luminous quasars (M₁₄₅₀ <  −25) typically lie above the local relation (Kormendy & Ho 2013) with BHs overmassive compared to local AGNs, the discrepancy becomes less evident at M_dyn ≳ 10¹¹ M_⊙ (see also e.g. Trakhtenbrot et al. 2017). On the other hand, most of the low-luminosity quasars show comparable or even lower ratios than the local one, particularly at a range of M_dyn ≳ 4 × 10¹⁰ M_⊙. Therefore, Izumi et al. (2018, 2019) show that, at least in this high M_dyn range, previous works based on sample of luminous quasars might have been biased toward more massive black holes. If these results were confirmed with future follow-up observations, the observed evolution of M_BH/M_dyn ratio out to z ∼ 2−3 could be explained as the result of sample selection bias only. However, Izumi et al. (2018, 2019) could be biased, due to the use of virial masses. As we illustrate in Sect. 9.3, the galaxy masses estimated in this work through a full kinematical model, are not correlated with those estimated through virial theorem. This fact suggests that the use of galaxy virial mass in studying the BH-galaxy relation, could be reflected in an increasing scatter in the observed relation.

9.2. Comparison with recent theoretical models and simulations

The benchmark correlation in local galaxies is based on bulge stellar masses (Kormendy & Ho 2013). In high-redshift quasars, the host galaxies are completely outshone by the central emission, and the limited angular resolution of the current UV-based observations does not allow us to easily decouple the quasar from its host. Therefore, estimating the stellar mass content in a galaxy out to z ∼ 2−3 is very challenging (see Sect. 1). However, at such high redshifts, the galaxies’ bulges may not have formed yet, or they cannot be detected. Therefore, in the high redshift studies, the dynamical mass of the host galaxy is estimated from the gas properties in the sub-mm observations, and it is usually used as a proxy of bulge stellar mass.

Beyond the observational biases affecting these studies (see Sect. 9.1), Lupi et al. (2019) recently pointed out that the high-z deviation from the local M_BH−M_dyn relation might be due to the different tracers used to estimate the mass of galaxies at high redshift, namely gas-based dynamical mass and the stellar mass. Lupi et al. (2019) performed a high-resolution cosmological zoom-in simulation in order to investigate the evolution of quasar hosts by properly resolving both the distribution of the cold gas phase (30 K <  T <  3000 K) traced by [CII]_158 μm line emission, and stars traced by far-UV flux. Their results show that the gas settles in a well-defined dense thin disc extending out to ∼1 kpc already at z  ∼  7. Adopting the techniques used in observational studies, they derived dynamical mass estimates through the virial theorem under the assumption of rotationally supported systems. Comparing these dynamical (virial) mass estimates with the masses of the central BHs, they obtained an average BH-dynamical mass ratio of M_BH/M_dyn ∼ 0.017, which is ∼15 times greater than local values (Decarli et al. 2010; Kormendy & Ho 2013; Reines & Volonteri 2015), in agreement with previous high-z observations and simulations (see e.g. Walter et al. 2004; Venemans et al. 2017b; Barai et al. 2018; Decarli et al. 2018), and roughly consistent with our results. Additionally, they also compared galaxy stellar masses with BH masses, showing no sign of clear deviation with the local relation (Reines & Volonteri 2015). This result implies that dynamical mass estimated using the virial theorem underestimates the dynamical mass of the system. A similar discussion was recently reported by Kohandel et al. (2019), who analysed the kinematical properties of a simulated star forming galaxy at z ∼ 7. They show that using the virial theorem in a rotationally supported system, the dynamical mass estimates suffer from large uncertainties depending on the disc inclination.

The approach proposed in our work, in which we perform a kinematical modelling of observed velocity fields, allows us to infer both disc inclination and dynamical mass, thus reducing the uncertainties and biases on our estimates. However, larger samples with higher angular resolution observations are required to finally assess whether a deviation from the local relation exists or not.

9.3. Comparison between virial masses and dynamical mass estimates

We tested the reliability of virial mass estimates by comparing them with our dynamical mass measurements. For this purpose, we made rough dynamical (virial) mass (M_vir) measurements of our galaxy sample following, e.g. Wang et al. (2013), Willott et al. (2015a), Decarli et al. (2018):

(8)

where R_em is the radius of the emitting region, and FWHM_line is the full width at half maximum of the line emission.

We performed 2D Gaussian fits, within CASA, of the flux maps obtained in Sect. 3, and we estimated the deconvolved major (a_maj) and minor axes (b_min) of the best model. Then, we computed R_em as a_maj/2 in physical length using the redshift estimates obtained in Sect. 3.1 and the disc inclination angle as sin²β = 1 − (b_min/a_maj)². Finally, we retrieved the FWHM_line from the Gaussian fit of the line spectra (see Table A.1). The results of the 2D Gaussian fits, the R_em values, and virial masses are listed in Table C.1. In Fig. 12, we compare dynamical virial masses estimated through Eq. (8) with the dynamical mass measurements obtained in this work through a full kinematical modelling (see Table 2).

Fig. 12. Comparison between dynamical virial mass (Mvir) computed using flux map properties and dynamical mass (Mdyn) obtained with the full kinematical model of the gas velocity field. The dashed line indicates the 1:1 relation, while dotted lines are 1 dex shifts. The objects represented with star symbols are upper limits on virial mass, since their emission appears consistent with a point-like source (as a result of 2D Gaussian fits within CASA). In these cases, following Willott et al. (2015a), we assume β = 55 deg as disc inclination angle. The black square indicates PSO J183−05, for which we retrieved a lower limit on Mdyn, but for which MBH is not available.

We conclude that M_vir and M_dyn are roughly in good agreement, but they appear not to be correlated, confirming that virial mass is not a reliable dynamical mass estimate of the host galaxy. We note that the errors on virial mass measurements are statistical errors, ignoring any intrinsic uncertainties and systematic biases associated to the virial assumption.

In order improve the galaxy virial mass estimate, we can use the spectroastrometry method by Gnerucci et al. (2011). With this method, is possible to probe spatial scales smaller than the angular resolution, thus allowing a more accurate measurement of the dimension of the line-emitting region. On the other hand, the mass estimates are affected by uncertainties associated with the measurement of the galaxy disc inclination angle, for which it is possible to use the axial ratio from galaxy morphology. In Appendix C, we compare the mass factor M_dynsin²β obtained through the full kinematical model, the virial formula, and the spectroastrometry method. The results show that spectroastrometry is a robust proxy for the galaxy dynamical mass in contrast with the “classical” virial estimates usually adopted.

10. Conclusions

In this work, we investigated the relation between supermassive black hole mass (M_BH) and the dynamical mass of the host galaxy (M_dyn) of a sample of 10 quasars at 2.3 ≲ z ≲ 6.5 targeted in either [CII]_158 μm or CO rotational transitions by ALMA. We then studied the evolution of Γ = M_BH/M_gal across cosmic time.

Previous works exploiting ALMA observations attempted to trace the M_BH−M_dyn relation at high-z by estimating the galaxy mass through virial theorem and thus possibly introducing significant uncertainties and biases. To avoid such large uncertainties, we performed a kinematical modelling of the cold gas in the hosts taking into account the beam smearing effect.

In summary, we conclude that:

• The galaxy mass estimated using the virial theorem combining the axial ratio of flux map to estimate the disc inclination angle, and line FWHM as a proxy of circular velocity, suffers from large uncertainties and could underestimate the dynamical mass of the system (Lupi et al. 2019; Kohandel et al. 2019).

• The beam smearing effect strongly affects the observed velocity field in the host galaxy, making the disc inclination angle and galaxy dynamical mass almost degenerate parameters. The more the angular resolution decreases, the more significant this effect becomes, and it should be taken into account in kinematical modelling.

• The dynamical masses estimated from the kinematical modelling highlight evidence of the evolution of the M_BH−M_dyn relation consistently with previous works (e.g. Walter et al. 2004; McLure et al. 2006; Maiolino 2009; Bennert et al. 2010, 2011; Decarli et al. 2010, 2018; Merloni et al. 2010; Wang et al. 2010, 2013, 2016; Canalizo et al. 2012; Targett et al. 2012; Venemans et al. 2012, 2016, 2017a; Bongiorno et al. 2014; Trakhtenbrot et al. 2015). In particular, we conclude that, on average, our sample is placed above the reference relation found for galaxies in the local Universe. The normalisation α of the BH-galaxy relation is such that, on average, for a given value of M_dyn, M_BH is ∼10× higher compared to that found by de Nicola et al. (2019), and ∼150× higher than what Reines & Volonteri (2015) found using the total stellar mass of local AGNs.

• The ratio Γ = M_BH/M_dyn appears to be ∼10× the local value at z ∼ 4−6, consistent with the result found by Decarli et al. (2010) extrapolated up to z = 7, except for four objects at z ∼ 4−6 that show Γ ratios consistent with the local one (Γ(z = 0)). Despite the low statistics, this is the first evidence of a Γ value decreasing at z ∼ 6. We are possibly witnessing the phase in which a black hole rapidly grows with respect to the galaxy mass.

• The observed spread in M_BH/M_dyn values at z ∼ 4−6 is much greater compared to galaxies in the local Universe. Given the accurate galaxy dynamical mass estimates obtained in this work, the observed spread could be due to physical factors, and not associated with the large uncertainties affecting the galaxy virial mass estimates usually adopted in high-z studies.

• The sources in our sample with M_BH to the order of 10⁸ M_⊙ are close to the relation found for galaxies in the local Universe (Kormendy & Ho 2013; de Nicola et al. 2019), while the most massive BHs (M_BH >  10⁹ M_⊙) lie above them, thus suggesting a faster evolution with respect to their host at z ∼ 6.

• Most of our sample represents the bulk of the quasar population at z >  4; thus, overall, the selection of our galaxy sample is not strongly affected by the “Lauer” bias (Lauer et al. 2007; Schulze & Wisotzki 2014). However, a wide range of BH masses and a larger sample is required in order to avoid the observational bias resulting from the intrinsic scatter in the M_BH−M_dyn relation.

Based on our blind search, we conclude that one third of high-z quasar hosts have gas kinematics consistent with rotating discs, but it is still very challenging to infer the dynamical mass due to the poor angular resolution and sensitivity of current observations. The typical angular resolution of the observations (∼0.5″) is frequently not good enough to constrain the dynamical parameters of the discs at z ≳ 5, and the fitting procedures cannot take into account possible distortions of the velocity field introduced by instrumental effects. As a result, we inferred the dynamical masses only for ten out of 72 quasars observed with ALMA so far.

On the other hand, for those quasars with deep ALMA observations and high angular resolution, this work shows that dynamical mass estimations are also feasible at z ∼ 6. Further ALMA high angular resolution observations of high-z quasars are crucial to studying the evolution of the M_BH/M_dyn ratio and verifying whether Γ(z) decreases at z ≳ 6 as suggested by our preliminary results.

Acknowledgments

We thank the anonymous referee for her/his careful reading of the manuscript and her/his comments which really helped us to improve the paper. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00206.S, ADS/JAO.ALMA#2012.1.00240.S, ADS/JAO.ALMA#2012.1.00882.S, ADS/JAO.ALMA#2013.1.00815.S, ADS/JAO.ALMA#2013.1.01153.S, ADS/JAO.ALMA#2015.1.00228.S, ADS/JAO.ALMA#2015.1.00399.S, ADS/JAO.ALMA#2015.1.00606.S, ADS/JAO.ALMA#2015.1.01115.S, ADS/JAO.ALMA#2015.1.01247.S, ADS/JAO.ALMA#2016.1.00544.S, ADS/JAO.ALMA#2016.A.00018.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. We acknowledge the support from the LBT-Italian Coordination Facility for the execution of observations, data distribution and reduction. The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporation partners are: The University of Arizona on behalf of the Arizona university system; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, the Astrophysical Institute Potsdam, and Heidelberg University; The Ohio State University, and The Research Corporation, on behalf of The University of Notre Dame, University of Minnesota, and University of Virginia. SC acknowledges support by the European Research Council No. 740120 “INTERSTELLAR”. MP is supported by the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. RM acknowledges supports by the Science and Technology Facilities Council (STFC) and from ERC Advanced Grant 695671 “QUENCH”.

References

Appendix A: Integrated spectra: line properties and derived quantities

From the best-fit of the integrated spectra (see Sec. 3.1), we directly retrieved the line FWHM, and the velocity-integrated flux of the line (F_line). Then, we also inferred line luminosity (L_line), [CII] mass (M_[CII]), total gas mass (M_gas), and the [CII]-based star formation rate (). The line luminosities were computed following Solomon & Vanden Bout (2005):

(A.1)

where F_line is in unit of Jy km s⁻¹; ν_rest in GHz, and D_L in Mpc.

Then, by analogy with Venemans et al. (2017b), assuming optically thin [CII] emission and local thermodynamical equilibrium (LTE) of the carbon line, we estimated the mass of singly ionised carbon in galaxies as:

(A.2)

where C is the conversion factor between pc² and cm², m_C the mass of a carbon atom, A_ul = 2.29 × 10⁻⁶ s⁻¹ the Einstein coefficient, Q(T_ex) = 2 + 4e^−91.2/T_ex the CII partition function, and T_ex the excitation temperature that we set equal to T_ex = 100 K (see Venemans et al. 2017b). Then, assuming that all carbon atoms are singly ionised, we also derived a lower limit on the total gas mass (M_gas) using the carbon abundance relative to hydrogen atom (Asplund et al. 2009) M_C/M_H = 3.54 × 10⁻³. Finally, we estimated the SFRs using the SFR − L_[CII] relation for high-redshift (z >  0.5) galaxies from De Looze et al. (2014):

(A.3)

with a systematic uncertainty of a factor of ∼2.5. In Table A.1, we list the results of spectral fits and the derived quantities for those sources with dynamical mass constrained. The reported quantities are consistent within ∼2σ to the estimates published in other works (e.g. Banerji et al. 2017; Wang et al. 2013; Decarli et al. 2017, 2018; Trakhtenbrot et al. 2017; Venemans et al. 2017c).

Appendix B: Mass support from random motions

We investigate here the turbulent pressure support term, which arises from non-rotational motions, on the total dynamical mass of our QSO host galaxies (see e.g. Epinat et al. 2009; Taylor et al. 2010). This term is not taken into account in our disc model since the gas is circularly rotating in a thin disc.

Following Epinat et al. (2009), we quantify the mass supported by random motions inside the galaxy through the virial theorem:

(B.1)

where R_D is the scale radius of the exponential brightness profile (see Sect. 4); σ₀ an estimate of constant velocity dispersion throughout the whole galaxy, and C is a parameter depending on the mass distribution and geometry. Here, we assume C = 2.25, which is the average value of known galactic mass distribution models (Binney & Tremaine 2008).

The intrinsic velocity dispersion can be estimated from the observed velocity-dispersion map after taking into account the angular and spectral resolution of observations. As a representative example, we choose SDSS J0923+0247. The velocity field of the latter shows a clear velocity gradient, and we estimated the disc inclination of ∼29 deg and dynamical mass of ∼6.0 × 10¹¹ M_⊙ (see Fig. 6). Although our rotating thin disc model reproduces the observed velocity map very well, the observed velocity dispersion is slightly (∼1.4×) higher than what was expected by our best-fit (see Fig. B.1). Indeed, the best-fitting velocity-dispersion field includes only the effect of the unresolved velocity gradients and enlargement of the emission line profile due to the beam smearing and the instrumental line-spread function (set to σ_LSF = 15 km s⁻¹). As previously mentioned, our model does not include random motions due to the physics of the gas. The intrinsic velocity dispersion of the gas can be estimated by quadratically subtracting the model from the measured velocity-dispersion map. The model and quadratic residuals are shown in Fig. B.1. Then, we computed the “1/errors”-weighted mean velocity dispersion (σ₀), where the aforementioned errors are associated with the observed velocity-dispersion values that are estimated from the spaxel-by-spaxel line fit of the data cube. The resulting value is σ₀ ∼ 54 km s⁻¹. Using this value in Eq. (B.1), we obtain M_σ ∼ 1.3 × 10⁹ M_⊙. This mass budget accounts for only ∼2% of the estimated dynamical mass obtained with kinematical modelling assuming rotating thin disc geometry. By comparing this value with the uncertainties of the dynamical mass estimate (∼0.5 dex), we conclude that the mass support from random motions is negligible in this system. We adopted the same technique for all the sources for which we obtained a dynamical mass estimate and verified that the eventual mass budget arising from random motions is always included within the dynamical mass uncertainties. Therefore, we conclude that the rotating gas disc model provides a good description for the gas kinematics of our QSO host galaxies.

Fig. B.1. Estimation of local velocity dispersion in J0923+0247. Left panel: observed velocity-dispersion field. Central panel: simulated flux-weighted velocity-dispersion map along the line of sight, corresponding to the best-fit model of the velocity field. Right panel: quadratic residuals representing the local velocity dispersion along the line of sight. The contours show the line-velocity integrated map. The contours correspond to [0.25I0, 0.5I0, I0/e, 0.68I0, 0.9I0], where I0 is the maximum value of the observed flux. In the bottom-left corner of each panel, the FHWM of the synthesised beam is shown.

We stress that, for systems that show complex or perturbed kinematics (e.g. due to the presence of outflow or a companion in the close environment of quasar, see Sect. 7), the assumption of rotating disc geometry is undoubtedly less meaningful than for relaxed systems. This is the case of, for example, SDSS J0129−0356 and ULAS J1319+0950, for which the best fits are consistent with face-on discs. On the other hand, the observed velocity dispersion is still high (∼100−200 km s⁻¹, but consistent within 2−3σ to the model), which is not expected for a face-on disc. Although this evidence supports the fact that, at least in these extreme cases, the hypothesis of a thin rotating disc is too simplified, and the dynamical mass estimates should be considered tentative; deeper observations are needed in order to properly describe the kinematics in these complex systems, rather than a more accurate kinematical modelling of the current observations. However, our fitting method enables us to homogeneously study the whole sample.

Appendix C: Comparison between dynamical masses and the spectroastrometric mass estimates

Since the dynamical masses estimated in this work are measured through a full kinematical modelling of the velocity field, they can be considered reliable mass estimates. Comparing them with mass measurements obtained with other methods enables us to test the consistency of the results. In this section, we compare the host galaxy’s dynamical mass measurements listed in Table 2 with the mass estimates obtained through the spectroastrometry method by Gnerucci et al. (2011) and the virial mass estimates obtained in Sect. 9.3.

Spectroastrometry is a technique that combines spatial and spectral resolution to probe spatial scales smaller than the angular resolution of the observations. We applied it to the case of our high-z quasar sample to estimate the product μ = M_dynsin²β. Following Gnerucci et al. (2011), we measured the FWHM and the central frequency of the line from the integrated spectra, (e.g. see Fig. 2), then we collapsed the redshifted and blueshifted channels obtaining “red” and “blue” maps, respectively (see Fig. C.1). If the galaxy disc is at least marginally resolved, the latter two maps are spatially shifted due to the rotation of the gas. Then, we performed a 2D fit of collapsed maps using an elliptical Gaussian function, and we determined the position of the two centroids. Thus, we computed the spectroastrometric radius (r_spec), as the half distance between “red” and “blue” centroids. Finally, we used the FWHM and r_spec measurements, and we estimated the spectroastrometric mass (see Eq. (2) in Gnerucci et al. 2011):

(C.1)

Fig. C.1. Integrated spectrum and red/blue maps of SDSS J0923+0247. Left panel: the data and the best-fit model is shown with black and green lines, respectively. The red dashed vertical line indicates the central frequency; the rms of the residuals is indicated by the dotted horizontal lines. Channels used to create the red and blue maps (show on the right panels) are filled with their respective colours. The central brown bin is added to the red and blue side of the collapsed images with a weight given by the fraction of red and blue bins. Right panels: the best 2D Gaussian model is shown with white contours. The blue and red circles indicate the centroid positions of the counterpart map. Here, pixels not defined in the maps are replaced with simulated noise in order to avoid numerical drawbacks in the fitting process.

where f_spec is the calibrator factor. Here, we used the value of Gnerucci et al. (2011); f_spec = 1.0 ± 0.1. The r_spec values and μ_spec are listed in Table C.1.

In Fig. C.2 (upper panel), we compare μ_spec from spectroastrometry with M_dynsin²β measured from the full kinematical modelling (see Table 2). We conclude that the two estimators are consistent within uncertainties. We also compare μ_vir = M_virsin²β from virial estimates (see Sect. 9.3) with M_dynsin²β (bottom panel of Fig. C.2). The results show that the virial mass factor μ_vir and M_dynsin²β have a non-linear relation with large dispersion. Therefore, we conclude that the classical virial method does not provide a reliable prediction of M_dynsin²β. On the other hand, although the spectroastrometry method also suffers from large uncertainties and biases, we conclude that μ_spec are in better agreement with dynamical mass factor obtained with the full kinematical modelling of galaxy discs presented in this work. In fact, as discussed by Gnerucci et al. (2011), the classical virial mass estimate can be biased by systematic errors mostly associated with the measurement of galaxy dimensions. This result confirms the reliability and usefulness of the spectroastrometry method, especially in the typical case of both poor spatial resolution and S/N ratio of the majority of the current available observations of high-z galaxies.

Fig. C.2. Comparison of galaxy mass factors obtained with spectroastrometry (upper panel) and virial estimates (bottom panel) with Mdynsin2β from full kinematical modelling of the host galaxy’s gas velocity field. The dashed black line represents the 1:1 relation. The object symbols are the same as in Fig. 12.

Appendix D: LBT observations and NIR spectra

The observations of the five quasars (SDSS J0129−0035, SDSS J2054−0005 and SDSS J2310+1855, PSO J308−21, and PSO J138+05) were executed between 2018 September and 2019 June (PI: G. Cresci) with LUCI in seeing-limited conditions using the standard strategy for near-infrared long-slit spectroscopic observations: we dithered the objects along the 1″ slit following an ABAB cycle in order to subtract the sky. We made use of the low-resolution grating (G200, λ/Δλ ≈ 2000) and the N1.8 camera (pixel size ∼0.25′ pix⁻¹) to maximise the signal-to-noise ratio. In order to obtain an accurate flux calibration of the spectra, that is required to estimate L_λ, we also obtained images of the QSOs with the J or K filter using the N3.75 camera (pixel size ∼0.12′ pix⁻¹). The total exposure time for spectroscopy is ∼2 h per target, and ∼100, 200 or 800 s for the imaging of the targets, depending on their apparent magnitude. The data were reduced and delivered by the LBT Imaging Data Center using the dedicated pipelines.

The requested time was derived assuming for the MgII line flux a typical value of 5 × 10⁻¹⁶ erg⁻¹ s⁻¹ cm⁻² at z ∼ 6 (Mazzucchelli et al. 2017), and for the CIV line a flux F_CIV ∼ 3 × F_MgII (Shen et al. 2011), and a line width of 4000 km s⁻¹. Unfortunately, we did not achieve the requested sensitivities due to bad weather conditions, and we detected BLR emission only in J2310+1855. We therefore present here the M_BH mass of J2310+1855 derived from the observations of CIV BLR line. We also note that both Feruglio et al. (2018) and Shen et al. (2019) obtained independent NIR spectra of this target with different facilities. In particular, Shen et al. (2019) published NIR spectra of a large sample z ∼ 5.7 QSOs, also providing an M_BH measurement for the J2310+1855 through virial relations based on CIV and MgII broad emission line.

Before modelling the CIV line in our LBT spectrum, we subtracted the continuum emission, fitting a power law at both sides of the ionised carbon line (in the two windows at 1450 Å and 1700 Å). Then, we used a single Gaussian model to reproduce the CIV BLR emission profile. In fact, the low S/N does not allow us to constrain the possible contribution from iron emission in the region around the CIV, which is expected to be negligible (see e.g. Shen et al. 2008, 2011), nor the possible emission from the CIV NLR line (e.g. Shen et al. 2011).

From the best fit of the CIV line, we derived a FWHM ∼ 12 500 km s⁻¹, and from the extrapolated continuum at 1350 Å, a flux of F₁₃₅₀ ∼ 1.1 × 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹, and a luminosity of ∼10^45.8 erg s⁻¹. Using the Vestergaard & Peterson (2006) relation, we obtained log M_BH ∼ 9.8, consistent with Shen et al. (2019). The CIV line is blueshifted with respect to the [CII]_158 μm systemic of Δv ∼ −7200 km s⁻¹, strongly suggesting the presence of outflows in this source. We note that the values of FWHM and Δv of CIV broad line in the J2310+1855 spectrum are consistent with the typical values estimated in high-z QSOs (see e.g. Vietri et al. 2018). The CIV-based M_BH estimate can be therefore strongly biased; by adopting the different calibrations introduced to correct for the outflow contribution in CIV lines (see e.g. Vignali et al. 2018, and references therein), we obtained mass estimates in the range log M_BH = 8.9−9.4, which is in agreement with the estimate by Feruglio et al. (2018). The uncertainties on these measurements are dominated by the intrinsic scatter (≈0.3 dex; see e.g. Vestergaard & Peterson 2006; Denney 2012; Park et al. 2017) in the single-epoch calibrations, which are much larger than the typical uncertainties ascribed to the measurements of the line widths and fluxes.

The latter values are consistent with the MgII-based M_BH reported in Shen et al. (2019), which is used in this paper to study the M_BH−M_dyn relation. In Fig. D.1, we show the NIR spectrum of J2310+1855 with the best fit of the CIV broad line.

Fig. D.1. Portion of LBT/LUCI zJ spectrum of SDSS J2310+1855 around CIV line for which the expected wavelength is indicated in the figure according to [CII]-based redshift. The orange and black dashed curves indicate the spectrum’s best fit (line + continuum and continuum, respectively).

Appendix E: Maps, integrated spectra, and the results of the kinematical modelling

Here, we report the integrated spectra, flux, velocity, and velocity-dispersion maps for objects in Table 2 (Fig. E.1). See Figs. 2 and 3 for the descriptions of each panel. We also report the 2D best fit of the flux and velocity maps (Fig. E.2). The different panels are labelled as they are in Fig. 6; we refer to the latter for a description of the figures.

Fig. E.1. Line integrated spectra, flux, velocity, and velocity-dispersion maps along the line of sight for objects listed in Table 2. See Figs. 2 and 3 for the description of each panels.

Fig. E.1. continued.

Fig. E.2. Results of 2D kinematical modelling for objects listed in Table 2. We refer to Fig. 6 for descriptions of individual panels.

Fig. E.2. continued.

Fig. E.2. continued.

Fig. E.2. continued.

Fig. E.2. continued.

Fig. E.2. continued.

Fig. E.2. continued.

All Tables

All Figures

	Fig. 1. Redshift distribution of 32 QSOs listed in Table 1. The histogram colours indicate the emission line detected in the ALMA data. The figure also shows the sample cumulative redshift distribution (dotted line and right axis). Individual sources are marked with blue circles (sizes are proportional to the ALMA synthesised beam FWHM).
In the text

Fig. 2. Integrated spectrum of source SDSS J0923+0247. Top panel: the observed data is shown in light yellow with error bars in grey (rms in each channel). The red solid curve is the best-fit Gaussian model. The velocity scale (top axis) has as its referecence the central frequency of the best fit. Bottom panel: the fit residuals (model-data), the grey filled area shows the rms along the spectral axis.

In the text

	Fig. 3. Observed maps of SDSS J0923+0247. From left panel to right: flux map, velocity map, and velocity-dispersion map along the line of sight, respectively. At the bottom-left corner of each panel, we report the ALMA synthesised beam FWHM. The coordinates indicate offsets with respect to the map centre.
In the text

	Fig. 4. Red and blue residual maps of SDSS J0923+0247. This test reveals spatially resolved kinematics consistent with what we would expect from a rotating disc.
In the text

Fig. 5. Simulated maps obtained with the kinematical model for a galaxy thin disc with an exponential brightness profile defined by a scale radius RD = 0.125″ and dynamical mass Mdyn = 5.0 × 1010 M⊙. The flux, velocity, and velocity-dispersion maps along the LOS are indicated by F(xp, yp′), V(xp, yp′), and σ(xp, yp′), respectively. From top to bottom panel: disc inclination and position angle (β, γ) are, respectively, equal to (40, −45) deg, (40, 0) deg, and (60, 45) deg. The other parameters defining the model are: the FWHM of the synthesised beam (0.4″ × 0.275″), the position angle of the beam (BPA = −50 deg), the FWHM of the LSF (σLSF = 20 km s−1), the bin size (0.05″ × 0.05″), and the angular radius of the disc (R = 1.25″). The number of clouds used in this model is Np = 7.5 × 106.

In the text

Fig. 6. Kinematical modelling performed on SDSS J0923+0247. Upper panels [a.1] and [a.2]: fit result of the flux map, bottom panels [b.1] and [b.2]: result of the 2D velocity field fit. [a.1]: curve of χ2 minima as a function of the disc inclination (see Sect. 4.2 for details). The magenta hexagon indicates the absolute minimum (=35 deg) of χm(β). The green area (see Sect. 4.3 for the definition) indicates the allowed inclination values used as a prior in performing a 2D velocity field fit. [a.2]: 2D best-fit model of the flux map. From left to right, we report the observed map, the model, and the residuals ([data-model]/error). Residuals values are color-coded according to the horizontal colorbar displayed at the bottom-left corner of the panel. [b.1]: the posterior probability distributions of the free parameters in 2D velocity field fits retrieved with the MCMC algorithm with the best values and their uncertainties. [b.2]: 2D best-fit model of the kinematical map. From left to right, we report the observed velocity map, the model, and the velocity curves extracted from a long-slit of two pixels in width aligned with the line of nodes. The slit is superimposed on maps; red circles and solid black lines in the right panel indicate the observed flux-weighted velocity values in each bin of the slit and the model, respectively.

In the text

Fig. 7. Simulated velocity field corresponding to four different couples of inclination angle β and dynamical mass Mdyn = 5.0 × 1010/​​sin2β. Increasing the dimension of the beam FWHM (shadowed area), the iso-velocity curves are increasingly smoothed, and different configurations of the disc appear almost indistinguishable. The simulated maps are obtained with the kinematical model for a galaxy thin disc with RD = 0.125″, γ = −90 deg. The other parameters defining the model are set as equal to the model in Fig. 5.

In the text

Fig. 8. Summary of the different sub-selections used throughout this work. Starting from an initial sample of 72 quasar host galaxy observations extracted from the ALMA data archive, the final sample is composed of ten high-z objects for which we study the BH-galaxy relation. The characters indicating the selection type are in accordance with Table 1. The number of objects rejected from the sample is also indicated at each step.

In the text

Fig. 9. High-redshift relation between the black hole mass (MBH) and the dynamical mass of the host galaxy (Mdyn). The dashed black line and the dotted red line represent the reference local relation inferred using samples of local galaxies (E = ellipticals, S/S0/Sb = spirals) shown as black dots (Kormendy & Ho 2013, also indicated as KH+13) and red triangles (de Nicola et al. 2019, or DN+19). The green line is the relation found by Reines & Volonteri (2015, or RV+15) by measuring the total stellar mass in a sample of the local AGN (green crosses). The solid blue line is the best fit to our data. The shadowed areas show the 1σ uncertainty. In the case of SDSS J0129−0035 and ULAS J1319+0950, we inferred a lower limit on the dynamical mass. We do not take these data into account in the fit. The circles of our data points indicate the sources for which the BH masses are estimated from bolometric luminosity assuming Eddington accretion.

In the text

Fig. 10. Evolution of Γ = MBH/Mdyn as a function of redshift z. The black dotted and dash-dotted lines represent, respectively, the relation found by Decarli et al. (2010) at z ≲ 3, and the corresponding ratio at z = 0. The inset panel shows the same plot at 2 <  z <  3. In the case of SDSS J0129−0035 and ULAS J1319+0950 we inferred lower limits on dynamical masses (i.e. an upper limit on the ratio Γ). The circles indicate those sources which BH masses are estimated from bolometric luminosity assuming Eddington accretion.

In the text

	Fig. 11. Comparison between our dynamical mass measurement distribution and galaxy stellar mass function. Upper panel: Galaxy stellar mass function at different redshifts (Song et al. 2016). The dashed vertical lines show the position of the M⋆ values that represent the “knee” of the mass function. Lower panel: stacked histogram of our dynamical mass estimates.
In the text

Fig. 12. Comparison between dynamical virial mass (Mvir) computed using flux map properties and dynamical mass (Mdyn) obtained with the full kinematical model of the gas velocity field. The dashed line indicates the 1:1 relation, while dotted lines are 1 dex shifts. The objects represented with star symbols are upper limits on virial mass, since their emission appears consistent with a point-like source (as a result of 2D Gaussian fits within CASA). In these cases, following Willott et al. (2015a), we assume β = 55 deg as disc inclination angle. The black square indicates PSO J183−05, for which we retrieved a lower limit on Mdyn, but for which MBH is not available.

In the text

Fig. B.1. Estimation of local velocity dispersion in J0923+0247. Left panel: observed velocity-dispersion field. Central panel: simulated flux-weighted velocity-dispersion map along the line of sight, corresponding to the best-fit model of the velocity field. Right panel: quadratic residuals representing the local velocity dispersion along the line of sight. The contours show the line-velocity integrated map. The contours correspond to [0.25I0, 0.5I0, I0/e, 0.68I0, 0.9I0], where I0 is the maximum value of the observed flux. In the bottom-left corner of each panel, the FHWM of the synthesised beam is shown.

In the text

Fig. C.1. Integrated spectrum and red/blue maps of SDSS J0923+0247. Left panel: the data and the best-fit model is shown with black and green lines, respectively. The red dashed vertical line indicates the central frequency; the rms of the residuals is indicated by the dotted horizontal lines. Channels used to create the red and blue maps (show on the right panels) are filled with their respective colours. The central brown bin is added to the red and blue side of the collapsed images with a weight given by the fraction of red and blue bins. Right panels: the best 2D Gaussian model is shown with white contours. The blue and red circles indicate the centroid positions of the counterpart map. Here, pixels not defined in the maps are replaced with simulated noise in order to avoid numerical drawbacks in the fitting process.

In the text

	Fig. C.2. Comparison of galaxy mass factors obtained with spectroastrometry (upper panel) and virial estimates (bottom panel) with Mdynsin2β from full kinematical modelling of the host galaxy’s gas velocity field. The dashed black line represents the 1:1 relation. The object symbols are the same as in Fig. 12.
In the text

	Fig. D.1. Portion of LBT/LUCI zJ spectrum of SDSS J2310+1855 around CIV line for which the expected wavelength is indicated in the figure according to [CII]-based redshift. The orange and black dashed curves indicate the spectrum’s best fit (line + continuum and continuum, respectively).
In the text

	Fig. E.1. Line integrated spectra, flux, velocity, and velocity-dispersion maps along the line of sight for objects listed in Table 2. See Figs. 2 and 3 for the description of each panels.
In the text

	Fig. E.1. continued.
In the text

	Fig. E.2. Results of 2D kinematical modelling for objects listed in Table 2. We refer to Fig. 6 for descriptions of individual panels.
In the text

	Fig. E.2. continued.
In the text

	Fig. E.2. continued.
In the text

	Fig. E.2. continued.
In the text

	Fig. E.2. continued.
In the text

	Fig. E.2. continued.
In the text

	Fig. E.2. continued.
In the text