A Spatially-Resolved Survey of Distant Quasar Host Galaxies: II. Photoionization and Kinematics of the ISM

We present detailed observations of photoionization conditions and galaxy kinematics in eleven z$=1.39-2.59$ radio-loud quasar host galaxies. Data was taken with OSIRIS integral field spectrograph (IFS) and the adaptive optics system at the W.M. Keck Observatory that targeted nebular emission lines (H$\beta$,[OIII],H$\alpha$,[NII]) redshifted into the near-infrared (1-2.4 \micron). We detect extended ionized emission on scales ranging from 1-30 kpc photoionized by stars, shocks, and active galactic nuclei (AGN). Spatially resolved emission-line ratios indicate that our systems reside off the star formation and AGN-mixing sequence on the Baldwin, Phillips $\&$ Terlevich (BPT) diagram at low redshift. The dominant cause of the difference between line ratios of low redshift galaxies and our sample is due to lower gas-phase metallicities, which are 2-5$\times$ less compared to galaxies with AGN in the nearby Universe. Using gas velocity dispersion as a proxy to stellar velocity dispersion and dynamical mass measurement through inclined disk modeling we find that the quasar host galaxies are under-massive relative to their central supermassive black hole (SMBH) mass, with all systems residing off the local scaling ($M_{\bullet}-\sigma~$,$M_{\bullet}-M_{*}~$) relationship. These quasar host galaxies require substantial growth, up to an order of magnitude in stellar mass, to grow into present-day massive elliptical galaxies. Combining these results with part I of our sample paper (Vayner et al. 2021) we find evidence for winds capable of causing feedback before the AGN host galaxies land on the local scaling relation between black hole and galaxy stellar mass, and before the enrichment of the ISM to a level observed in local galaxies with AGN.


INTRODUCTION
Today, feedback from supermassive black holes (SMBH) is an integral part of galaxy evolution models. It is commonly used to explain the lack of observed baryons in local massive galaxies (Behroozi et al. 2010), the enrichment of the circumgalactic medium with metals (Prochaska et al. 2014) and the observed local scaling relation between the mass of the galaxy/bulge and the SMBH (Ferrarese & Merritt 2000;Gebhardt et al. 2000;McConnell & Ma 2013).
The latest observational and theoretical results point to a critical question; at what points does the AGN drive an outflow powerful enough to clear the galaxy of its gas into the surrounding CGM? (King & Pounds 2015) According to theoretical work, this typically happens once the galaxy reaches the M • − σ relationship (Zubovas & King 2014). However, there has been growing evidence for galaxies with massive SMBH and powerful outflows that are offset from the local scaling relationship (Vayner et al. 2017). The origin and evolution of the local scaling relationships with redshift have been an active debate topic over the last decade. When are the local scaling relations established? Are the local scaling relationships the end product of galaxy evolution? Meaning, as galaxies form and evolve, do they fall in and out of the relationships due to rapid growth or feedback process? Do galaxies eventually end up on the local scaling relations once the galaxy or SMBH catch up and finish growing (Volonteri 2012)? Alternatively, is there an inherent evolution in the scaling relationship with redshift and a symbiosis between the galaxy and SMBH growth? (i.e., evolution in slope, offset, and scatter). Finally, there is still an open question regarding the role of quasar feedback in establishing the relationship and whether the merging of galaxies can produce the M • − σ relation following the central limit theorem (Jahnke & Macciò 2011). From a sample of AGN in the COSMOS field (Merloni et al. 2010) finds an offset in the local scaling relationship between redshift 0 and 2. These authors use SED decomposition with numerous spectral bands to measure the stellar mass of the AGN host galaxy in the redshift range of 1 < z < 2.2. From a sample of lensed quasars at 1 < z < 4.5 and broadband HST imaging, Peng et al. (2006) finds an offset in the local scaling relationship. While Sun et al. (2015) using multi-band SED fitting of galaxies in the COSMOS field finds that z∼ 0.2 − 2 galaxies are consistent with being on the local scaling relationship. Schramm & Silverman (2013) using HST imaging in the Chandra Extended Deep Field also finds that galaxies at z∼ 0.6 − 1 are also consistent with being on the local scaling relationship. In the nearby Universe, there is tentative evidence that all of the most massive black holes (> 10 9 M ) are systematically more massive relative to their host galaxies (Martín-Navarro et al. 2018). Fields such as COSMOS or the Extended Chandra Deep Field-South are relatively small in the sky; hence, the number of luminous quasars with massive SMBH is small. Studies that explored the evolution of the local scaling relationships have generally focused on lower-mass black holes with masses < 10 9 M . Furthermore, a large fraction of these studies used broadband HST imaging to study the host galaxies of their quasars/AGN. It is often difficult to disentangle the bright AGN emission from the host galaxy at smaller angular separations (<0.5 ). These studies have a limited number of filters to measure the stellar population's age and the mass to light ratio. Alternatively, mm-interferometry observations have become an essential tool in measuring the dynamical masses of quasar host galaxies across different redshift ranges. At the highest redshifts (z> 4), the [CII] 158µm line has been the most commonly used tracer of the dynamics of the ISM. There is growing evidence that the most massive (> 10 9 M ) SMBH in the highest redshift quasars known to date (z> 6) appear to be over massive for the mass of their host galaxies (Wang et al. 2013;Venemans et al. 2016;Decarli et al. 2018), indicating that the most massive SMBHs residing in high redshift quasars grow first relative to their host galaxies. At more intermediate redshifts 1<z<3, some systems also appear to have overly massive SMBH relative to their stellar/dynamical mass (Shields et al. 2006;Trakhtenbrot et al. 2015;Vayner et al. 2017). While a significant fraction of galaxies with lower SMBH < 10 9 M appear closer or within the scatter of the local scaling relations, galaxies with the luminous quasars and massive SMBH appear to be under massive relative to the mass of their SMBH. As outlined by (Lauer et al. 2007;Schulze & Wisotzki 2014), the offset from the local scaling relations for the systems with more massive black holes is biased due to the steep decline in the galaxy mass function at the massive end.
Integral field spectroscopy (IFS) behind adaptive optics is another method with which it is possible to disentangle the bright quasar emission from the extended emission of the host galaxy. A point spread function can be constructed using data channels confined to the broad emission line of the quasar. After the point spread function is normalized, it is subtracted from the rest of the data channels in the cube. This technique was first shown to be able to resolve host galaxies of low redshift (z < 0.2) luminous type-1 quasars in seeing limited observations (Jahnke et al. 2004) and extended Lyα emission around high redshift quasars (Christensen et al. 2006). Later, when the first near-infrared IFS came online along with their own adaptive optics system, this technique was expanded to samples of higher redshift quasars in search for their host galaxies (Inskip et al. 2011;Vayner et al. 2016) and has shown to work on all the 8-10m class near-infrared IFS (e.g., SIN-FONI, NIFS, and OSIRIS). This technique has shown continued success in seeing limited optical IFS data as well (Herenz et al. 2015;Arrigoni Battaia et al. 2019). This PSF subtraction routine provides better contrast at smaller angular separations than HST, with an inner working angle of 0.1-0.2 , compared to ∼ 0.5 for HST (Vayner et al. 2016). Although today's near-infrared IFSs are not sensitive enough to detect the stellar continuum from the quasar/AGN host galaxies, they can still detect extended ionized emission, enabling us to extract the dynamical properties of the galaxy (Inskip et al. 2011;Vayner et al. 2017) and compare systems to the local scaling relation. However, today, the largest fraction of quasar host galaxy masses still come from HST and mm-interferometric observations. Most likely, selection effects play an important role in determining whether there is a systematic offset from the local scaling relations among the different studies.
Besides measuring the host galaxies and SMBH masses, there are vital open questions regarding the gas phase properties. Galaxies exhibit a correlation between the stellar mass and metallicity across a wide redshift range (Erb et al. 2006a;Sanders et al. 2015). It is often difficult to place galaxies with bright AGN on the mass-metallicity relationship due to limited contrast and the fact that the AGN has a strong impact on the ISM's ionization state. What are the metallicities of the gas in quasar hosts? How does the metallicity in quasar host galaxies evolve with redshift? What is the dominant source of ionization in quasar hosts? What are the star formation rates? One of the best ways to measure the ionization properties of the gas in galaxies is through the BPT (Baldwin, Phillips & Terlevich) diagram (Baldwin et al. 1981;Veilleux & Osterbrock 1987). The traditional BPT diagram plots the ratio of log([OIII]/Hβ ) vs. log([NII]/Hα) and contains two clearly defined sequences: the star-forming sequence and the mixing sequence. The star-forming sequence provides information about the metallicity of HII regions, the stellar ionizing radiation field as well as information on the gas condition in star-forming regions. On the other hand, the mixing sequence consists of gas photoionized by hot stars, AGN, and shocks. It can potentially provide information on the hardness of the AGN ionizing radiation and the metallicity of the gas photoionized by the quasar/AGN (Groves et al. 2006) and shocks. Studies of high redshift star-forming galaxies have shown evidence for elevated line ratios relative to low redshift galaxies. At z∼2, the observed elevated line ratios have been attributed to denser ISM conditions (Sanders et al. 2016) and harder ionizing radiation fields at fixed N/O and O/H abundances relative to typical z=0 galaxies (Strom et al. 2017). Evolutionary BPT models by Kewley et al. (2013a) are consistent with these observations. The evolutionary BPT models also provide a prediction on the evolution of the mixing sequence between z=0 and 3. The location of the mixing sequence moves to lower log([NII]/Hα) value at a relatively fixed log([OIII]/Hβ ) value, primarily due to lower on average gas-phase metallicity at higher redshift (Groves et al. 2006;Kewley et al. 2013a). There is tentative evidence that gas photoionized by AGN is consistent with this picture, as there are several galaxies with AGN, which have emission line ratios offset from the local mixing sequence (Juneau et al. 2014;Coil et al. 2015;Strom et al. 2017;Nesvadba et al. 2017b;Law et al. 2018). Given the presence of the AGN, young stars and shocks in quasar host galaxies, it is crucial to spatially resolve the quasar host galaxy to understand the various contributions to gas ionization. In the distant Universe, this generally requires observations with an IFS and adaptive optics. Resolved BPT diagnostics in both nearby and distant AGN/quasar host galaxies have found regions with distinct photoionization mechanisms (Davies et al. 2014;Williams et al. 2017;Vayner et al. 2017). The question remains whether the ISM condition in the most luminous high redshift quasar host galaxies is different from local AGN and where they lie relative to the mass metallicity relationship.
We have begun a survey to study the host galaxies of z= 1.4−2.6 radio-loud quasars, which are likely to evolve into the most massive elliptical galaxies in the nearby Universe. The sample consists of eleven objects, selected to have young-jets with sizes up to 20 kpc in order to study their impact on galactic scales at early times. The observations consist of near-infrared IFS observation behind laser-guide-star adaptive optics (LGS-AO) at the W.M. Keck Observatory with the OSIRIS instrument. The survey aims to understand the gas phase conditions and ionization mechanisms in high redshift quasar host galaxies and search for evidence of quasar feedback and weighing the masses of the quasar hosts. The observations target nebular emission lines (Hβ ,[OIII] ,Hα , [NII] , [SII] ) redshifted in the near-infrared bands (1 − 2.4 µm), at the distance of our sample, the angular resolution of the OSIRIS/LGS-AO mode corresponds to approximately 1.4 kpc in projection. This paper is part two of two papers focusing on understanding the photoionization mechanisms of gas in radio-loud quasar host galaxies and weigh the mass of the galaxy and SMBH to compare them to the local scaling relations. Refer to paper I (Vayner et al. 2021) for details on the sample selection, properties, and data reduction. Details on archival HST imaging data set are presented in §2. Blackhole masses are presented in §3, we describe how we identify spatially-resolved dynamically quiescent regions in each quasar host galaxy in §4.1, resolved BPT diagrams and our interpretation of the line ratios are present in §5, dynamical masses of the quasar host galaxies and their place relative to the local scaling relations is presented in §7 & §6, we discuss our results in broader context of massive galaxy evolution in §8 and present our conclusions in §9. Notes on individual sources are presented in §9. Throughout the paper we assume a Λ-dominated cosmology (Planck Collaboration et al. 2014) with Ω M =0.308, Ω Λ =0.692, and H o =67.8 km s −1 Mpc −1 . All magnitudes are on the AB scale unless otherwise stated.

