Resolving Disks&Mergers in $z\sim2$ Heavily Reddened Quasars and their Companion Galaxies with ALMA

We present sub-arcsecond resolution ALMA imaging of the CO(3-2) emission in two $z\sim2.5$ heavily reddened quasars (HRQs) - ULASJ1234+0907 and ULASJ2315+0143 - and their companion galaxies. Dynamical modeling of the resolved velocity fields enables us to constrain the molecular gas morphologies and host galaxy masses. Combining the new data with extensive multi-wavelength observations, we are able to study the relative kinematics of different molecular emission lines, the molecular gas fractions and the locations of the quasars on the M$_{\rm{BH}}$-M$_{\rm{gal}}$ relation. Despite having similar black-hole properties, the two HRQs display markedly different host galaxy properties and local environments. J1234 has a very massive host, M$_{\rm{dyn}} \sim 5 \times 10^{11}$M$_\odot$ and two companion galaxies that are similarly massive located within 200 kpc of the quasar. The molecular gas fraction is low ($\sim$6%). The significant ongoing star formation in the host galaxy is entirely obscured at rest-frame UV and optical wavelengths. J2315 is resolved into a close-separation major-merger ($\Delta$r=15 kpc; $\Delta$v=170 km/s) with a $\sim$1:2 mass ratio. The total dynamical mass is estimated to be $\lesssim$10$^{11}$M$_\odot$ and the molecular gas fraction is high ($>$45%). A new HSC image of the galaxy shows unobscured UV-luminous star-forming regions co-incident with the extended reservoir of cold molecular gas in the merger. We use the outputs from the Illustris simulations to track the growth of such massive black holes from $z\sim6$ to the present day. While J1234 is consistent with the simulated $z\sim2$ relation, J2315 has a black hole that is over-massive relative to its host galaxy.