ARCHIVAL HST IMAGING
The sources within our sample have a rich set of multiwavelength space and ground-based data sets. To assist in our analysis and interpretation of distinct regions within these quasar host galaxies, we utilize high angular resolution images from the Hubble Space Telescope. We download fully-reduced data from the Barbara A. Mikulski Archive for Space Telescopes (MAST). Table  1 list the archival HST observations used in this study.
We construct a model of the PSF using stars in the vicinity of the quasar within the FOV of each instrument. Images centered on each star are extracted in a box region of roughly 5 x 5 . We then subtract the local background for each star and median combine the stellar images into a normalized "master" PSF. This PSF is then re-scaled to the quasar's peak flux and subtracted out at the spatial location of the quasar. In cases where the quasar was saturated, we scale the flux in the diffraction pattern of the PSF.

BLACK HOLE MASS MEASUREMENT
Blackhole masses are calculated using the broad-Hα line width and luminosity using the scaling relation from Greene & Ho (2005) for a single epoch SMBH mass estimate. We describe the details of the nuclear spectrum fitting in Vayner et al. (2021), which comprises of multi Gaussian models with a broad component for the BLR emission, a narrow Gaussian for the narrow-line region, and an intermediate width Gaussian for the case where there is an outflow. We use the flux and width of the broadest Gaussian to compute the black hole mass. For 3C9, 3C298, there are strong telluric/filter transmission issues that prevent accurate measurement of the FWHM for the emission line. For these targets, we use the Mg II single epoch black hole mass estimate from Shen et al. (2011). The black hole masses are provided in Table  2. We assume an uncertainty of 0.4 dex on the SMBH masses.

DISTINCT REGIONS WITHIN EACH QUASAR HOST GALAXY
In this section we outline how we define various regions within the data cube of each individual object.

Spatially-Resolved Dynamically "Quiescent" Regions
In the first survey paper, we outline our methodology for fitting the emission lines in individual spaxels of our data cubes. From these fits, we derive integrated intensity, velocity, and velocity dispersion maps. The errors on the radial velocity and dispersion maps come directly from the Least-squares Gaussian model fit. The flux map's errors come directly from integrating a noise spectrum in quadrature over the same wavelength range where the emission line is integrated. The fits are presented in the appendix of (Vayner et al. 2021). Here we utilize the radial velocity and dispersion maps to select regions with low-velocity dispersion to search for gas in gravitational motion and search for regions where star formation may have recently happened.
We define a dynamically "quiescent" region of our data set that contains gas with a velocity dispersion (V σ ) less than 250 km s −1 . A quiescent region that belongs to the host galaxy of the quasar must have a radial velocity < 400 km s −1 as we expect the maximum rotational velocity for a given host galaxy to be at most 400 km s −1 . The maximum rotational velocity found for the most massive galaxies studied with IFS at z∼2 is about 400 km s −1 (Förster Schreiber et al. 2018). We define gas with V r > |400| km s −1 and V σ < 250 km s −1 belonging to a merger system. A system is defined as a merger if there are components with V r > |400| km s −1 or more than one distinct kinematic component. For example, in the 3C298 system, two galactic disks are found to be offset by less than 400 km s −1 . All radial velocity and velocity dispersion measurements are relative to the redshift of the quasar. The redshifts for the individual quasars are calculated in Vayner et al. (2021) and are taken from the fit to the narrow-line region. For sources with no spatially unresolved narrow emission, we use the redshift of the broad-line region. We label quiescent regions in the following manner: source name + direction We follow these with a one or two-word comment about the region. Examples of description words are clump, diffuse, or tidal feature. Where clump referrers to a typical few kpc in size compact ionized emission typically seen in high redshift star-forming galaxies. Diffuse referrers to gas that has a surface density of less than typical clumpy star-forming regions. A tidal feature refers to ionized gas associated with a tidal tail in a merging system, containing both diffuse and clumpy ionized gas morphology. For each dynamically quiescent region, we construct a 1D spectrum by integrating over its spaxels. We show an example of this for 4C09.17 in Figure 1, spectra of distinct regions for the rest of the sources are presented in the appendix (Figures 11-18). The emission lines in each spectrum are fit with multi-Gaussian profiles. In these plots, we also present the outflow regions from (Vayner et al. 2021), to illustrate the location of dynamically quiescent regions relative to turbulent regions in these quasar hosts. From these fits, we derive integrated intensity and velocity dispersion that are presented in Tables 3 and 5.

Spatially unresolved narrow-line regions
We search for narrow spatially unresolved emission in each object. To do so, we first subtracted a model of the extended emission from our fits to each emission line in individual spaxels. We then perform aperture photometry on the spatially unresolved emission and extract a spectrum. The emission lines are fit with multiple Gaussian profiles. The fluxes of the narrow emission (σ < 250 km s −1 ) lines from unresolved regions are presented in Table 4. For sources where no unresolved narrow emission line is detected, we place a 1 sigma upper limit on the line flux. Based on the average angular resolution of   Figure 1. On the left, we present the spectra of distinct regions and fits to individual emission lines for the 4C09.17 system. On the right, we present the three-color composite where Hα is color-coded to red, [OIII] to green and [NII] to blue. The contour outlines the spatial location of the region. The bar on the right in each stamp represents 1 or approximately 8.6 kpc at the system's redshift.
about 0.1 , the unresolved narrow line emitting regions' sizes are < 1 kpc.