INTRODUCTION
In our current paradigm of galaxy formation, accreting supermassive black-holes play a fundamental role in shaping the evolution of the most massive galaxies. Gas-rich major mergers have for some time been invoked in theoretical models to fuel both star formation and accretion onto the black hole in massive galaxies (e.g. Di Springel et al. 2005;Hopkins et al. 2008;Narayanan et al. 2010;Gabor et al. 2016;Steinborn et al. 2018). Idealised as well as cosmological simulations of merging galaxies now have the requisite spatial and temporal resolution to be able to demonstrate that the merger leads to enhanced star formation (e.g. Moreno et al. 2015), triggering of an active galactic nucleus (AGN; e.g. Capelo et al. 2015;Angles-Alcazar et al. 2020), and an alteration ★ E-mail: M. Banerji-Wright@soton.ac.uk in the interstellar medium properties of the galaxy (e.g. Moreno et al. 2019). Energy injection from the quasar during a feedbackdominated phase (Fabian 2012) is then thought to be responsible for quenching star formation and curtailing the growth of these massive galaxies. Cosmological simulations routinely implement such feedback in order to explain the observed distributions of massive galaxies (e.g. Sĳacki et al. 2015;Beckmann et al. 2017). A major aim of current observational research is to find compelling evidence in favour of this standard paradigm. Recently several aspects of this simple picture have been confronted by new observational data. For example, the existence of massive gas-disks in the high-redshift Universe (e.g. Tacconi et al. 2013;Förster Schreiber et al. 2018) provides evidence that processes such as smooth accretion via cold gas flows (Dekel et al. 2009) or minor mergers may play a more prominent role in assembling galaxies even at the high-mass end of the galaxy mass function, therefore challenging the major-merger paradigm. The effects of AGN feedback too can be more nuanced than previously thought. For example, theoretical models predict that compression of gas in galactic outflows can indeed trigger rather than inhibit star formation (e.g. Ishibashi et al. 2017). There have been some claims of observations of such positive feedback effects (e.g. Cresci et al. 2015;Carniani et al. 2016) with evidence now emerging that stars can in fact be formed in-situ within a galactic outflow Gallagher et al. 2019).
In the context of the merger-driven evolutionary scenario for the formation and co-evolution of massive galaxies and their supermassive black-holes, populations of obscured quasars could represent an important phase as a galaxy is transitioning from a dustobscured starburst to an ultraviolet luminous quasar (e.g. Hickox & Alexander 2018). The idea has its origins in the seminal work of Sanders et al. (1988). With the advent of wide-area surveys in the X-ray, near and mid infrared wavelengths, population studies of the luminous, obscured quasar population are now possible (e.g. Banerji et al. 2012Banerji et al. , 2013Banerji et al. , 2015Glikman et al. 2012;Assef et al. 2015;Hamann et al. 2017;LaMassa et al. 2017;Temple et al. 2019). The most luminous obscured quasars detected in these wide-area surveys are often seen to be hosted by merging galaxies (e.g. Urrutia et al. 2008;Glikman et al. 2015;Fan et al. 2016;Bischetti et al. 2018;Díaz-Santos et al. 2018) with significant amounts of ongoing star-formation (e.g. Banerji et al. 2014Banerji et al. , 2017Urrutia et al. 2012;Wethers et al. 2018). Powerful AGN-driven outflows have also been detected in several of these systems (e.g. Urrutia et al. 2009;Perna et al. 2015;Cresci et al. 2015;Zakamska et al. 2016;Temple et al. 2019;Kakkad et al. 2020), with the observations suggesting that these outflows could have both a positive and negative feedback effect on star formation.
Luminous quasars were among the first high-redshift ( 2), extragalactic sources to be targeted in molecular gas emission (e.g. Solomon & Vanden Bout 2005 and references therein) and submillimeter and millimeter observations provide us with an opportunity to study host galaxy emission from these quasars over a wavelength range where the emission from the quasar itself is sub-dominant. Molecular gas has been detected in quasars out into the Epoch of Reionisation (e.g. Walter et al. 2003Walter et al. , 2004Maiolino et al. 2005;Wang et al. 2010Wang et al. , 2011Venemans et al. 2017) and the new capabilities of the Atacama Large Millimeter Array (ALMA) are now making it possible to spatially resolve both the dust and gas emission in high-redshift quasar hosts (e.g. Carilli et al. 2013;Wagg et al. 2014;Salomé et al. 2012;Carniani et al. 2013;De Breuck et al. 2014;Venemans et al. 2017;Trakhtenbrot et al. 2017;Shao et al. 2017;Brusa et al. 2017;Fogasy et al. 2020). Nevertheless, the number of high-redshift quasars with such resolved measurements of the molecular gas emission are still relatively small. Resolved measurements of the host galaxies in the millimeter wavelengths enable unique constraints on the dynamical masses, gas morphologies and gas fractions in such systems.
Observational evidence is mounting that many of the major physical processes thought to govern massive galaxy formation -e.g. mergers, starbursts, high-accretion black-hole growth and outflowsare synchronously occurring during the dust-obscured quasar phase. How are these processes triggered and what effect do they have on the quasar host galaxies? We address these questions here via detailed, spatially resolved observations of the CO(3-2) gas emission in the quasars ULASJ1234+0907 ( ∼ 2.50) and ULASJ2315+0143 ( ∼ 2.56). The two quasars are very similar in terms of the AGN properties. Both quasars are extremely luminous (L 6 m ∼ 10 47 erg/s; L X ∼ 10 45 erg/s -see Lansbury et al. 2020) with massive black-holes (M BH > 10 9 M ) obscured by significant columns of dust and gas (N H ∼ 10 22 cm −2 ; A V ∼ 5 − 6 mag; Banerji et al. 2012Banerji et al. , 2015Banerji et al. , 2017Wethers et al. 2018;Lansbury et al. 2020). Both quasars are known to reside in highly star-forming host galaxies (Banerji et al. 2014;Wethers et al. 2018Wethers et al. , 2020. The first detections of molecular gas via the CO(3-2) emission in the two quasars was presented in Banerji et al. (2017) (B17 hereafter). In Banerji et al. (2018) (Paper I hereafter) we combined the CO(3-2) detections with detections of CO(1-0), CO(7-6) and CI(2-1) to better constrain the interstellar medium properties of the quasar host galaxies. While Paper I focused on the integrated line emission properties of our quasar host galaxies, in this paper we concentrate on studying the resolved emission to gain a better understanding of the physical processes that govern the formation of these galaxies. Throughout the paper we assume a flat, ΛCDM cosmology with H 0 =70 km s −1 Mpc −1 , Ω M =0.3 and Ω Λ =0.7.

ALMA BAND 3 OBSERVATIONS
Our previous work (B17) presented detections of the CO(3-2) line in a sample of four heavily reddened quasars (HRQs) at ∼ 2.5 that were initially discovered in Banerji et al. (2012Banerji et al. ( , 2015. These four quasars are among the reddest quasars found in our nearinfrared search for luminous, obscured quasars that are not present in wide-field optically selected spectroscopic quasar samples such as those from the Sloan Digital Sky Survey (SDSS) and BOSS. Briefly, the quasars have extremely red near infra-red colours of ( − ) Vega > 2, corresponding to dust extinctions towards the quasar continuum of ( − ) 1 (A V 3 mag for R V = 3.1 mag). The bolometric luminosities and black-hole masses are comparable to the most luminous quasars selected in optical surveys such as the Sloan Digital Sky Survey. As a detection experiment, the Cycle 3 ALMA observations in B17 made use of the most compact available configuration of the telescope leading to an angular resolution of ∼2.5-3 . Two of the quasars from the sample -ULASJ1234+0907 ( ∼ 2.50; J1234 hereafter) and ULASJ2315+0143 ( ∼ 2.56; J2315 hereafter) -were marginally resolved in the Cycle 3 data and were therefore re-observed during ALMA Cycle 4 at ∼3× better angular resolution with the aim of spatially resolving the gas emission in both sources and constraining the gas dynamics. J1234 has two nearby gas-rich galaxies that were first discovered via their CO(3-2) detections in B17. Both companion galaxies (named G1234N and G1234S) were also imaged at higher resolution as part of the observations presented in this work. Further details of the new ALMA Cycle 4 observations can be found in Table 1.
The spectroscopic observations of both quasar systems made use of a correlator configuration with four dual polarisation bands, each of bandwidth 1.875GHz. In the case of J1234, observations made use of the Time Division Mode (TDM) of the correlator, providing a channel width of 15.625 MHz. In the case of J2315, the Frequency Division Mode (FDM) of the correlator was instead used, providing a much finer channel width of ∼488 kHz across the entire bandwidth, which we then binned up to 15.625 MHz resolution.
The raw data were processed using (v4.7.0), by executing the appropriate calibration scripts for each set of observations as provided by the ALMA Observatory. Line cubes were produced after u-v plane continuum subtraction employing the line-free channels. We then used the clean algorithm with a natural weighting of the visibilities to transform to the image plane. Continuum images were also produced by running clean in mfs mode on the line-free spectral windows, and once again using a natural weighting. All im-ages have had a primary beam correction applied to them using the task pbcor. Images have been generated with a pixel-scale of 0.1 per pixel. The beam sizes as well as RMS noise levels reached by the data for each of the sources is summarised in Table 1. The spectra have all had Hanning smoothing applied to them as part of the pipeline processing.

OBSERVED MOLECULAR GAS AND DUST CONTINUUM PROPERTIES
In Fig. 1 and 2 we present the CO(3-2) images as well as integrated spectra for all galaxies in the J1234 and J2315 systems respectively obtained via our ALMA Cycle 4 program. We infer initial source sizes by fitting a 2-dimensional Gaussian to the zeroth moment image using . The spectra have been extracted using apertures of variable sizes (diameter 1-2 ) centred on the positions of each galaxy stated in Table 2. The aperture sizes chosen depend on the size of the CO emitting region for each galaxy, which is also quoted in Table 2 and discussed further in Sections 5 and 6. These new observations provide higher signal-to-noise ratio (S/N) spectra for all sources compared to the Cycle 3 data presented in B17, as well as higher S/N measurements of the continuum flux densities at observed frequencies of ∼86-103 GHz. These measurements are summarised in Table 2.
All three galaxies in the J1234 system are spatially resolved in our data based on the fact that CASA is able to deconvolve the intrinsic emission from the convolved image. We have also checked that the radial profiles of all three galaxies are more extended than would be expected for a point source. The spectra exhibit clear double-peaked CO(3-2) emission line profiles. Such clear doublepeaked profiles were not evident in the original Cycle 3 spectra in B17, and indicate that we have now resolved multiple kinematic components in the gas emission. The J2315 system, which was already spatially resolved in CO emission in B17, now breaks up into two distinct galaxies in a close-separation major merger in our new data. We are therefore able to extract CO(3-2) spectra for both the quasar host and the companion galaxy lying to the North of the quasar. Table 2 presents the line and continuum flux densities for both our quasars and their companion galaxies. These measurements have been analysed in more detail in Paper I, where we combine them with observations of the CO(1-0), CO(7-6) and CI(2-1) lines as well as the far infra-red to radio continuum in order to infer the properties of the interstellar medium in these galaxies. In Paper I we find that the quasar J1234 as well as its companion galaxy G1234S are comprised of highly excited, dense gas consistent with what is seen in luminous AGN. J2315 on the other hand has excitation properties that are more similar to that observed in star-forming galaxies. Both the quasars J1234 and J2315 as well as their companion galaxies were inferred to be highly starforming based on their dust continuum detections in Paper I.
The focus of the current paper is on a detailed kinematic analysis of the CO(3-2) emission and inferred constraints on the gas dynamics and gas morphologies in our two quasar systems. Before turning to such an analysis we first present a description of the dynamical models we will employ.

3D Barolo
We model the CO(3-2) traced gas dynamics in our quasars and their companion galaxies using the tilted ring model fitting code 3D Barolo -3D-Based Analysis of Rotating Objects from Line Observations (Di Teodoro & Fraternali 2015). The approach is similar to our previous work (B17) employing the GIPSY rotcur task, but with the key improvement that 3D Barolo fits models to the entire data cube rather than just the velocity field. By convolving the model with the restoring beam and comparing data and model on a channel-by-channel basis, the effects of beam smearing are reduced, and additional dynamical aspects may be examined. Such an approach has been used previously to study the rotation curves of nearby galaxies (e.g. de Blok et al. 2008;Iorio et al. 2017;Shelest & Lelli 2020) and is now being routinely applied at high-redshifts (e.g. Jones et al. 2017;Xue et al. 2018;Tadaki et al. 2020).
The three-dimensional model fitting considers some of the same variables that are included in two-dimensional fits: inclination, i, position angle, PA, rotational velocity, rot , systemic redshift, sys and the spatial centroid, [ 0,0]. By expanding to the spectral dimension, we can also consider velocity dispersion , radial brightness profile ( ), and scale height as free parameters. While 3D Barolo has been well-tested on low-S/N and low-resolution datasets (e.g., Di Teodoro & Fraternali 2015), the best-fit results are strongly dependent on both the model geometry (i.e., number of rings, ring width) and initial estimates for each parameter. To reduce any implicit bias, we automate the fitting process by using signal isolation and moment maps.
We first utilize the 3D Barolo SEARCH utility (based on DUCHAMP; Whiting 2012) to create a signal mask. This utility automatically finds the RMS noise level of an input data cube and identifies signal peaks above a user-provided signal-to-noise threshold (S/N upper = 3.5). A search is then conducted around these peaks for emission above a second user-provided S/N lower = 3.0, creating a three-dimensional mask of significant emission. This approach includes fewer non-signal pixels than a 2-D spaxel mask or a mask based on a single S/N threshold.
Next, we use the CASA task immoments to create moment 0 (integrated intensity) maps using all channels containing significant line emission. By fitting this map with a two-dimensional Gaussian (CASA task imfit), we find a best-fit beam-deconvolved FWHM of the major and minor axes, position angle, and morphological central position. We also apply immoments to all pixels identified as signal to create moment 1 (line of sight velocity) and moment 2 (velocity dispersion) maps. These maps can be seen in the left most panels of Fig. 3, 4, 5 and 8.
With this mask and these maps in hand, we are able to provide physically motivated estimates for most model parameters. We base the maximum model radius on the 2-D Gaussian fit to the moment zero map: 1.0×(the deconvolved FWHM of the major axis). In order to not over-resolve the data, we adopt a minimum ring width of (FWHM of the minor axis of the restoring beam)/2.5 -see also Talia et al. (2018); Fan et al. (2019). The number of model rings is then found by dividing the maximum model radius by the minimum ring width and rounding down to the nearest integer, while the ring width is taken to be the maximum model radius divided by this number. For all galaxies modelled in this work, this process returns models with three rings. The central position of the model is fixed to be the centroid of the 2-D Gaussian fit to the moment zero map. Since each emission line is well-detected, the systemic velocity is expected to be 0 km s −1 , with respect to the redshifts listed in Table  Table 1  The right-hand panels show the integrated CO(3-2) spectra for the three galaxies together with the best-fit double Gaussian fit to these line profiles. Double-peaked profiles are evident for all three galaxies. Zero velocity, indicated by the dotted vertical line corresponds to the systemic redshifts of these galaxies obtained from our ALMA Cycle 3 observations in B17.

Figure 2.
Top: Map of the CO(3-2) emission in the J2315 system. North is up and East is to the right. The source is clearly resolved into two spatial components and the quasar position is indicated by the star. The synthesized beam size is indicated by the blue ellipse at the bottom-right of the image. Bottom: CO(3-2) spectra extracted over the entire system (left), at the quasar position (middle) and at the position of the northern companion galaxy (right) together with the best-fit single Gaussian used to fit these line profiles. Zero velocity, indicated by the dotted vertical line corresponds to the systemic redshift obtained from our ALMA Cycle 3 observations in B17.
3. 3D Barolo provides an estimate of the rotational velocity based on the moment 1 map, while we provide the value of the moment 2 map at the morphological center as an estimate of the velocity dispersion. The inclination estimate is set to 45 • , and allowed to vary between 10 − 80 • , while the position angle is provided by the user. We assume a thin disk with a scale height of 0.01 ∼ 100 pc.
Using this set of initial guesses for each of these variables and a list of ring radii and widths, 3D Barolo first constructs a model of the first ring by randomly populating a space with discrete emitting gas clouds. These clouds are then given velocities to replicate the estimated rotation curve and velocity dispersion. The resulting data cube is convolved with the synthesized beam of the observation, and the summed intensity in each pixel along the velocity axis of the model cube is rescaled to equal that of the data cube. The pixels of the data and model cubes are then compared, and each parameter is adjusted. When the residual is minimized, the fitting stops, and the spatial parameters (i, PA, sys ) are averaged across all rings. The fitting is then repeated for the remaining parameters. When this reaches a minimum, the next ring is analyzed. The final outputs for this model fitting process are a model data cube and best-fit morphological (i.e., , ) and kinematic (i.e., rotational velocity curve, velocity dispersion curve, systemic redshift) parameters. The output model cube encodes information about the rotation curvei.e. rot as a function of radius, , from which we can derive the dynamical mass of the galaxy as follows: The uncertainties on the dynamical mass only include the output rotational model uncertainties from 3D Barolo ( rot,max ) and the width of the model ring ( max ) and the significant uncertainty on the inclination is not included. Since the fitting process is threedimensional and multi-staged, the inclusion of the inclination uncertainty in the error budget is non-trivial and we simply note that the formal uncertainties on the dynamical masses from 3D Barolo are likely to be under-estimated.

Exponential Disk Modelling
3D Barolo makes no assumption regarding the underlying density distribution and mass profile of the galaxy and simply fits the CO emitting material across the extent over which it is seen in our ALMA data. We therefore use an independent method to check the mass estimates by approximating the kinematics of the galaxies as a thin rotating disk (Carniani et al. 2013;Pensabene et al. 2020). Such a model assumes that the surface brightness of the disk has an exponential profile Σ ∝ − / , where is the scale radius that corresponds to a half-light radius of 1/2 = 1.67 × . By assuming that the mass surface density has the same profile as the gas emission, the circular velocity associated with the exponential profile can be calculated following equation 2.157 of Binney & Tremaine (2008). Both intensity and velocity model maps are then matched to the spatial resolution of the observations by taking into account the observed synthesized beam. The free parameters of this dynamical model are the spatial coordinates of the disk centre, the half-light radius (or scale radius), inclination, position angle, dynamical mass, and systemic velocity. The dynamical mass is the disk mass enclosed within a radius of 10 kpc. Because this radius is much larger than the typical disk scale radius at ∼ 2 (Cresci et al. 2009), this mass corresponds to the total baryonic mass in the disk. Note also that this radius is much larger than the radius over which the CO emission is visible, and therefore the radius over which 3D Barolo fits are conducted. The full parameter space is explored by using the Markov Chain Monte Carlo (MCMC) algorithm emcee (Foreman-Mackey et al. 2013)), which allows one to obtain a probability distribution of all free-parameters along with marginalised uncertainties. The disk centre coordinates and halflight radius are derived from modelling the intensity maps, while the probability distributions of the other parameters are driven from the fitting of the velocity maps.

Dynamical Masses, Sizes & Gas Fractions
The quasar host galaxy J1234 ( = 2.501) and its two CO-luminous companions -G1234N ( = 2.514; ∼170 kpc projected distance from the quasar) and G1234S ( = 2.497; ∼90 kpc projected distance from the quasar) -are from hereon referred to as the J1234 system. The CO images shown in Fig. 1 have been constructed by collapsing all the spectral channels containing the line in each source. Initial source sizes for the gas emission were derived by fitting a single 2-dimensional Gaussian to the CO images using the task imfit. Using these constraints as an initial guess, the fitting was then re-done in the u-v plane using uvmodelfit. The results from the two fitting approaches are consistent for all the galaxies and the full-width-half-maximum (FWHM) of the major axis is (0.48±0.14) , (0.57±0.18) and (0.53±0.24) for ULASJ1234, G1234N and G1234S respectively, after deconvolving from the beam. The FWHM values can be converted to a half-light radius using the equation, R 1/2 =0.5×FWHM for a Gaussian profile, and correspond to physical scales of ∼2 kpc at the redshifts of these galaxies. We summarise the half-light radius constraints for all three galaxies in Table 3. A second estimate of the half-light radius can be estimated from fitting the exponential disk models (Section 4.2) to the intensity maps. The half-light radii derived from the exponential disk models are also specified in Table 3 and are generally larger than those derived from the Gaussian fits -∼3-4 kpc at the redshifts of these galaxies. This is to be expected given the 2-D Gaussian profiles are not able to reproduce the full spatial extent of the CO emission seen from these galaxies.
We have already noted the double peaked line profiles for all three sources seen in Fig. 1, which provide evidence for kinematically distinct components of the gas emission. Two Gaussian components are required to adequately model the 1-dimensional spectral line profiles and these best spectral-fit models are over-plotted on the observed spectra in Fig. 1. The velocity separation between the two Gaussian components is 510±160 km/s, 340±50 km/s and 350±40 km/s for the quasar, G1234N and G1234S respectively. Such double-peaked profiles could arise from the approaching and receding sides of a rotating gas disk, a coalescing merger of two distinct galaxies, or large-scale inflows and outflows affecting the molecular gas.
As can be seen in the moment 1 maps in Fig. 3, 4 and 5, a velocity gradient is evident across the full spatial extent of all three galaxies that comprise the J1234 system. The best-fit 3D Barolo models together with the fit residuals can also be seen in the same figures. In general the fit residuals are small, at least in the central portions of these galaxies, suggesting that a smooth velocity gradient provides a reasonable approximation to these data. The resulting dynamical masses are summarised in Table 3. With dynamical masses of ∼ 10 11 M , all three galaxies in the J1234 system can be considered among the most massive systems seen at these epochs.
We also derive dynamical masses for all three galaxies using the exponential disk modelling. The resulting dynamical masses summarised in Table 3 are larger than those obtained from 3D Barolo by 0.9 dex, 0.1 dex and 0.3 dex for the three galaxies. As stated earlier, the exponential disk models are sensitive to mass over a much larger physical radius than the 3D Barolo models. Furthermore,  the inclination of these galaxies is poorly constrained and adds significant additional uncertainty to the mass estimates. Given these factors, the discrepancies in the dynamical mass estimates from the two methods is not surprising.
In Fig. 6 we show the dynamical mass versus the half-light radius of the CO emission from the exponential disk models for all three galaxies in the J1234 system and compare to both field and cluster galaxy samples from the literature at ∼ 1 − 2.5. The field galaxies from the PHIBSS sample correspond to CO(3-2) measurements (Tacconi et al. 2013) whereas the ∼ 1.6 cluster galaxy measurements are taken from the CO(2-1) observations presented in Noble et al. (2018). We take the galaxy mass as the sum of the stellar mass and molecular gas mass in the case of the literature samples and as the dynamical mass in the case of our quasar host galaxy and its companions. In general, all three galaxies in the J1234 system are at the more massive and more compact end of the distribution of ∼ 2 galaxies from the literature. We also show the disk half-light radii for ∼ 4 − 6 quasars from Pensabene et al. (2020), which uses the same disk modelling as we have employed in this work. The higher redshift quasars are generally found to be more compact in size but this likely reflects the differing spatial distribution of the [CII] gas in the case of the ∼ 4 − 6 QSOs and the relatively low-J CO gas traced by the ∼ 2 samples.
In Paper I, we derived gas masses for all galaxies in the J1234 system using several different estimators, thus mitigating against the uncertainties implicit in using a single molecular tracer of the total gas mass. The resulting gas masses were found to be ∼ 3 × 10 10 M , ∼ 2 × 10 10 M and ∼ 1 − 7 × 10 10 M for ULASJ1234, G1234N and G1234S respectively. In the case of G1234S the large range results from a much larger CI-derived gas mass compared to the CO-derived gas mass (Paper I). Combining these gas mass estimates with the dynamical masses and assuming that the inner parts of the disks traced by our CO observations are indeed baryondominated, we can infer the molecular gas fractions. The gas fraction is traditionally defined as the ratio of the gas mass to the total stellar mass. In the absence of stellar mass estimates for our highluminosity quasars, we infer the stellar masses from the dynamical masses via M * =M dyn −M gas −M BH , as has been done in previous studies of high-luminosity quasars (e.g. Bischetti et al. 2021). In the case of the quasar, very different gas fractions are implied depending on whether we consider the dynamical mass from 3D Barolo, which represents the mass in the innermost regions of the galaxy, which are seen in CO emission, or whether we consider the exponential disk model mass, which extends out to 10kpc. In the former case, we obtain gas fractions of ∼85% whereas the gas fraction relative to the extended disk model mass is more like ∼6%. It is reasonable to assume that the galaxy is baryon dominated even out to 10kpc, in which case there must be a substantial stellar component already in place with M * ∼ 10 11 M in this quasar. G1234N has a gas fraction of ∼14-18% whereas in G1234S the gas fraction is more uncertain owing to the discrepant gas masses from the CO and CI observations, and could be in the range 7-90%.
The star formation rates for the two companion galaxies are poorly constrained due to a lack of photometric data points sampling the peak of the galaxy dust SED. However, for the quasar host galaxy we find an SFR of >1000M yr −1 based on Herschel and ALMA observations (Paper I; Banerji et al. 2014). The well-established star-forming "main-sequence" of galaxies at high-redshifts (Daddi et al. 2007) can be used to predict how much star formation a typical star-forming galaxy of a certain stellar mass is expected to have at these redshifts. We find that a M * ∼10 11 M galaxy at ∼ 2 would have a star formation rate of ∼200M yr −1 were it to lie on the main-sequence and the quasar host galaxy is therefore clearly experiencing a prolific starburst. Both starburst galaxies as well as the most massive main-sequence galaxies in the high redshift Universe have molecular gas fractions of ∼30-50% Tacconi  If we consider the inner region of the J1234 quasar, the gas fraction is very high and comparable to what is seen in the most gas-rich starburst galaxies. However, over more extended scales the galaxy would appear to be substantially depleted of gas. In the case of the quasar host galaxy, it is plausible that feedback from the AGN has begun to deplete the molecular gas reservoir e.g. as suggested by Brusa et al. (2017). However the two companion galaxies also exhibit somewhat depleted gas reservoirs with currently no evidence for an actively accreting AGN within them. Tadaki et al. (2019) suggest that the most massive galaxies with masses of > 10 11 M in the highest density environments at these redshifts could have gas accretion suppressed due to inefficient cooling of massive haloes, which could explain the fact that both G1234N and G1234S have gas fractions that are at the lower end of what is seen in typical starforming galaxies at these epochs. We caution however that there are considerable uncertainties in both the dynamical and gas masses as illustrated by the discrepancies in these estimates using different dynamical modelling methods and different gas mass tracers.

CO Line Profiles
In Fig. 7 we compare the line profiles for the various different CO transitions presented here and in Paper I for all three galaxies in the J1234 system. Progressively higher J transitions of CO trace warmer, denser gas and a comparison of the line profiles therefore helps illustrate any potential differences in the kinematics of these different gas components. For the quasar J1234 and G1234N all the CO lines appear to be well matched in terms of their velocity and line shape suggesting that the cooler lower density, and warmer, denser gas are experiencing similar large scale motions in these systems. A double-peaked line profile is evident in the CO(7-6) emission from G1234S, which matches what is seen in the CO(3-2) emission. However for this particular galaxy we also see a very broad wing in the CO(7-6) emission extending out to velocities of >1000 km/s. Note -the blue side of the broad-wing is not traced by our observations as the CI(2-1) emission line affects this portion of the spectrum. Curiously, no such broad wings are seen in the CO(3-2) line profile of this galaxy. If the broad wing is ascribed to a molecular outflow, it appears to only be affecting the dense gas in the galaxy. Molecular outflows with similar observational properties to what is seen for G1234S, have previously been observed in luminous, highredshift quasars (Maiolino et al. 2012;Cicone et al. 2015), where the AGN is responsible for driving the gas to velocities of ∼1000 kms −1 . Currently we have no evidence for an AGN being present in G1234S. However, no optical/infra-red counterpart to the galaxy exists in our imaging data of this field down to -band magnitudes of < 25.1 (Banerji et al. 2014), which suggests that the galaxy is highly obscured. We therefore speculate that G1234S with its broad wings seen in CO(7-6) emission and high excitation CO SLED (Paper I) may well host an obscured AGN although high-resolution, deep Xray or mid infra-red data would be needed to confirm or refute this claim. It is also curious that only the dense gas is being affected by any putative outflow in this galaxy. Such an observation would be consistent with recent claims that star-formation can happen in-situ within galactic outflows due to the dense outflowing gas cooling and fragmenting. Indeed recent observations of the prototypical local obscured quasar Mrk231 have shown that the fast outflowing gas in this galaxy shows an unusually high fraction of dense gas, which is conducive to forming stars (Cicone et al. 2020)

Dynamical Mass & Gas Fraction
Our compact configuration ALMA observations for J2315 in B17 had already resolved the CO(3-2) emission in this quasar over physical scales of >20 kpc. In those observations we found that the velocity field could be satisfactorily modelled as a single (albeit very large) rotating gas disk. As can be seen in Fig. 2, the source breaks up into two distinct components in the higher spatial resolution observations presented in this paper, which clearly highlights that many high-redshift, massive galaxies that only have single component gas emission seen in low angular resolution observations, could well be merging systems. The two components are spatially separated by 15.0±0.5 kpc and very close together in velocity, with the northern companion galaxy blueshifted relative to the quasar host by only 170±20 km/s. We fit two 2-D Gaussians to the zeroth moment CO image shown in Fig. 2, finding major axis FWHM of (1.0±0.2) and (1.4±0.2) for the quasar host galaxy and the northern companion respectively, after deconvolving from the beam. This corresponds to half-light radii of ∼4-6 kpc at the redshifts of these galaxies, which is more than 2× the source sizes inferred for all three galaxies in the J1234 system based on the 2D Gaussian fitting (Table 3).
As was done for the galaxies in the J1234 system, we construct intensity, velocity and dispersion maps for both galaxies using the CASA task immoments. These maps can be seen in Fig. 8. The two galaxies of the merger are clearly resolved in the maps with comparable CO line intensities in each merging component. No velocity structure is evident in the individual galaxies at the current resolution, although we clearly see that the northern companion galaxy is blueshifted relative to the quasar host.
For completeness, we also show the best-fit 3D Barolo models fit to the J2315 system as well as the fit residuals in the middle and right-most panels of Fig. 8. The fits are poor as expected as the system cannot be well modelled as a rotating disk and we do not consider the 3D Barolo fitting results in detail in the paper as we do not expect them to be robust. A simple estimate of the total mass of the system can instead be obtained using the virial theorem assuming a projected distance between the galaxies of 15 kpc and a velocity separation of 170 km/s. This leads to M tot =1.0×10 11 M for both galaxies assuming that the galaxies are moving along the line-of-sight.
Alternatively, we can assume that both the galaxies in the J2315 system are individual rotating disks and use the FWHM of the CO line profiles shown in Fig. 2 together with the rotating disk estimator (Neri et al. 2003) to calculate dynamical masses for each galaxy. Inclination angles of 52 degrees and 61 degrees are estimated for the quasar host galaxy and the northern companion respectively by considering the inverse cosine of the minor-to-major axis ratio derived from the CO image. Correcting for these inclination values we find that the quasar host galaxy has M dyn =2.3×10 10 M and the northern companion has M dyn =1.0×10 10 M , corresponding to a total mass of 3.3×10 10 M . Thus the two galaxies appear to be undergoing a major merger with a ∼1:2 mass ratio.
In Paper I we used various different estimators to derive gas masses for this system. The total gas mass is ∼3-4×10 10 M . Combining this with the dynamical mass estimates above and estimating the stellar mass as M * =M dyn −M gas −M BH as was done in the case of J1234, we find very high gas fractions of >45%. The lower limit assumes that the merging galaxies are moving along the line of sight and the gas fraction could be substantially higher if the inclination angle is sufficiently different from face-on.

CO Line Profiles
In Fig. 9 we compare the CO line profiles for the various different CO rotational transitions detected in the J2315 galaxies. The CO(1-0) and CO(7-6) profiles are integrated over the entire merging system as the individual components of the merger are not resolved in these data. We do however show separate CO(3-2) profiles extracted for the quasar host galaxy and the northern companion galaxy. These profiles are considerably less well matched compared to the J1234 quasar host galaxy (Fig. 7) and the kinematics of the gas in this major-merger system is clearly messier as already indicated by the poor 3D Barolo fits. The CO(7-6) emission is significantly broader than the CO(3-2) emission with a FWHM of twice that of the CO(3-2) line (see Table 3 in Paper I). However, as can be seen in Fig. 9, the broadening almost certainly arises because there isn't sufficient angular resolution in the ALMA Band 6 data to disentangle the companion galaxy from the quasar host. The separate CO(3-2) line profiles extracted for the quasar and the companion galaxy, straddle the broad CO(7-6) emission seen in the total system. As already noted in Paper I, the CO(1-0) emission is significantly redshifted by ∼ 500 kms −1 relative to the other CO lines. Thus we have clear evidence for different gas dynamics affecting the cold/diffuse and warm/dense gas components in this major merger. In Paper I we also saw tentative evidence for the CO(1-0) emission originating from closer to the northern companion galaxy rather than the quasar, which could explain the velocity differences observed. However, it is interesting to note that the CO(3-2) emission in the companion galaxy is actually blueshifted rather than redshifted relative to the CO(3-2) in the quasar host. Higher angular resolution observations of all the CO lines are necessary to be able to locate co-spatial components of the gas emission and therefore disentangle the origin of the complicated velocity shifts observed in Fig. 9.

MULTI-WAVELENGTH PROPERTIES OF THE HRQS
In this Section we briefly review existing observational constraints on the multi-wavelength properties of the two HRQs and combine the information provided by these multi-wavelength data with the new constraints on the properties of the quasars from this work, in order to place the quasars and their host galaxies in context within an overall picture of galaxy black-hole co-evolution. Both sources were selected as red, infra-red luminous quasar candidates with broad emission lines present in the near infra-red spectra (Banerji et al. 2012(Banerji et al. , 2015. X-ray observations (Lansbury et al. 2020) place constraints on the accretion properties of the black hole and the line-of-sight gas column densities. Both quasars have similarly high X-ray luminosities -log 10 (L 2−10keV )> 45 -and Eddington ratios of Figure 8. Line intensity (top), velocity (middle) and line dispersion (bottom) maps for ULASJ2315. The source is clearly resolved into a major-merger of two gas-rich galaxies with comparable line intensities. The northern companion galaxy is blueshifted relative to the quasar host. The 3D Barolo best-fit dynamical models and residuals shown for completeness in the middle and right-most panels, do not provide a good description of these data as expected for a merging system.
>0.1, indicating that we are seeing the black holes during a phase of rapid accretion. The gas column densities are also significant -∼10 22 cm −2 , consistent with the quasars being in an obscured phase. In terms of the properties of the central engine, both HRQs are therefore very similar.
Deep optical photometry from current wide-field imaging surveys allows us to trace the rest-frame ultraviolet emission from the HRQs. As the quasar is highly obscured, in the majority of HRQs the rest-frame UV emission likely arises from unobscured star-forming regions in the quasar host galaxy (Wethers et al. 2018). J1234 and its companion galaxies are undetected in the -band down to a 5 magnitude limit of ∼25.1 (Banerji et al. 2014). Far infrared to millimeter observations (Banerji et al. 2014(Banerji et al. , 2018 however suggest significant levels of obscured star formation in the host galaxy. J2315 is detected in data from the Dark Energy Survey with a measured -band magnitude of 22.5 (Wethers et al. 2018).
The source is morphologically extended in the Dark Energy Survey imaging confirming that the UV emission arises from unobscured regions within the quasar host galaxy rather than the quasar itself. The UV-derived star formation rate for the galaxy is estimated to be ∼100-200 M yr −1 (Wethers et al. 2018). Recent, deeper and higher spatial resolution observations from the HyperSuprimeCam Survey Data Release 2 (Aihara et al. 2019) clearly shows resolved emission spatially co-incident with the location of the quasar ( Fig. 10; left panel) as well as more extended emission in the direction of the northern companion galaxy. The rest-frame UV observations of the quasar host galaxies therefore indicate very different levels of obscuration towards the star-forming regions despite similar levels of obscuration towards the black holes, with the host galaxy of J1234 being completely obscured in the UV whereas J2315 shows evidence for unobscured star-forming regions, albeit the main bulk of Figure 9. Comparison of the CO line profiles in J2315. The CO(1-0) and CO(7-6) line profiles are derived over the entire merging system as the individual components of the merger are not resolved in these data. Separate CO(3-2) line profiles for the QSO host and the northern companion galaxy are however shown. The CO(1-0) line profile has been scaled up by a factor of 5 for ease of visualisation. star formation activity in this galaxy is also likely obscured (Banerji et al. 2018;Wethers et al. 2020).
In Paper I we presented constraints on the interstellar medium properties of the two HRQs using detections of multiple molecular emission lines of CO and CI. The excitation conditions in J1234 are consistent with those seen in other luminous AGN with the CO spectral line energy distribution (SLED) peaking at high J values. Photo-dissociation region modelling suggests high radiation field strengths in the interstellar medium of the galaxy. J2315 on the other hand appears to be dominated by lower excitation gas. The CO SLED peaks around J upper =3 and the radiation field strength is ∼1.5-2 dex lower than that inferred for J1234. The ISM conditions in J2315 are therefore more reminiscent of star-forming galaxies rather than AGN. Temple et al. (2019) study narrow-line region outflows in HRQs and find that [OIII] is undetected in J1234, which, even after correcting for the significant levels of obscuration in the galaxy, suggests that the narrow-line region has either been over-ionised or physically removed from the galaxy. This is also the case for ∼10 per cent of very high-luminosity optically selected unobscured quasars ) and consistent with AGN feedback having had a significant effect on the narrow line region gas. In the case of J2315, broad blue wings are seen in the [OIII] emission line profile extending to velocities of ∼2000 km s −1 , providing direct evidence for high-velocity outflowing gas being powered by the luminous AGN. The [OIII] and H region from a recent XShooter spectrum of J2315 are shown in the right-hand panel of Fig. 10 (paper in  preparation).
While we should be cautious about drawing conclusions from a sample of only two quasars studied here in detail, the comprehensive multi-wavelength observations that have been assembled for the two HRQs paint a consistent picture of them being seen in different phases during their evolution. J2315 is conceivably an earlier phase, where two galaxies are viewed in the process of undergoing a major merger. While the quasar itself and a significant amount of the ongoing star formation is obscured by dust, regions of unobscured star formation remain visible in the rest-frame UV with their location broadly consistent with the extended reservoir of molecular gas. Evidence for quasar-driven outflows is provided by the broad blue wings seen in the [OIII] emission line profiles. Nevertheless, the molecular gas fraction is high and the excitation conditions of the interstellar medium are more consistent with those of star-forming galaxies rather than luminous quasars, which might suggest that the effects of AGN feedback are yet to have a dramatic impact on the properties of the host galaxy.
J1234 can putatively be assumed to be seen at a more evolved stage. The quasar together with both companion galaxies are completely obscured at rest-frame UV and optical wavelengths. A significant starburst appears to be occurring in the quasar host galaxy based on the ALMA dust continuum detections (Paper I), which may or may not have been triggered by a merger. While the velocity fields of the CO(3-2) molecular gas observed in this work are consistent with rotational dynamics, the relatively modest resolution of the observations do allow the possibility of a merging system seen close to or soon after final coalescence when the gas is perhaps settling into a rotation pattern. Unlike in the case of J2315, the star formation is entirely obscured in the rest-frame UV. No [OIII] emission is detected in the galaxy, suggesting that powerful AGN driven winds have already over-ionised or cleared the galaxy of its [OIII] gas. Further evidence for efficient feedback is provided by the low molecular gas fraction observed in this system. The ISM properties are very similar to those observed in other high-luminosity quasars with highly excited warm and dense phases of the molecular gas being dominant.

THE M BH -M GAL RELATION
Our spatially resolved observations of the HRQ host galaxies have allowed us to estimate their dynamical masses. Such dynamical mass estimates allow luminous quasars to be placed on the M BH -M gal relation (Willott et al. 2015(Willott et al. , 2017Shao et al. 2017). The majority of optically selected quasars with M BH 10 9 M have black holes that are significantly over-massive relative to their host galaxies, and compared to the local relation (Kormendy & Ho 2013). While this is largely a selection effect of optical surveys targeting the most massive black holes at the rapidly declining end of the mass function (e.g. Lauer et al. 2007), it is also symptomatic of a much larger scatter in the M BH -M gal relation at high-redshifts relative to the local Universe. On the other hand, the population of SCUBA selected submillimeter galaxies at ∼ 2 were inferred by Alexander et al. (2008) to have black holes that are under-massive relative to the local relation. Dust-obscured quasars are a plausible transition population between the starburst and optical quasar phases. The emergence of wide area surveys at mid and far-infrared wavelengths such as the WISE All-Sky Survey and the Herschel Extragalactic Surveys in the last decade, has allowed a wider population of obscured galaxies and AGN to be assembled, with higher dust temperatures than the initial SCUBA detected galaxies. At least some of these also have black holes that are over-massive relative to their hosts (e.g. Assef et al. 2015;Matsuoka et al. 2018) although the number of obscured galaxies with both host and black hole mass measurements remains limited. The location of our heavily reddened quasars on the M BH -M gal plane could therefore place interesting constraints on the relative timescales for black-hole and galaxy growth as these massive galaxies are being assembled.
Before turning to the M BH -M gal ratio in our quasars we first re- evaluate the black-hole mass estimates for our two quasars that were presented in Banerji et al. (2012) and Banerji et al. (2015). As noted in those papers, there are large systematic uncertainties associated with the black-hole mass derivations, not least because the use of a continuum luminosity at 5100Å in the black-hole mass estimate requires significant dust-corrections given the typical extinction levels in the HRQs. These corrections are particularly large for the ALMA sample, which correspond to the reddest quasars in the population. Recent analysis of deep optical photometric data for these quasars also suggests that there could be a significant contribution from the host galaxy to the rest-frame 5100Å continuum fluxes for the reddest quasars (Wethers et al. 2018). We therefore take this opportunity to derive new, more conservative estimates of the black-hole masses for the two quasars in this paper. In particular for ULASJ2315 we utilise a higher signal-to-noise, more red-sensitive XShooter spectrum (paper in prep) compared to the SINFONI spectrum presented in Banerji et al. (2015) which allows a more robust estimate of the H linewidth. We also make use of X-ray observations of both our quasars from Banerji et al. (2014) and Lansbury et al. (2020) to derive hard X-ray luminosities, which suffer significantly less from the effects of extinction compared to the optical continuum. These X-ray luminosities are L 2−10keV = 1.1 × 10 45 erg/s and 3.0 × 10 45 erg/s for J1234 and J2315 respectively. We then use the calibration in Bongiorno et al. (2014) to derive the new black-hole masses: log 10 (M BH ) = 7.11 + 2.06log 10 FWHM H 1000kms −1 + 0.693log 10 L 2−10keV 10 44 ergs −1 The resulting black-hole masses for J1234 and J2315 are 4.1×10 9 M and 2.8×10 9 M respectively, which can be consid-ered to be conservative estimates and are lower than the original black-hole masses presented in Banerji et al. (2012) and Banerji et al. (2015) for these quasars, although still within the ∼0.5 dex systematic uncertainty expected in these black hole mass measurements (Shen et al. 2011). Using these new black-hole masses we show the location of our two ∼ 2.5 dusty quasars on the M BH -M gal relation in Fig. 11 and compare to a higher redshift quasar population where host galaxy dynamical masses have been estimated in a similar way using resolved ALMA observations (Pensabene et al. 2020). We use the dynamical masses inferred from the exponential disk modelling as a proxy for the galaxy mass, as this is more representative of the total mass of the galaxy rather than the mass of the CO bright region (See Table 3). As a reminder, the exponential disk modelling evaluates the mass enclosed within a radius of 10 kpc by assuming a surface brightness profile that extends beyond the CO-bright regions traced by the observations. The Pensabene et al. (2020) high-z sample, which we use as our comparison sample has had galaxy dynamical masses inferred in exactly the same way as our sample -i.e. using an exponential disk model to model the dynamics of the spatially resolved [CII] gas in these quasars and evaluating the mass enclosed within a radius of 10 kpc. While inferences regarding the locations of different populations of quasars on the M BH -M gal plane are complicated by the large systematic uncertainties on single-epoch black-hole mass measurements as well as galaxy mass estimates, conclusions from the comparison between our sample and the z∼ 4 − 6 sample should nevertheless be robust given that both black-hole masses and galaxy masses are evaluated in a very similar way for both samples and make the same underlying assumptions. Our dusty quasars are similar to the luminous ∼ 4 − 6 quasars in that their black-hole masses put them on or above the local relation Thus, even if these sources are seen during a transition phase, their black-holes appear to already be fully grown. Figure 11. Black hole mass versus galaxy mass. In the case of observational data, the galaxy masses all refer to dynamical masses inferred using resolved observations of molecular emission lines arising from the ISM in the quasar host galaxy. In the case of our HRQs the uncertainties on the black hole masses are shown to be 0.5 dex, which encodes the level of systematic uncertainty expected in these masses due to various factors. The resolved dynamical mass estimates for the z∼4-6 quasars are taken from Pensabene et al. (2020) and the black hole masses come from De Rosa et al. (2011Rosa et al. ( , 2014; Trakhtenbrot et al. (2011);Venemans et al. (2015). The grey points show the total baryonic mass contained within 2× the stellar half mass radius for Illustris galaxies (Sĳacki et al. 2015) at ∼ 2. The solid line denotes the Kormendy & Ho (2013) relation.
In Fig. 11 we also show the total baryonic mass contained within 2× the stellar half mass radius, versus black hole mass for ∼ 2 galaxies in the Illustris simulations (Sĳacki et al. 2015), which is described in more detail in the next Section. It is interesting to note that J1234, which is the more "evolved" of our two systems appears to be broadly consistent with an extrapolation of the Illustris galaxies to the highest masses, whereas J2315, the major merger, lies further away from the relation suggesting some further evolution of the total galaxy mass is necessary in this system in order to bring it on to the M BH -M gal relation predicted at these redshifts.

Comparison to Illustris
It is interesting to use the outputs of the Illustris simulations (Vogelsberger et al. 2014;Nelson et al. 2015;Sĳacki et al. 2015) to consider both the progenitors of our ∼ 2.5 quasars as well as the kinds of systems they would evolve into at the present day. The Illustris-1 simulations used here simulate a total comoving volume of (106.5 cMpc) 3 and self-consistently evolve five different types of resolution elements from = 127 to the present day. The public data release paper  contains more details regarding the simulation outputs that we make use of. At a ∼ 2.5, corresponding to the redshift of our quasars, the simulated volume contains 22079 black holes, 345 of which have M BH > 10 8 M and 29 of which have M BH > 10 9 M . While the simulation volume is not large enough to contain > 10 8 M black holes at the highest redshifts of ∼ 6, the simulations can nevertheless be used to track back the ∼ 2.5 population to understand the progenitors of these massive black holes at higher redshifts. DeGraf et al. (2015) and Ricarte & Natarajan (2017) both found that the most massive black holes with M BH ∼ 10 9 M assembled a larger fraction of their black hole mass at earlier times. In the left panel of Fig. 12 we show the evolution of the most massive black holes in Illustris at ∼ 2.5 tracked forward to ∼ 0.5. This confirms that the bulk of the black hole mass is already assembled at ∼ 2.5 and the evolution happens in a mostly horizontal direction from this point onwards with the galaxy stellar mass growing significantly more than the black hole mass. To quantify the direction of the evolution, we define the mean slope, of the evolutionary tracks in the M BH -M gal relation as: = log 10 (M BH,z2 ) − log 10 (M BH,z1 ) log 10 (M * ,z2 ) − log 10 (M * ,z1 ) Selecting the 10 most massive black holes in Illustris at 1 = 2.5 and tracking them forward to 2 = 0.5 results in a slope of 0.64 ± 0.28. This is also broadly consistent with the putative direction of evolution from the J2315 major merger system to the J1234 system discussed in Section 7 and seen in Fig. 11. Given the location of the ∼ 4 − 6 quasars in Fig. 11 relative to our ∼ 2.5 quasars, it is also interesting to ask whether the evolution from ∼ 6 to ∼ 2.5 also occurs in a predominantly horizontal direction in the M BH -M gal plane. As can be seen in the right-hand panel of Fig. 12, this is not the case. The Illustris volume is too small to contain a typical high-luminosity ∼ 6 quasar, but the most massive black holes in Illustris at ∼ 6 with M BH ∼ 10 7 M could putatively evolve to M BH ∼ 10 9 M quasars at ∼ 2.5. The evolution at these epochs happens in a diagonal direction -i.e. between ∼ 6 and ∼ 2.5 black holes tend to grow more rapidly than their host galaxies. Selecting the 10 most massive black holes at 1 = 6 and tracking them forward to 2 = 2.5 results in a mean slope in the M BH -M gal plane of 1.88 ± 0.4, considerably steeper than the slope of the evolution from = 2.5 to = 0.5. The highest redshift quasars shown in Fig. 11 are therefore unlikely to be the progenitors of our ∼ 2.5 quasar population. Indeed, in Illustris the z∼6 progenitors of the most massive z∼2.5 black holes lie ∼ 0.13 ± 0.18 dex below the best-fitting scaling relation, suggesting that although the z∼6 quasars will continue to grow along with their host galaxies, the z∼2.5 quasars likely originated from z∼6 black holes which were under-massive relative to their host galaxies.

CONCLUSIONS
We have presented spatially resolved observations of the CO(3-2) emission from ALMA in two high-luminosity, heavily reddened quasar systems seen at the peak epoch of galaxy formation at ∼ 2.5. Both quasars have very similar luminosities, black-hole masses and dust obscurations and the new ALMA observations clearly demonstrate that there is a range in gas morphologies, gas fractions and dynamical masses seen in these systems.
• The quasar J1234 has two massive, CO-luminous companion galaxies. The CO(3-2) emission shows a clear velocity gradient in the quasar host as well as the two companion galaxies located ∼ 170 kpc and ∼ 90 kpc (projected distance) from the quasar. By modelling these velocity gradients using the tilted-ring model 3D Barolo, we infer dynamical masses of ∼10 11 M in all three galaxies. Dynamical masses estimated using an exponential disk model are larger in all cases, presumably because they trace mass over a much larger physical extent of ∼10 kpc. • The gas fractions in all three galaxies in the J1234 system are somewhat lower than expected from typical star-forming galaxies at these redshifts, particularly when deriving these gas fractions using the full dynamical mass encompassed within ∼ 10 kpc. This is consistent with the process of AGN feedback and inefficient cooling in the most massive halos suppressing the accretion of gas into the system.
• All three galaxies in the J1234 system are at the more massive and more compact end of the distribution of both field and cluster galaxies at ∼ 2 on a mass-size relation.
• We compare the velocity profiles of the CO(3-2) emission to the line profiles of other CO lines detected in the J1234 system galaxies. For the galaxy G1234S, we uncover very broad wings in the CO(7-6) gas emission that are not seen in the CO(3-2) line profiles. This could indicate the presence of a molecular outflow with a much higher fraction of warm, dense gas in the outflow. The implied outflow velocities are high and lead us to speculate about the presence of a "hidden", highly obscured AGN within G1234S, although this needs to be confirmed with deep, high resolution X-ray and/or mid infra-red data.
• We resolve the quasar J2315 into a close separation major merger with a ∼1:2 mass ratio. The two galaxies in the merger have a physical separation of ∼15 kpc and a velocity separation of ∼170 km/s providing evidence that major mergers are associated with luminous accretion and prolific star formation in high-redshift, obscured quasar systems. We also present new rest-frame UV and optical images of the system from HyperSuprimeCam where we see unobscured UV emitting regions spatially co-incident with the extended molecular gas reservoir.
• The gas fraction in J2315 is 45%, which is consistent with or higher than other star-forming galaxies and much higher than in the J1234 system. This could suggest that AGN feedback is yet to fully kick in and that this system is viewed at an earlier evolutionary stage compared to J1234.
• The various CO line profiles in the J2315 merging system show significant velocity shifts relative to one another suggesting that the warm/dense and cool/diffuse gas phases are undergoing different motions and are likely not co-spatial. In particular, the CO(7-6) emission is much broader than the CO(3-2) emission most likely as the two merging galaxies are not resolved in these data, and the CO(1-0) emission is redshifted by ∼500 kms −1 relative to the other CO lines.
• We use the dynamical mass estimates from our resolved observations to place our heavily reddened quasars on the M BH -M gal relation. J1234 is consistent with the predictions from the Illustris simulations for ∼ 2 galaxies. J2315 on the other hand, has a black hole mass that is larger than would be implied by the host galaxy mass. This is consistent with J1234 representing a later phase in the evolution of these massive galaxies compared to J2315.
• We use the outputs from the Illustris simulations to understand both the progenitors and likely descendants of our ∼ 2.5 quasars. We find that most of the mass in the black hole for such massive quasars is expected to be assembled by ∼ 2.5 consistent with very massive black holes already being in place at these redshifts. From this epoch the galaxies mainly grow in stellar mass until the present day. We also demonstrate that the high luminosity ∼ 6 quasar population is unlikely to evolve into our ∼ 2.5 quasars. Instead our massive ∼ 2.5 quasars likely evolved from black-holes that were under-massive relative to their host galaxies at ∼ 6.
Overall these observations highlight the diversity in gas morphologies, environments and gas fractions even within a uniformly selected quasar sample with very similar black-hole properties and AGN obscuration. The differences could be consistent with the two quasars being seen in different evolutionary stages along a majormerger sequence. During final approach, when the galaxies are separated by a few tens of kpc (as is the case for the J2315 major-merger system), the gas fractions are still elevated as the effects of AGN feedback are yet to fully kick in. These feedback processes likely rapidly become more efficient following final coalescence, when the AGN begins to deplete the gas supply leading to the low gas fractions we observe in J1234. It is also important to point out however that there are numerous examples of unobscured, UV-luminous quasars with similar host galaxy properties to both J1234 (e.g. Coppin et al. 2008) and J2315 (e.g. Wagg et al. 2014) suggesting that AGN continuum obscuration is not necessarily a pre-requisite for gas-rich, starburst/merger host galaxies. With a sample of only two obscured quasars analysed in the present work, it is difficult to draw any firm conclusions regarding how they fit into a standard, merger-driven paradigm of massive galaxy formation. Crucially, matched control samples of unobscured, UV-luminous, high-redshift quasars with similar resolved observations of the gas emission are still missing. Assembling these large samples is costly, but nevertheless essential if we are to put the obscured quasar population in context within galaxy formation scenarios.