NEBULAR EMISSION LINE DIAGNOSTICS AND SOURCES OF GAS EXCITATION
In this section, we explore the photoionization mechanism in all distinct regions of each quasar host galaxy. The Baldwin, Phillips & Terlevich (BPT) diagram is used to differentiate between different gas photoionization sources (Baldwin et al. 1981). Here, we use the log([OIII]/Hβ ) and log([NII]/Hα) line flux ratios to distinguish heating from young stars, AGN, and shocks.
To construct the BPT diagram for our sources, we integrated each emission line over the same velocity width (∆V) and velocity offset relative to the redshift derived from the [OIII] emission line at each spaxel. We integrate the maps relative to [OIII] since it is typically the brightest emission line in any given spaxel. The higher signal-to-noise [OIII] emission line leads to a smaller spaxel-spaxel variation in the radial velocity and dispersion maps, creating a more consistent log([OIII]/Hβ ) and the log([NII]/Hα) ratio between neighboring spaxels. We find that for the entire sample, the standard deviation on the log([OIII]/Hβ ) ratio decreases by 0.2 dex compared to when integrating the cubes relative to the Hα line.
A resolved BPT diagram allows us to investigate the source of ionization throughout each quasar host galaxy. Due to sensitivity and, in some cases, wavelength coverage, we cannot create an integrated emission-line map for Hβ on a similar scale to Hα , [OIII] , or [NII] maps. For our BPT diagrams, we construct our Hβ map by assuming case B recombination (Hβ =Hα /2.86) with a gas temperature of 10 4 K and an electron density of 10 2 cm −3 . Assuming other recombination cases and ISM conditions with reasonable temperatures and densities would not change our results by a significant amount as the ratios between Hβ and Hα would only change at most by a factor of ∼1.3 (Osterbrock & Ferland 2006). For sources with the brightest extended emission and wavelength coverage of both Hα and Hβ we find a maximum V band extinction of 1 mag, however in most cases, line ratios consistent with case B recombination. In regions where gas extinction is present, the log([OIII]/Hβ ) ratios are preferentially lower.
Only spaxels where at least Hα and [OIII] were detected are analyzed and presented here. Typically [NII] (Vayner et al. 2021) we over plot their line ratios and label them with a star. Individual spaxels typically have high uncertainties in their ratios but tend to cluster together on the BPT diagram. Integrating over distinct regions and re-calculating the ratios from a high SNR spectrum confirms that region's true line ratio.
To conserve space, we do not over-plot the error bars on points from individual spaxels in Figure 2, we only show the error bars of ratios computed for integrated values of the distinct regions. In Figure 3, we plot points of individual spaxels along with the error bars.

Ionization Diagnostic Models
We find that a large portion of our line ratios values lies outside the two typical sequences of the BPT diagram (Figure 2). At a fixed log([NII]/Hα) nearly, all values are above both the local mixing and star-forming sequence. At a fixed log([OIII]/Hβ ) value, nearly all values are outside the local mixing sequence. A large portion of points falls between the star-forming and mixing sequence, with relatively high log([OIII]/Hβ ) values. Metallicity, electron density, the hardness of the ionization field, and the ionization parameter determines the location of the galaxy/region on a given sequence. With changing conditions in the ISM between z=0 and the median redshift of our sample, the locus of both the star-forming and mixing sequence can change locations (Kewley et al. 2013a).
Galaxies at a fixed stellar mass have lower metallicities at high redshift compared to galaxies today (Yuan et al. 2013). Near the peak of the star formation rate density at z∼ 1 − 3, the ISM conditions and star formation histories of star-forming galaxies may differ from local systems. Star formation appears to happen in denser environments in the distant Universe with higher electron densities (Sanders et al. 2016), akin to conditions seen in local ULIRGs. According to Steidel et al. (2014); Strom et al. (2017) ISM in high redshift galaxies experiences harder ionization at a fixed N/O and O/H abundances than z=0 star-forming galaxies. On the other hand, galaxies at higher redshift have elevated N/O rations . Taken together, Kewley et al. (2013a) shows that such changes in ISM conditions can alter the location of the star formation sequence between z=0 and z=2.5. Notably, the combination of harder ionization, electron density, and ionization parameter can shift the locus of the star-forming sequence to higher log([NII]/Hα) and log([OIII]/Hβ ) values. It appears that UV/emission line selected galaxy samples tend to show a more significant shift from the SDSS star formation locus, as evident in a large sample of 377 starforming galaxies explored by Strom et al. (2017). Nearly all their galaxies have a higher log([OIII]/Hβ ) value at a fixed log([NII]/Hα) compared to local galaxies studied in SDSS. Various galaxy selection techniques may lead to samples of galaxies with inherently different ionization properties. However, the overall conclusion from studying star-forming galaxies in the distant Universe is that the line ratios of these systems lie on different star formation locus compared to the local Universe.
Changes in the ISM conditions of distant galaxies may also lead to changes in the location of the mixing sequence. Kewley et al. (2013a) and Groves et al. (2006) show that for galaxies with lower metallicities, the mixing sequences shifts to lower log([NII]/Hα) values with relatively small changes in the log([OIII]/Hβ ) value. We explore the various evolutionary models of the starforming and mixing sequence with redshift and ISM conditions proposed by Kewley et al. (2013a). The best fit model to our sample is the one where the ISM of high redshift galaxies have more extreme condition (higher electron density, harder ionization field, and larger ionization parameters) and the metallicity of the gas photoionized by the quasar is at a lower metallicity compared to the gas ionized by local AGN in the SDSS sample. The median log([NII]/Hα) value is about 1.0 dex lower than that of the mixing sequence at z=0. If the primary source in the shift of the mixing sequence from z=0 to z=1.5-2.5 is a change in the gas phase metallicity, then the gas photoionized by the quasar in our sample has a metallicity a 1/2-1/5 of that in narrow line regions of z=0 AGN on the Kewley & Dopita (2002) metallicity scale. One of the consequences of the shift in the mixing sequence is that it becomes harder to distinguish between gas photoionized by AGN vs. star formation, especially in systems with potentially multiple ionization sources. Changes in the photoionization condition likely also play a role in the offset from the local mixing sequence. In a sample of local type-1 quasars, (Stern & Laor 2013) found that systems with higher bolometric luminosities and higher Eddington ratios are systematically offset to lower log([NII]/Hα) ratios.   Figure 2. Line ratio diagnostics of individual resolved distinct regions. In grey, we plot the line ratios of individual spaxels where at least [OIII] and Hα was detected at an SNR>3. Uncertainties on these line ratios are generally large; hence, we also integrate over all spaxels in individual regions to increase the SNR and lower the uncertainties on the line ratios. We show region-integrated line ratios with the colored symbols where each object has the same symbol, and each region has a different color. The names of the distinct region are present in the lower-left corner, and these match the names given in Table 3. We present the evolutionary models of the mixing and star-forming sequence with red and green curves from Kewley et al. (2013a). We show the upper limit of a sequence with a straight line and the lower boundary of each sequence with a dashed curve. Teal curves represent the bounds of the two sequences where the majority of the line ratio in low redshift galaxies fall. Our line ratios are consistent with a model where gas photoionized by the quasar is denser, has lower metallicity, and experiencing harder ionization compared to the gas photoionized by AGN in nearby galaxies.
Most of the gas in our quasar host galaxies lies in the mixing sequence where the gas is photoionized by a combination of quasar ionization and radiative shocks. In Figure 3, a significant fraction of the points in individual objects match the predicted location of radiative shocks on the BPT diagram. The radiative shock models assume solar metallicity and a preshock density of 215 cm −3 . With the present data, it is difficult to distinguish the percentage of photoionization from shocks vs. AGN. However, given the overlap of both line ratios and kinematics with shock models, we cannot rule them out to contribute to gas photoionization.
A number of our distinct regions appear to have low log([NII]/Hα) values (<0.5) with low-velocity dispersion gas (V σ < 250 km s −1 ). Morphologically these regions appear to be clumpy in their Hα maps reminiscent of typical star-forming regions in galaxies at z > 1. The line ratios of these points do not coincide with regions of fast or slow shocks photoionization on the BPT diagram (Allen et al. 2008;Newman et al. 2014). Archival HST data of 3C9, 3C298, 4C09.17, 3C268.4, 3C446 all showcase that the dynamically "quiescent" regions in these galaxies have clumpy morphology in rest-frame UV continuum data, similar to those of star-forming galaxies at these redshifts. In Figure 4, we overlay the Hα emission from dynamically quiescent regions onto archival HST observations at rest-frame UV wavelength. Combining these clues suggests that the quasar does not entirely photoionize the gas in these regions. The elevated log([OIII]/Hβ ) in these regions compared to local and distant star-forming regions may be from the mixing of photoionization from massive stars and the quasar. There is some evidence for this based on the morphology of the ionized gas and their respective log([OIII]/Hβ ) ratios. For example, in 4C09.17, we see that more diffuse emission with low-velocity dispersion tends to have a higher log([OIII]/Hβ ) value compared to clumpier regions where there is evidence for recent star formation activity and potentially more shielding from the quasar radiation field. Using the empirical star formation rate Hα luminosity relationship from Kennicutt (1998), we convert the Hα luminosity of the distinct extended quiescent regions to star formation rates. Most likely, the majority of these regions are photoionized by a combination of AGN and star formation, hence the derived star formation rates are upper limits. Regions "3C298 SE component B Tidal feature" and "4C09.17 W component B clumps" have line rations most consistent with photoionization by O/B stars, the star formation rates derived in these regions are closer to their actual value. To partially address the contribution from AGN to photoionization in dynamically quiescent regions, we also derive star formation rates only using spaxles that fall within the line ratio error inside the star-forming region of the BPT diagram based on the Kewley et al. (2013b) maximum separation between star formation and AGN photoionization. Generally, we find lower (1/2 -1/10) star formation rates when using the BPT diagram criteria. We also measure the metallicities of these regions using the Pettini & Pagel (2004) empirical gas-phase metallicitylog([NII]/Hα) relationship. Given that log([NII]/Hα) is elevated in the presence of an AGN/quasar ionization field, the metallicities for the majority of the regions are also upper limits. We also calculate the metallicity using theoretically derived chemical abundance calibration for narrow-line regions of AGN (Storchi-Bergmann et al. 1998). We present quantitative values of these regions in Table 5, where we show the Hα luminosity of each quiescent region, along with the star formation rate, metallicities, and velocity dispersion. Since nearly all of the unresolved narrow line regions are consistent with quasar photoionization, we do not include them in Table 5 with the exception of 3C318. In this object, the line ratios are consistent with star formation, indicating a nuclear starburst on scales <1 kpc. For cases where we do not detect any unresolved narrow emission, we place an average upper limit on the star formation rate. We do so if there is an ongoing nuclear starburst with a star formation rate below the sensitivity of OSIRIS at our given exposure times. We find an average star formation rate limit of 9 M yr −1 for the four objects (4C05. 84, 4C04.81, 4C57.29, and 7C1354) where no strong narrow nuclear emission was detected.

SMBH-GALAXY SCALING RELATIONSHIPS
In this section, we place our galaxies on the velocity dispersion and galaxy mass vs. SMBH mass plots, comparing their locations to the local scaling relations (M • −σ and M • −M * ). We calculate the SMBH masses from the broad Hα luminosity and line width using the methodology presented in Greene & Ho (2005). The SMBH masses span a range of 10 8.87−9.87 M . The velocity dispersions are taken from dynamically quiescent regions, while the galaxy masses are calculated from the virial equation and from modeling the radial velocity of targets with rotating disks and extracting a dynamical mass.

Host Galaxy Velocity Dispersion
We identify several dynamically quiescent regions within most of the quasar host galaxies in our sample. These regions show lower log([NII]/Hα) line ratios and typically have clumpy morphology, reminiscent of the general star-forming regions seen in nebular emission and UV continuum in high redshift galaxies. In most galaxies, these regions lie away from any galacticscale outflows. Hence their observed dynamics could be a probe of the galactic gravitational potential. These regions can be used to measure the velocity dispersion of our quasar host galaxies. In combination with the measured black hole masses, we can compare them to the local scaling relation between the SMBH mass and the velocity dispersion of the galaxy/bulge. In Figure  5, we plot the mass of the SMBH presented in Table 2 against the velocity dispersion of distinct quiescent regions measured with the Hα line. Also, we include the velocity dispersion measured from CO (3-2) emission for 3C 298 from Vayner et al. (2017). We find a significant offset from the local scaling relation between the SMBH mass and the velocity dispersion of the galaxy/bulge (M • − σ ) (Gültekin et al. 2009;McConnell & Ma 2013). To address the issue that the velocity dispersion may be systematically lower in dynamically quiescent regions offset from the quasar (3C446) or regions where the surface area of the narrow emission is significantly lower than the outflow (4C09.17, 3C298), we recalculate the velocity dispersion in a larger aperture that includes outflows and narrow emission. We see no strong systematic difference in the velocity dispersion of the narrow gas. The source integrated narrow velocity dispersion for 3C298, 3C446 and 4C09.17 are 100.7 ± 16.6, 187.5 ± 1.0 and 144.0±10.0 km s −1 , respectively. In the nearby universe, the velocity dispersion is typically measured inside the effective radius of the galactic bulge. The difference within our observations is that we do not know the bulges' true sizes for our galaxies as we have no way to measure them with current data and telescope/instrument technology. However, the extents of the dynamically quiescent regions are in the range of the effective radii for bulges in massive (10 10.5−11.5 M ) galaxies studied in the CANDELS survey (Bruce et al. 2014).
There have seen numerous discussions in the literature about whether the velocity dispersion measured from gas traces the stellar velocity dispersion. The gas and stars might not have the same uniform distribution, and winds can broaden the nebular emission lines. Furthermore, the line of sight absorption and emission lines from which the velocity dispersion is calculated are luminosity weighted subject to galactic dust extinction. Because of the different light distribution between stars and gas, the measured velocity dispersion can differ. These differences can lead to increased scattering in any correlation between σ * and σ gas . Data-sets that spatially resolve the gas and stellar components and have enough resolving power to separate multi-component emission from different regions (e.g., outflowing/in-flowing gas vs. galactic disk rotation) are important when making a comparison between σ * and σ gas . In Bennert et al. (2018), for a large sample of local AGN, they find when fitting a single Gaussian component to the [OIII] emission line, they can overestimate the stellar velocity dispersion by about 50-100%. Only by fitting multiple Gaussian components to account for both the narrow core and the broader wings of the [OIII] line profile can they adequately match the velocity dispersion from the narrow component of the [OIII] line to that of the stellar velocity dispersion. For their entire sample, the average ratio between the velocity dispersion of narrow Gaussian component and the stellar velocity dispersion is ∼1. The 1 σ scatter on the ratio between σ [OIII],narrow and σ * is about 0.32 with a maximum measured ratio of about a factor of 2 which translates to a scatter in ∆σ = σ [OIII] − σ * of 43.22 km s −1 with a maximum difference of about ±100 km s −1 . However, only a few sources show such drastic velocity differences (∼ 2.5% of the entire sample, 82% of the sources show |σ [OIII] − σ * | < 50 km s −1 ). When fitting for the M • − σ relationship with the narrow [OIII] emission a Star formation rate derived using the Hα luminosity of the entire distinct quiescent region b Star formation rate derived using spaxels that fall within the star formation sequence on the BPT diagram.
c Value indistinguishable from the integrated value over the entire dynamically quiescent region d Tidal tail as a proxy for stellar velocity dispersion, the resultant fit agrees with that of quiescent galaxies reverberationmapped AGNs. These results indicate that for the sample as a whole Bennert et al. (2018) finds that both the stars and gas follow the same gravitation potential. Given the Bennert et al. (2018) results that demonstrates that the gas velocity dispersion can be used as a proxy for the stellar velocity dispersion, we follow a similar analysis using our IFS data sets to explore the location of our galaxies relative to the M • −σ relation at high redshift. We attempted to the best of our ability to separate regions that contain galactic scales winds from those with more quiescent kinematics both spectrally and spatially with OSIRIS. Hence similar to Bennert et al. (2018) we think that the measured velocity dispersions in quiescent regions are good tracers of the galactic potential on average. Throughout the paper we use the narrow velocity dispersion of [OIII] and Hα emission lines of dynamically quiescent regions as a proxy for the stellar velocity dispersion. We still find a significant offset for our sample after applying the observed scatter in the difference between σ * and σ gas . This is also true when applied to the more distant quasar host galaxies studied with 158 µm [CII] emission. In nearby galaxies, there is a dependence on the velocity dispersion with the radius from the galaxy center (Bennert et al. 2018;Ene et al. 2019). However, based on the local galaxy observations, the velocity dispersion is unlikely to increase by ∼ 200 km s −1 that is necessary to bring the galaxies onto the local scaling relations.
Using N-body smoothed-particle hydrodynamics simulations Stickley & Canalizo (2014) examines how the stellar velocity dispersion evolves in a binary galaxy merger. At various stages in the merger (e.g., a close passage, nucleus coalescence), they measure the stellar velocity dispersion along 10 3 random lines of sight. Near each close passage and during coalescence, they find that the scatter on the velocity dispersion significantly increases from ∼ 5 − 11 km s −1 to about 60 km s −1 with the average velocity dispersion a factor of ∼1.7 higher than after the galaxies have finished merging. For several sources in our sample (3C9, 3C298, and 3C446), the measured velocity dispersion might be higher than what it will be once the galactic merger is complete adding uncertainty due to projection effects. Following the simulations results, we add in quadrature an additional un-

3C9
4C09.17 3C446 3C268.4 3C298 certainty on the velocity dispersion of 60 km s −1 given that the majority of our mergers are near coalescence or a close passage (∆R < 10 kpc).
It should be noted that this is near the maximum scatter seen in the simulations on σ. These simulations also find that for merging galaxies at their maximum separation, the measured velocity could be a factor of ∼ 1.7 times smaller compared to the final system. They find that for a 1:1 merger, the maximum separation after the first passage is 10-100 kpc, which is much larger than any separation that we find in our systems from observed projected separations and measured relative velocities. No obvious merging companions are found for 3C318, 4C22.44, or 4C05.84 hence for these systems, the mergers might be past their coalescence stage where the measured velocity dispersion is close to its final value, and the scatter due to the line of sight effects is minimal (∼ a few km s −1 ). However, we still apply an additional 60 km s −1 uncertainty in these regions.
Even after these corrections are made to approximate the stellar velocity dispersion from the [OIII] emission lines in our sample, we still find that all of our systems are offset from the local scaling relation between the mass of the SMBH and the velocity dispersion of the bulge/galaxy. Given that we are dealing with relatively small sample size, we performed statistical tests to confirm the offset between the local scaling relation and our sample. We measure the offset between the observed and predicted velocity dispersion for the SMBH mass of our systems for each object. We use the local scaling relation fit from (McConnell & Ma 2013), and Hα measured SMBH masses. We construct a data set consisting of velocity differences. From bootstrap re-sampling of the velocity difference data set, we find that the average offset of 188.7 km s −1 is significant at the 3.25σ level. Using Jackknife re-sampling similarly, we find that the offset is significant at the 3.3σ level with the 95% confidence intervals of 154.4 km s −1 to 223.0 km s −1 on the velocity dispersion offset. Performing similar statistical tests on the Decarli et al. (2018) sample, we find an average offset of 178.8 km s −1 with a significance of the shift at 2.7σ and 2.8σ for Jackknife and bootstrap resampling, respectively from the local relationship. We also measure the offsets of massive BCGs in the local Universe from the M • − σ relationship. Using a twosided Kolmogorov-Smirnov test, we can ask if the observed offsets of the local and high redshift data sets are drawn from the same continuous distribution. We find a Figure 5. The location of our galaxies on the velocity dispersion vs. SMBH mass plot compared to the local M• −σ relationship. We use the narrow Hα emission line velocity dispersion of dynamically quiescent regions as a proxy for the stellar velocity dispersion. Red stars are the measurements from our sample, where we measure the velocity dispersion from the narrow Hα emission line. We measure the black hole masses using the broad Hα line from the broad-line-region discussed in section 3. The two blue stars represent the velocity dispersion measured in the disk of the host galaxy of 3C 298 and the tidal tail feature 21 kpc away from the quasar (Vayner et al. 2017). Blue circles are quasars from the Shields et al. (2006) sample, where they measure the velocity dispersion from CO emission lines. The yellow points are from quasars at z> 6, where they measure the velocity dispersion from the 158 µm [CII] emission line (Decarli et al. 2018). Green points represent the local sample of all the bright cluster galaxies with SMBH greater than 10 9 M taken from McConnell & Ma (2013). The blue curve represents the best fit to the entire galaxy sample from McConnell & Ma (2013) with the blue shaded region represents the intrinsic scatter on the M• − σ relationship from the fit while the green curve is the fit to the sample studied in Gültekin et al. (2009). We find a significant offset between the galaxies in our sample and local BCG and the general local M• − σ relationship. This large offset indicates that the host galaxies appear to be under-massive for their SMBHs.
p-value of 5.7×10 −9 , indicating that the two populations are not drawn from the same distribution. Applying the Kolmogorov-Smirnov test to the velocity dispersion offsets from our sample and in the higher redshift quasars, we find a p-value of 0.84, indicating that these two data sets could be drawn from the same continuous distribution. We find similar results by comparing the Shields et al. (2006) sample at z ∼ 2 to our own and that of Decarli et al. (2018).

DYNAMICAL MASS MEASUREMENTS
We can also test whether these systems lie off the local scaling relationship between the SMBH mass and the dynamical mass of the bulge/galaxy. First by using a virial estimator for the dynamical mass of the galaxy M virial = Cσ 2 r G where C=5 for a uniform rotating sphere (Erb et al. 2006b). We assume 7 kpc for the radius, which is the median effective radius of massive quiescent galaxies in the local Universe (Ene et al. 2019). Here σ is derived from a Gaussian fit to the integrated spectra over the distinct region. For galaxies with multiple distinct regions, we derive two or more dynamical masses as there may be a dependence on the velocity dispersion as a function of position with the galaxy. For systems in a clear merger, the galactic component belonging to the quasar is used to estimate the dynamical mass since we are interested in the correlation between the SMBH and the velocity dispersion of the quasar host galaxy.
For systems with velocity shear in the 2D radial velocity map, we fit a 2D inclined disk model to the kine- Figure 6. We present the location of individual galaxies compared to the local scaling relation between the mass of the SMBH and mass of the galaxy/bulge shown with a blue curve. Blue points represent systems with virial dynamical masses. Red points represent systems where we calculate the dynamical mass by modeling the radial velocity maps with an inclined disk model. Gray points show the location of galaxies at z> 0.5, with lower SMBH masses and lower AGN luminosity compared to our sample. The blue curve represents the local scaling relationship as measured in McConnell & Ma (2013), with the shaded region representing the intrinsic scatter. We find the majority of our points are offset from the local scaling relationship, outside the observed scatter. matics data to measure the dynamical mass. The model is a 2D arctan function where V(r) is rotation velocity at radius r from the dynamical center, V max , is the asymptotic velocity, and r dyn is the radius at which the arc-tangent function transitions from increasing to flat velocity. The measured line-of-sight velocity from our observations relates to V(r) as where cos θ = (sin φ(x 0 − x)) + (cos φ(y 0 − y)) r . ( Radial distance from the dynamical center to each spaxel is given by where x 0 , y 0 is spaxel location of the dynamical center, we quote the value relative to the centroid of the quasar, V 0 is velocity offset at the dynamical center relative to the redshift of the quasar, φ is position angle in spaxel space, and i is the inclination of the disk. V max is not the true "plateau" velocity of the galaxy's disk. V max can have arbitrarily large numbers, especially when r dyn is very small (Courteau 1997). To fit the data we use the MCMC code emcee. We construct the model in a grid with a smaller plate scale than the observed data, which gets convolved with a 2D Gaussian PSF with an FWHM measured from the quasar PSF image. The image is then re-sized to the plate scale of the data. We construct the priors on each of the seven free parameters. The prior on V max is 300 < V max < 1000 km s −1 the prior on both x 0 , y 0 is the boundary of the FOV of the imaged area, the prior on the position angle is 0 < φ < 2π, the prior on the inclination angle is 0 < i < π/2, the prior on the radius is 0.5 < r dyn < 10 pixels and the prior on V 0 is −100 < V 0 < 100 km s −1 . We then sample this distribution with emcee. We initialize 1000 walkers for each free parameter using the best fit values from leastsquares fitting as the starting point, with a small random perturbation in each walker. We run MCMC for 500 steps starting from the perturbed initial value. The best-fit parameters, along with their confidence intervals, are presented in 6 for the quasar host galaxies of 7C 1354+2552, 3C9. For 3C 298 we do not see the disk in the ionized emission with the OSIRIS data, it is solely detected in CO (3-2) observations from ALMA, here we present the best fit values from Vayner et al. (2017). Also, we present ∆v obs /2, the average between the maximum and the minimum velocity along the kinematic major axis as determined by the position angle (φ). We also present the intrinsic velocity dispersion (σ 0 ), measured along the kinematic major axis, towards the outskirts, away from the steep velocity gradient near the center of the disk. In addition, we measure V rot /σ 0 to gauge whether these systems are dynamically supported by rotation or dispersion. We measure a value of 6.8±1, 2.7±0.6, and 4.4±1.5 for 7C 1354, 3C9, and 3C298, respectively. In all systems, rotation dominates over the velocity dispersion for the dynamical support according to the criteria outlined by Förster Schreiber et al. (2018), and henceforth the systems can be classified as true disks.
Assuming a spherically symmetric system, we can compute the total enclosed mass using the following formula: M (R) = 2.33 × 10 5 rV 2 r / sin(i) 2 Where V r is the radial velocity, i is the inclination angle from the disk fit. For the radial velocity we use ∆v obs /2. Similarly, we assume a radius that is the median value of nearby BCGs (7.1 kpc). The selected radius should give us an absolute upper limit on the dynamical mass of the galaxy/bulge as this radius is much larger than the typical size of a galactic bulge at this redshift and is larger than the observed extent of the galactic disks. The reason for choosing a larger radius is to address the case where the quasar host galaxy extends to a larger radius and is not captured in our OSIRIS observations because they are not sensitive enough to low surface brightness emission at larger separation from the quasar. Virial and dynamical masses are presented in 7. However, it is not guaranteed that the extent of the ionized gas will match the stellar. We attempted to measure the size of the stellar continuum from the HST observations but were unsuccessful. Using the Galfit package, we were unable to constrain the radius due to the sources' complex morphologies and the increased inner working angle due to the quasars' brightness and saturated counts in the HST observations. Due to the limited sensitivity of OSIRIS to lower surface brightness emission, we are missing an accurate measurement of the plateau velocity for the galactic disks at large separations from the quasar. Hence, our fitting routine is unable to constrain V max for 3C9 and 7C 1354. Also, it appears that the turn over radius is very small for these two systems, smaller than the resolution element of our observations. For this reason, we are unable to constrain the turn over radius, and we only provide a limit. For 7C 1354, there is a degeneracy between the maximum velocity, turn over radius, and inclination; thus, the values that we provide are those that give the smallest velocity residual.
Using the measured virial and disk fit dynamical masses and the SMBH masses, we can now compare our galaxies to the local M • − M * relationship. Not only are these galaxies offset from the local M • − σ relationship, but we also find that these galaxies are on an average offset from the local M • − M * relationship. The galaxies need about an order of magnitude of stellar growth if they are to evolve into the present-day massive elliptical galaxies.
We note that we have used two different methods for exploring the scaling relationship for galaxy mass vs. SMBH. Both the gas velocity dispersion method and dynamical measurement imply that the SMBH is over- massive compared to their host galaxies when exploring the local scaling relationship. It will be important to further compare these methods with larger samples, as well as future observations with the James Webb Space Telescope that will be able to directly measure the stellar velocity dispersion.

DISCUSSION
Our survey aimed to study host galaxies of redshift 1.4 -2.6 radio-loud quasars through rest frame nebular emission lines redshifted into the near-infrared.
We place distinct regions of each quasar host galaxy on the traditional BPT diagram (log([OIII]/Hβ ) vs. log([NII]/Hα) ). The majority of the points for our  Figure 9. Fitting an inclined disk model to the radial velocity map of the 3C298 quasar host galaxy. Far-left we plot the isolated radial velocity structure belonging to the quasar host galaxy of 3C298, middle left shows the best fit model overlaid as contours on top of the radial velocity map, middle right is the best fit model, and on the right, we plot the residuals.
sources lie outside the two local sequences (the mixing and star-forming sequence). In section 5, we introduce evolutionary BPT models from Kewley et al. (2013a) that indicate changes in the photoionization and metallicity conditions of the gas can shift both of the starforming and mixing sequences. We fit these models to our data and find that the best-fitting model is the one where the gas in our quasar host galaxies is at least two to five times less metal-rich compared to the narrow line regions of nearby (z<0.2) AGN. The best fit model also indicates that the gas is ten times denser compared to nearby galaxies. In Figure 2, we show all of our points on the BPT diagram along with the best fit model. Kewley et al. (2013b) studied a sample of star-forming galaxies and galaxies with AGN in the redshift range of 0.8<z<2.5. They also find that galaxies at z>2 show elevated line ratios on average outside the local star formation and mixing sequences. They find that normal ISM conditions similar to the SDSS sample transition to the more extreme conditions with elevated line ratios somewhere between redshift z=1.5 and z=2. This is an agreement with our results as the majority of our targets are at z > 1.5.
High redshift radio galaxies also appear to show ISM conditions with metallicities that are lower compared to local AGN. In a study of a large sample of distant radio galaxies, Nesvadba et al. (2017a) finds that their gas-phase metallicities are at least half of that seen in local AGN. Nesvadba et al. (2017a) finds the same bestfitting model from Kewley et al. (2013a) as we do for our sample to explain their observed nebular line ratios. The average log([NII]/Hα) value of our sample seems to be lower than that of Nesvadba et al. (2017a); this could be due to the lower metallicity of our sample. On the other hand, a different approach to how we compute our line ratios can cause the discrepancy. Nesvadba et al. (2017a) only presents source integrated line ratios, while we explore ratios of distinct regions because we typically have a factor of 5-10 better angular resolution due to adaptive optics and hence can resolve the different ionized/kinematics structures of our galaxies. In the majority of our sources, we see significant variations in log([NII]/Hα) and log([OIII]/Hβ ) values across each system, hence why we explore distinct regions. Line ratios from integrated spectra that include regions with various ionization sources and from multiple components of a merger system may shift towards higher log([NII]/Hα) , and log([OIII]/Hβ ) values as the regions photoionized by the quasar/AGN tend to be brighter. Line ratios of galaxies with lower luminosity AGN compared to quasars/radio galaxies studied in Strom et al. (2017) are nearly all outside the local mixing sequence. These points overlap with the location of our line ratios and that of the radio galaxy sample. The MOSDEF survey finds similar results for their AGN sample at a range of bolometric luminosities (Coil et al. 2015). The ubiquity of elevated line ratios in host galaxies of AGN, meaning they are typically above the local mixing or star-forming sequence on the traditional BPT diagram (log([OIII]/Hβ ) vs. log([NII]/Hα) ), indicates that regardless of the active galaxy population selected at z∼2 the conditions of the gas that is photoionized by an AGN are different from those in the local Universe.
Overall, this suggests that the ISM conditions in high redshift galaxies with AGN at a range of bolometric luminosities are different from those in local systems. The ISM conditions appear to be far more extreme with gasphase metallicity lower than that of local AGN, suggesting an evolution in the ISM gas that is photoionized by AGN from z=0 to z=2.5. 8.1. Star formation and dynamically "quiescent" regions in the host galaxies In 9/11 quasar host galaxies within our sample, we see the morphology of clumpy star-forming regions seen in other galaxies at these redshifts. These regions also typically show lower velocity dispersion and lower log([NII]/Hα) values. We described them in more detail in section 4.1. These regions lie 1 -21 kpc from the quasar and generally do not coincide with the location of galactic outflows. For sources with available HST imaging of rest-frame UV continuum, these regions appear bright and clumpy (see Figure 4). Taking these two results together indicates that O and B stars could photoionize a non-negligible fraction of the gas in these clumpy regions. In section 5, we derive an upper limit on their star formation rates and gas-phase metallicities.
Taking this together, there is evidence for very recent star formation activity in 9/11 quasars within our sample. We find an average star formation rate of 50 M yr −1 for the star-forming regions within our sample. The average dynamical mass of our quasar host galaxies of ∼ 10 11 M , indicates that the galaxies sit near the galaxy star formation rate -stellar-mass sequence at z∼ 2 (Tomczak et al. 2016). Using the average metallicity of 8.5 measured in dynamically quiescent regions and the average stellar mass of our sample indicates that our galaxies sit on the mass-metallicity relationship at z∼2 (Sanders et al. 2015).
Quasars at z∼ 2 are found to reside in galaxies with a broad range of star formation rates, spanning from quiescent to star-bursting galaxies. However, our sample preferentially contains quasar host galaxies in a starburst phase. High specific accretion rate AGN are more likely to be found in star-bursting galaxies with rates on or above the star formation rate -stellar-mass sequence in the distant Universe (Aird et al. 2019). We selected to observe compact steep spectrum radio-loud quasars, this class of objects tend to contain younger AGN. One of the mechanisms to trigger a luminous AGN is through a massive gas-rich galaxy merger (Treister et al. 2012). During the ongoing merger, the loss of angular momen-tum feeds gas into the centers of galaxies, providing fuel for both star formation and SMBH growth. Since we selected AGN that may have recently triggered, they are more likely to be in an ongoing merger, where star formation activity is enhanced. Indeed, about 7/11 of the quasar host galaxies in our sample are mergers. This can explain why our sample preferentially contains galaxies with active or recent star formation and rapid accretion onto the SMBH.
The measured star formation rates within our sample are significantly lower than those measured through dust emission in the far-infrared by the Herschel Space Observatory (Podigachoski et al. 2015;Barthel et al. 2017) for 4C04.81, 4C09.17, 3C318, and 3C298. The most likely explanation is that the quasar itself could partially heat the dust, Hα misses a significant fraction of the obscured star formation, or the dust traces a different star formation history. Interestingly for 3C298 and 3C318, where both high spatial resolution imaging of the dust and Hα emission is present, there is a significant misalignment between the maps. In places where we see evidence for recent star formation based on nebular emission-line ratios in 3C 298 and 3C 318, Barthel et al. (2018); Barthel & Versteeg (2019) does not see any dust emission. For the case of 3C 298 in the location where we see recent star formation traced by Hα , we also detect a molecular reservoir; however, no dust emission is present there. Furthermore, in the places where dust emission exists in the case of 3C 298, the molecular gas at that location is stable against gravitational collapse and has been on a time scale longer than the propagation of the observed outflow. For the case of 4C09.17 and 4C04.81, no high-resolution dust maps are available. The dust emission could originate at any location within the ∼ 17 Herschel SPIRE beam, which translates to a physical scale of about 150 kpc. Future high spatial resolution dust and molecular gas emission maps are necessary for proper comparison between the obscured and unobscured star formation traces and the molecular gas dynamics.

Offset from local scaling relations
The majority of our systems appear to be offset from both local scaling relationships between the mass of the SMBH and mass and the velocity dispersion of the bulge (see Figures 5, 6). To explain the large offset from the local M • − σ and M • − M * relationship, we could invoke a significant error in the estimated SMBH masses. The bolometric luminosities of some of our quasars are far greater than those used for reverberation mapping in the nearby Universe, which is used in calibrating the single epoch SMBH mass (Greene & Ho 2005). The SMBH masses would have to be off by 2-3 orders of magnitude to explain the observed offsets. By assuming that the SMBH grows primarily through gas accretion, we can use the Eddington luminosity formula to estimate the SMBH mass. Given that our quasars are most likely not all accreting at or close to the Eddington limit, this derived mass is effectively a lower limit.
For the derived bolometric luminosities in Table 2 we find a range of minimum SMBH of 10 7.5−9 M , consistent with what we measure from single epoch SMBH masses using the Hα emission line. Hence, ther is likely no significant error in our measured black hole masses.
In Figure 10, we plot the offset from the local scaling relation against the redshift of each object from our sample, the local galaxies sample with SMBH > 10 9 M and higher redshift quasars. Quasars with SMBH > 10 9 M appear to be offset from the local scaling relationship, which indicates that SMBH growth appears to outpace that of stars in these systems. The SMBHs may grow rapidly up to a mass of several times 10 9 M as early as 690 Myr after the Big Bang (Bañados et al. 2018), matching in mass to some of the most massive SMBH seen today. Some galaxies with lower luminosity AGN and lower mass SMBH also appear to be offset from the local scaling relation at z> 1 (Merloni et al. 2010;Bennert et al. 2011). Given the typically large uncertainty on the measured values and generally small sample sizes, it is difficult today to say whether a different population of AGN/galaxies are offset differently from the local scaling relationships at z> 1.
Under the assumption that SMBH primarily grows through Eddington-limited gas accretion, the growth is expected to be exponential. The e-folding or "Salpeter" time scale is about 50-300 Myr, depending on the SMBH spin. At the mean redshift of our sample (z=1.87), the SMBHs are expected to experience 30-200 e-folds in mass growth. However, for a duty cycle of around 10% (Wang et al. 2006) the expected number of e-folds drops down to about 3-20. Furthermore, the quasars in our sample are not accreting near the Eddington limit and can eventually switch from high to low accretion-rate mode, further decreasing the Eddington ratio. Hence, the SMBHs in our sample have nearly finished forming and will only further grow by a factor of 1.2-7 under the assumption of an Eddington ratio of 10%, and a duty cycle of 10%. If these galaxies are to assemble onto the local scaling relation and to evolve into the most massive early-type galaxies that we see today, then the rapid SMBH growths at early times in the Universe must be followed by significant stellar growth. On average, the galaxies within our sample need to grow the stellar mass within a radius of 7 kpc at a constant rate of 100 M yr −1 from z=2 to z=0.
In the host galaxy of 3C 298, there is currently insufficient molecular gas for the galaxy to grow in stellar mass to match the mass predicted by the local scaling relationship. Furthermore, the quasar 3C 298 does not appears to live in an over-dense environment based on the number count of galaxies seen with the Spitzer space telescope imaging data (Ghaffari et al. 2017). The open question is, how do these galaxies obtain the stellar mass necessary to grow into the massive galaxies we see today? Are minor mergers responsible for growing these galaxies? Alternatively, is the accretion of cool gas from the CGM responsible for providing the fuel necessary for future star formation? The results we find for the host galaxy of 3C 298 favor the scenario where cold accretion flows from the CGM will supply most of the fuel necessary for future star formation. Another scenario could be that the Spitzer observations are too shallow to see lower mass galaxies. If these systems are gas-rich, they can supply future fuel for star formation from merging the gas in their CGM and ISM with the quasar's host. Indeed in recent hydrodynamical simulation (Anglés-Alcázar et al. 2017) found that for dark matter halos with masses > 10 12.5 M majority of the mass build up happens from gas accreted from the CGM and transfer/exchange of gas from CGM and ISM of cannibalized low mass galaxies. These simulations also find that stellar build-up from dry mergers and just accretion of stars from merging galaxies is not significant to grow the stellar mass of galaxies in massive halos. If this is the case for the majority of our galaxies, it implies that they have enormous amounts of gas inside their CGM.
Our results can be in stark contrast to the predicted evolutionary paths of massive galaxies. In today's theoretical framework (Di Matteo et al. 2005;Hopkins et al. 2008;Zubovas & King 2012, feedback from the SMBH is predicted to happen once the galaxy reaches the local M • − σ relationship. However, our systems are showcasing outflows that are capable of causing major feedback when the mass of the galaxies is a fraction of their predicted final mass from the local scaling relations. Also, the gas-phase metallicities are far lower than those observed in nearby AGN. The kinetic luminosities for half of the outflows in our sample are far lower than the values predicted in simulations for the bolometric luminosities of our quasars (Vayner et al. 2021). Ionized outflows in other samples show similar results, where about half the objects lie below the predicted minimum energy-coupling between the quasar and the outflow of  (2013), log 10 (MBH /M ) = 8.32 + 5.64 log 10 (σ/200 km s −1 )). On the y-axis, we quantify the offset as the difference between the observed and predicted velocity dispersion from the local scaling relation based on the observed SMBH mass. We plot the observed offset from the local scaling relation against the redshift for individual targets. The labels are similar to 5. The shaded blue region represents the intrinsic scatter in the M• − σ relationship for black holes with a mass of 10 9.5 M . There is an overall offset for galaxies with massive SMBH at z>1 from the local M• − σ relationship. We find no statistically significant difference in the offset between any of the high redshift samples, while there is a statistically significant offset from the local BCG points (green).
0.1% at z∼ 2 (Vayner et al. 2021). If all these systems are offset from the local scaling relationship, it would be easier to launch the outflows because their masses are smaller compared to if they were on the local scaling relations. This could lead to lower energy coupling efficiency. On the other hand, we might be missing a significant fraction of the gas within the outflows because a large portion of the gas could be in either a molecular or neutral phase.
In the quasar host galaxy of 3C298, we find the majority of the gas in the outflow is in a molecular state, and once combined with the ionized kinetic luminosity we find values closer to those predicted in simulations. The kinetic luminosity in 3C298 is close to 1% of the quasar's bolometric luminosity. Regardless if we are accounting all the gas in the outflow, outflows capable of causing feedback are occurring before the galaxies are on the M • − σ relationship. We might need to reconsider our theoretical framework for massive galaxy formation, where the gas is not cleared from the galaxy in a single "burst" of feedback once the galaxies reach the M • − σ relationship. Instead, the SMBH grows first in massive dark matter haloes, followed by a delayed growth of the host galaxy with regulatory feedback from the SMBH and near-continuous accretion of gas from the CGM and nearby satellite galaxies. In such a scenario, the coupling efficiency might be lower per outflow event, compared to a single burst model where a single outflowevent clears all the gas. At later times, maintenance mode feedback from jets can heat the CGM, preventing gas from cooling and accreting onto the galaxy, keeping the galaxies "quenched".

CONCLUSIONS
We have conducted a near diffraction-limited survey of 11 quasar host galaxies to study the distribution, kinematics, and dynamics of the ionized ISM using the OSIRIS IFS at the W.M. Keck Observatory. This survey paper aimed to understand the source of gas ionization, the physical and chemical conditions of the ISM and to estimate the masses of radio-loud quasar host galaxies at z∼2. We detected extended emission in all objects on scales from 1-30 kpc and found that: • The AGN photoionizes the majority of the extended gas. A significant fraction of emissionline ratios are found to reside between the two sequences on the traditional BPT diagram. By applying evolutionary models of the mixing and star-forming sequence from z=0 to z=2.5, we find that the gas within our systems is denser and has lower metallicity compared to the gas photoionized in local AGN.
• In 9 objects, we find dynamically quiescent regions, with lower average log([OIII]/Hβ ) ratios. For systems where Hubble Space Telescope imaging is available, their morphologies are consistent with clumpy star-forming regions commonly observed in the distant Universe, indicating the presence of recent star formation. We find these systems to be forming stars at a maximum rate of 9-160 M yr −1 based on the Hα luminosity.
• For nine objects, we are able to measure the mass of the SMBH, the stellar velocity dispersion using the narrow component of Hα emission line as a proxy, and galaxy mass. We compare these nine objects to the local scaling relation between the mass of the SMBH and the mass or velocity dispersion of the galaxy. Our systems are both offset from the M • − σ and M • − M * relationship. Substantial growth is still necessary if these systems are to evolve into the present-day massive elliptical galaxies. Gas accretion from the CGM and gas-rich minor mergers are necessary to grow the stellar mass and increase the metallicity of the ISM. On average, the galaxies need to grow by at least an order of magnitude in stellar mass if they are to assemble onto the local scaling relations. A near-constant mass growth rate of ∼100 M yr −1 is necessary within a radius of 10 kpc from the quasar from z∼ 2 to 0.
• Combining the results of this paper with (Vayner et al. 2021) we find evidence for the onset of feedback before the galaxies are on the local M • − σ relationship. Luminous type-1 quasars are not the end phase of massive galaxy formation.
The authors wish to thanks Jim Lyke, Randy Campbell, and other SAs with their assistance at the telescope to acquire the Keck OSIRIS data sets. We want to thank the anonymous referee for their constructive comments that helped improve the manuscript. The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
Software: OSIRIS DRP (Larkin et al. 2013), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020 3C9 is a luminous quasar at z = 2.019922 with a prominent blue rest-frame UV continuum. For this source, we identify three distinct regions. "SE-SW component A" is a region with a ring-like morphology associated with the 3C9 quasar host galaxy. We measure a velocity dispersion from a Gaussian fit to the nebular emission lines of 407.6±12.9km s −1 and the kinematics resembling a rotating disk. "SE component A" is classified as an outflow region with a very high emission line FWHM of 1362.7±60.5 km s −1 and elevated log([OIII]/Hβ ) and log([NII]/Hα) ratios relative to the rest of the system. "N component B" is the merging galaxy in the 3C9 system showcasing a line FWHM of 472.15±11.8 km s −1 and a velocity offset of ∼200 km s −1 from the quasar. The projected spatial separation between the two apparent nuclei is 9 kpc. The quasar lies in the galaxy with a ring-like morphology showing the kinematic structure of a disk. Archival HST imaging of rest-frame UV continuum shows the ring morphology as well (see Figure Figure 11. Spectra of distinct regions along with fits to individual emission lines for the 3C 9 system.   Figure 14. Spectra of distinct regions along with fits to individual emission lines for the 3C 298 system. 4), indicating very recent star formation activity in the ring. The merging galaxy "N component B" appears to be a dispersion dominated system with active star formation and appears in rest-frame UV emission. The 3C9 system best resembles the local galaxy merger system Arp 148 (z=0.036), also known as Mayall's Object. The outflow in this system appears to be emanating from the ring of the galaxy with the quasar.  Figure 15. Spectra of distinct regions along with fits to individual emission lines for the 3C 446 system. 4C 09.17 is a luminous quasar at z=2.117 with a blue UV continuum. For this source, we identify four distinct regions. "SW component A" is a star-forming clump associated with the quasar host galaxy. The spectrum of this region shows a single narrow emission line with an FWHM of 312.0±7 km s −1 . "S/E component A" is an outflow region driven by the quasar, the nebular emission lines for this region have an FWHM of 887.2±22.4 km s −1 . A second narrow component is required for a good fit for each emission line in this region, with a line FWHM of 290.4±29.9 km s −1 . "W component B clumps" is a region part of the merging galaxy within the 4C09.17 system. The region consists of clumpy emission selected by isolating spaxels with an Hα line surface density > 6 × 10 −16 erg s −1 cm −2 arcsec −2 . "W component B diffuse" is emission associated with "diffuse" ionized emission in the merging galaxy selected by isolating spaxels with an Hα spatial line surface density < 6 × 10 −16 erg s −1 cm −2 arcsec −2 . The diffuse region shows higher log([OIII]/Hβ ) and log([NII]/Hα) line ratios associated with both AGN and star formation photoionization while the clumpy regions of the merging galaxy showcase lower ionization levels consistent with photoionization by star formation. This region is associated with bright UV emission in HST imaging of this object (Lehnert et al. 1999). "S/E component A outflow" shows high log([NII]/Hα) and log([OIII]/Hβ ) values relative to the rest of the system, indicating this region is predominantly photoionized by the quasar. The 4C09.17 system is a merger of two galaxies with velocity offsets of ∼1000 km s −1 and a projected separation of ∼ 4 kpc. HST imaging of rest-frame UV continuum (see Figure 4) shows evidence for a population of young hot stars indicating recent star formation activity. The majority of the star formation activity is confined to the merging galaxy.
C. 3C 268.4 3C 268.4 is a luminous quasar at z=1.39, with a slightly reddened UV continuum compared to the average type-1 quasar. For this target, we identified two distinct regions. "SW component A" is an outflow driven by the quasar. The FWHM of the emission lines is 2075±354 km s −1 as measured from the Gaussian fit to the [OIII] line. The spectrum extracted over this region also shows a narrow component with an FWHM of 603.7±54.9 km s −1 , most likely signaling emission from an extended narrow-line region close to the quasar. Because of issues with miss-assignment of flux in the OSIRIS pipeline (Lockhart et al. 2019), the rows below and above the centroid of the quasar do not have properly extracted spectra in the H band observations of this object. Hence we do not have a good spectrum of the extended emission in a 0.2-0.3 radius around the quasar in the H band, which covers the Hα and [NII] emission lines of the ionized outflow. "SW component B" is a region associated with the merging galaxy, showcasing clumpy morphology in ionized gas emission. The emission lines have an FWHM of 367.7 ± 3.9 km s −1 and an offset of −300 km s −1 relative to the redshift of the quasar. The log([OIII]/Hβ ) line ratios are lower for this region compared to the rest of the system, consistent with a mixture of AGN and star formation photoionization. This region is also associated with bright rest-frame UV continuum emission, seen with HST observations of this target Hilbert et al. (2016).

D. 7C 1354+2552
7C 1354+2552 is a luminous quasar at z=2.0064 with a blue UV continuum. For this target, we identify two distinct regions. "Component A" is the extended emission associated with the quasar host galaxy. The kinematics show a smooth velocity gradient, indicating the presence of a galactic disk. The size, morphology, and kinematics of the disk are similar to that of star-forming galaxies on the more massive end of the star formation main sequence at z ∼2 (Förster Schreiber et al. 2018). We measure an emission line FWHM of 357.2±2.0 km s −1 on the redshifted side of the disk and 497.7±6.5 km s −1 on the blue-shifted side of the disk. Although this region only has a single label ("component A"), in Figure 13 rows one and two show the fits to the red and blue-shifted sides of the disk that are part of this region. This region is selected based on the location where Hα emission is detected. This is done to boost the SNR in the Hα line as it appears to be clumpier, more compact, and less extended than [OIII] . In Table 3 we provide values integrated over the entire galactic disk. "E component B " is a region associated with the merging galaxy at a projected separation of 6-7 kpc. The kinematics are consistent with a dispersion dominated galaxy. The entire "component A" is consistent with quasar photoionization. The gas in "E component B" is photoionized by star formation.
E. 3C 298 3C298 is a luminous quasar at z=1.439 with a slightly reddened UV continuum. For this target, we identify five distinct regions. We present a detailed discussion of each region in Vayner et al. (2017). "W/E component A" are outflow regions with a bi-conical morphology, where the western (W) is the redshifted receding cone, and the eastern (E) is the approaching cone. In Vayner et al. (2017), they are referred to as the red(blue) shifted outflow region. The emission lines over the outflows are very broad, with FWHM up to ∼1500 km s −1 . A combination of shocks and quasar activity is likely responsible for photoionizing the gas. "SE component B outflow" is an outflow region belonging to a merging galaxy. "SE component B ENLR" is an extended narrow-line region belonging to the disk of the merging galaxy, with gas photoionized by the quasar or secondary AGN. "SE component B Tidal feature" is a region of the merging galaxy with active/recent star formation as evident by lower log([NII]/Hα) and log([OIII]/Hβ ) values compared to the rest of the regions.

F. 3C318
3C318 is a luminous quasar at z=1.5723 with a reddened UV continuum. There is evidence for a spatially unresolved nuclear star-burst with an upper limit on the star formation rate of 88±9M yr −1 . This star formation rate is far lower than the far infrared derived rate of 580M yr −1 . The extinction towards the nuclear region measured from Willott et al. (2000) alone cannot explain the mismatch between the Hα and far-infrared derived SFR. Either a larger fraction of the far-infrared emission is from dust that is being heated by the AGN, or the far-infrared emission traces a different star formation history than Hα (Calzetti 2013). No narrow extended emission is detected in this object.
The merger status of this object is unclear. Two nearby galaxies to the north and west of the quasar are visible in archival HST imaging (Willott et al. 2000). We do not detect the western object that is 2 away from the quasar in our OSIRIS observations in any emission line. Willott et al. (2007) studied this object with PdBI through CO 2-1 emission at a fairly coarse (∼8 arcseconds) resolution. There appears to be CO emission that could plausibly be associated with the western object. We have recently obtained a much higher angular resolution CO 3-2 spectroscopy of this target that will be discussed in detail in a forthcoming paper. We confirm the existence of CO 3-2 emission associated with the CO 2-1 emission. We resolve the molecular emission into multiple components. However, the CO 3-2 emission is not associated with either one of the galaxies seen in the HST data. We obtained a wide field of view IFS observations of this target with KCWI aimed at attempting to measure the redshifts of the nearby galaxies and to confirm the merger scenario of this object. We detect both the northern and western objects in the continuum. We confirm that the northern target is at a different redshift than the quasar from the detection of [OII] emission, while for the western object, a reliable redshift is challenging to determine with the current data set. Hence no clear evidence of a companion galaxy that is part of a merger is detected for this quasar associated with the brightest galaxies seen in optical imaging within a few arcseconds from the quasar.
G. 4C 57.29 4C 57.29 is a luminous quasar at z=2.1759 with a blue UV continuum. For this target, we identify two regions. Region "NE component A" belongs to the host galaxy of the quasar. The relatively high log([OIII]/Hβ ) value indicates that this region is consistent with being photoionized by the quasar. The 500.7 nm [OIII] is the only emission line detected for this region. The region is marginally resolved, making it hard to measure the kinematic structure. We require a double Gaussian fit to the [OIII] emission in this region to obtain a good fit, and we measure an FWHM of 474.3 and 502.5 km s −1 with offsets of 35.0 km s −1 and -1050.1 km s −1 relative to the redshift of the quasar. We identify a second region north of the quasar. It is unclear if it belongs to a merging galaxy or the quasar host galaxy. There is a ∼ 100 km s −1 offset from the quasar, and the line has an FWHM of 550.13∼19 km s −1 . This region is also only detected in [OIII] . The SNR is too low to measure any kinematics structure.
H. 4C 22.44 4C22.44 is a luminous quasar at z=1.5492 with a reddened UV continuum. Similar to 3C318, we do not detect any evidence for a merging galaxy for this system. For this target, we identify a single region, "N, S component A". The kinematics of this region may be consistent with a galactic disk belonging to the quasar host galaxy. We see evidence for a smooth gradient in the radial velocity map, however, the region is marginally resolved. We measure an emission line FWHM of 434.8 km s −1 . The region is consistent with being ionized by star formation with some contribution from quasar photoionization.
I. 4C 05.84 4C05.84 is a luminous quasar at z=2.323 with a slightly reddened UV continuum. For this target, we identify three distinct regions. Regions "S component A" and "NE component A" are the blue(red) shifted outflow regions resembling a bi-conical outflow. They showcase broad extended emission with a line FWHM of ∼800 km s −1 . The quasar photoionizes these regions. Region "SW component A clump", shows a line FWHM of 467.9±3.0 km s −1 and is photoionized by a combination of star formation and the quasar. This clump is also detected in NIRC2 imaging of this object studied by Krogager et al. (2016), where they consider this clump to be associated with a damped Lyα system. However, here we confirm that this objected is part of the quasar host galaxy. We find no evidence for a merging galaxy within our OSIRIS observations. J. 3C 446 3C446 is a quasar at z = 1.404. For this target, we identify two regions, "N component A tidal feature" is a region belonging to the quasar host galaxy, resembling a tidal feature that is most likely induced by the merger. We measure an FWHM of 395.14±2.0 km s −1 for this region. "E-W component B" belongs to the merging galaxy, a portion of it resembles a tidal feature, counter to the tidal arm of "N component A tidal feature." For this region, we measure a line FWHM of 558.5±63 km s −1 however, it appears to be a blend of two velocity components. It is unclear where the nucleus of the merging galaxies resides. It could be that it has already merged with that of the quasar host galaxy. The two galaxies appear to be offset by at least 500 km s −1 in velocity, and there is a possibility that a portion of the merging galaxy lies on top of the quasar host galaxy.
K. 4C 04.81 4C04.81 is a luminous quasar at z=2.5883 with a reddened UV continuum. For this target, we identify a single region, "E component A outflow". The kinematics show blue and red-shifted broad (FWHM∼800 km s −1 ) emission. The quasar mainly photoionizes the gas. We do not identify any narrow extended emission in this object.