Powerful Radio-loud Quasars Are Triggered by Galaxy Mergers in the Cosmic Bright Ages

While supermassive black holes are ubiquitous features of galactic nuclei, only a small minority are observed during episodes of luminous accretion. The physical mechanism(s) driving the onset of fueling and ignition in these active galactic nuclei (AGN) are still largely unknown for many galaxies and AGN-selection criteria. Attention has focused on AGN triggering by means of major galaxy mergers gravitationally funneling gas toward the galactic center, with evidence both for and against this scenario. However, several recent studies have found that radio-loud AGN overwhelmingly reside in ongoing or recent major galaxy mergers. In this study, we test the hypothesis that major galaxy mergers are important triggers for radio-loud AGN activity in powerful quasars during cosmic noon (1 ≲ z ≲ 2). To this end, we compare Hubble Space Telescope WFC3/IR observations of the z > 1 3CR radio-loud broad-lined quasars to three matched radio-quiet quasar control samples. We find strong evidence for major-merger activity in nearly all radio-loud AGN, in contrast to the much lower merger fraction in the radio-quiet AGN. These results suggest major galaxy mergers are key ingredients in launching powerful radio jets. Given many of our radio-loud quasars are blue, our results present a possible challenge to the “blowout” paradigm of galaxy evolution models in which blue quasars are the quiescent end result following a period of red quasar feedback initiated by a galaxy merger. Finally, we find a tight correlation between black hole mass and host galaxy luminosity for these different high-redshift AGN samples that is inconsistent with those observed for local elliptical galaxies.


INTRODUCTION
Essentially all massive galaxies seem to harbor a supermassive black hole (SMBH, taken to be black holes with masses M BH ≳ 10 5 M ⊙ ) at their dynamical centers (Kormendy & Ho 2013).Yet only a small fraction of these SMBHs are typically observed during a phase of luminous accretion characteristic of active galactic nuclei (AGN).The unified model for AGN stipulates that accretion of matter onto the SMBH is the primary energy source for the observed radiative output from these systems (see e.g.reviews by Antonucci 1993;Netzer 2015).The accretion of gas and dust (either directly, or indirectly through the tidal stripping of stars), and subsequent viscous dissipation within the developing accretion disk, allows for the conversion of gravitational potential energy into thermal optical-UV disk radiation.The accretion disk can outshine the host galaxy at visible wavelengths and its ionizing continuum can excite atomic tranisitions in gas clouds both near and far from the SMBH, giving rise to the so-called broad and narrow emission lines seen in type I and II AGN, respectively1 .
However, uncertainty persists on how AGN are produced in the context of galaxy evolution and dynamics, where there are various physical processes which might contribute to supplying SMBHs with sufficient gas densities for the initiation of AGN activity.A large focus in the literature is on whether major galaxy mergers might be the dominant AGN triggering mechanism (especially at high AGN luminosities) through galaxy-wide, gravitationally-induced torques which drive gas towards the galactic center (as bolstered by various semi-analytic and numerical works, e.g., Hopkins et al. 2006).Curiously, there have been a substantial number of studies which both seem to support (e.g., Koss et al. 2010;Ellison et al. 2011;Hong et al. 2015;Weston et al. 2017;Goulding et al. 2018;Gao et al. 2020) and oppose (e.g., Cisternas et al. 2011;Kocevski et al. 2012;Karouzos et al. 2014;Villforth et al. 2019;Lambrides et al. 2021) the relevance of major galaxy mergers in triggering AGN.Collectively, these conflicting results suggest that the physical properties of both AGNs and their host galaxies, and extragalactic environments, may be pertinent when considering the role of galaxy mergers towards their initial triggering (or even sustained fueling).Additionally, they highlight the need to look at alternatives to major galaxy mergers for transporting gas from kpc to sub-pc scales in galactic nuclei and ulti-mately triggering AGN, namely: internal "secular" processes, minor/satelite-galaxy accumulation (e.g., Hernquist & Mihos 1995), or gas inflows induced from other interactions within a larger-scale galaxy cluster environment besides galaxy mergers.
For AGN residing in clusters, ram-pressure induced triggering (Marshall et al. 2018), "galaxy harrassment" from high-speed galaxy encounters (Moore et al. 1996), and tidal interactions between the host galaxy and gravitational potential of the cluster (Byrd & Valtonen 1990) are all potential AGN-triggering mechanisms apart from major galaxy mergers.Among secular processes, common channels for SMBH growth include stellar winds (Davies et al. 2007;Vollmer et al. 2008, although see Dubois et al. 2015 who showed supernovae can inhibit black hole growth in low-mass systems), chaotic and cold accretion streams driven by turbulence in the host galaxy intersellar medium (ISM) (Hobbs et al. 2011;Gaspari et al. 2013), and gravitational disk instabilities mediated by stellar bars (and spiral density waves) (e.g., Shlosman et al. 1989), or dense, inhomegeneous clumps formed through cold accretion streams from the larger-scale halo environment (the so-called "violent disk instabilities" which seem to be more efficient at higher redshifts, see e.g., Dekel et al. 2009;Bournaud et al. 2011;Gabor & Bournaud 2013).
Because the triggering of AGN is intimately linked to the growth of SMBHs, studies of AGN triggering will help us better understand the tight scaling relations observed between SMBH mass and various host galaxy properties, typically of the bulge or "spheroid" component (e.g., Magorrian et al. 1998;Gebhardt et al. 2000;Ferrarese & Merritt 2000;Marconi & Hunt 2003;Häring & Rix 2004;Gültekin et al. 2009;Graham et al. 2011;McConnell & Ma 2013).The origin of these correlations is still debated, but leading non-mutually exclusive hypotheses include the hierarchical buildup of both black holes and bulges through galaxy mergers (Peng 2007;Jahnke & Macciò 2011), stellar feedback through starbursts (Norman & Scoville 1988;Davies et al. 2007;Wild et al. 2010), and AGN feedback ( e.g., Fabian 2012;Heckman & Best 2023).AGN feedback can take the role of self-limiting SMBH growth through the ejection or heating of gas by radiatively-coupled outflows or jets (e.g., Nesvadba et al. 2010;Couto & Storchi-Bergmann 2023), or positive feedback where SMBH growth is coupled to host star formation by compressing the ISM (e.g., Silk 2013).Perhaps the most striking example of AGN feedback is the heating of the intracluster medium (ICM) by large-scale radio jets hosted by massive cD galaxies (see e.g. the famous example of the Perseus cluster, Fabian et al. 2003Fabian et al. , 2006) ) Several previous studies have found strong evidence for enhanced incidence of major galaxy mergers among samples of radio-loud AGN (Heckman et al. 1986;Colina & de Juan 1995;Ramos Almeida et al. 2012;Ivison et al. 2012;Kaviraj et al. 2015;Chiaberge et al. 2015;Noirot et al. 2018), and in constrast to the less conclusive results for AGN samples selected by other criteria (e.g., Grogin et al. 2005;Gabor et al. 2009;Georgakakis et al. 2009;Rosario et al. 2015;Lambrides et al. 2021;Sharma et al. 2021).Understanding why galaxy mergers seem to be ubiquitous among radio-loud AGN, and not among other AGN samples, may help give us better insight into the dyanamical (co)evolution between SMBHs, AGN, and their host galaxies.Furthermore, studying the link between galaxy mergers and radio-loud AGN may help give clues to the jet formation physics operating in these systems.
The main goal of this study is to test whether major galaxy mergers are an important process for the triggering of radio-loud AGN.To this end, we analayzed the host galaxy stellar morphologies of the z > 1 (all but two 1 < z < 2) 3CR sample of broad-lined radio-loud quasars with Hubble Space Telescope (HST) WFC3/IR images.This work is a follow-up study to that of Chiaberge et al. (2015), who analyzed the host galaxies of twelve type II radio galaxies drawn from the same z > 1 3CR parent sample, observed with HST WFC3/IR SNAP observations (program SNAP13023).Chiaberge et al. (2015) found much higher merger fractions among radio-loud AGN in comparison to matched samples of inactive galaxies and radio-quiet AGN − essentially all z > 1 3CR radio galaxies were mergers in comparison to a merger fraction of ∼ 40% for the control samples.In this work, we focus our analysis on the 292 type I (broad-lined) quasars from the z > 1 3CR sample, in order to complement the type II radio galaxies studied in Chiaberge et al. (2015).The observations require HST point spread function (PSF) subtractions of the unresolved quasar light in order to examine the host galaxy stellar morphologies, as we discuss further in section 3.2.We then used human experts to blindly classify HST WFC3/IR images of our 3CR sample in addition to three matched radio-quiet quasar3 control samples in order to investigate the role of galaxy mergers in triggering the formation of powerful radio jets.The shaded black histogram represents the 3CR sample, orange represents the M19, purple the V17, and cyan the M16 sample.See section 2 for a description of each of these control samples.

SAMPLE PROPERTIES
The z > 1 3CR sample we consider comprise 64 objects observed in the northern hemisphere in Cambridge, United Kingdom (Spinrad et al. 1985).The 3CR sample is a low-frequency, flux-limited sample (F 178 MHz > 9 Jy) of extragalactic objects (|b| > 10 • and declinations above 10 • ), unbiased in jet orientation because the radio emission at such low frequencies is dominated by the isotropic lobe contribution.These 3CR galaxies host relatively young radio jets (Murgia et al. 1999;Dallacasa et al. 2021;O'Dea & Saikia 2021), with edge-brightened FR type II morphologies (Fanaroff & Riley 1974), and include some of the most powerful radio-loud AGN known.Furthermore, the z > 1 3CR quasars reside in the so-called "Bright Ages" (1 ≲ z ≲ 2, an epoch also referred to as "cosmic noon"), thought to correspond to the cosmic peak of AGN activity, star formation, and galactic bulge assembly (e.g., Alberts et al. 2016;Hopkins et al. 2006), making them an indispensible sample for examining the relationship between SMBHs and their host galaxies.
Our radio-quiet quasar control samples were constructed from three separate archival AGN samples of ∼ 20 sources observed with similar depth and filter HST WFC3/IR observations, as shown in Table 1.All control quasar samples have published papers assessing the merger fraction of the quasar samples relative to matched samples of inactive galaxies (Mechtley et al. 2016;Villforth et al. 2017;Marian et al. 2019).None of these studies found evidence for an enhanced merger fraction among the quasars, thus casting doubt on the scenario of merger-driven triggering for those AGN.The selection criteria of our three control samples are as follows: • The first control quasar sample is the high-redshift sample of Mechtley et al. (2016) (1.9 ≤ z ≤ 2.1), referred to as "M16" for the rest of the paper.This sample was selected from the Sloan Digital Sky Survey (SDSS) DR5 quasar catalog (Schneider et al. 2007).This sample of 20 broad-line quasars was selected to have high black-hole masses, spanning the range 9.3 ≤ log (M BH ) ≤ 9.7 (as determined from virial estimates based upon the Mg II broad line from Shen et al. 2011), and a uniform color-selection algorithm as described by Richards et al. (2002).Bolometric luminosities from this sample are from Shen et al. (2011).
• The second sample is that of Villforth et al. (2017), which we will refer to as "V17".This sample is comprised of 20 X-ray selected AGN (ROSAT Xray detections) crossmatched with SDSS DR5 optical quasar broad-line detections (Anderson et al. 2007).The redshift range spans 0.5 ≤ z ≤ 0.7, and the quasars were also selected to have bolometric luminosities in excess of 10 45 erg s −1 .Black hole masses and bolometric luminosities for this sample are from Shen et al. (2011).
• Our final control sample is composed of the 1.8 ≤ z ≤ 2.2 broad-line quasars from Marian et al. (2019), which we refer to as "M19".These were sub-selected from the SDSS DR7 quasar catalog (Schneider et al. 2010;Shen et al. 2011).These 21 AGN were selected to have 8.5 ≤ log (M BH ) ≤ 8.7 (as determined by the Mg II broad-line virial estimates from Shen et al. 2011), and more impor-tantly Eddington ratios > 0.7.Here the Eddington ratio, λ, is defined as λ ≡ L bol /L Edd , where L bol is the bolometric luminosity, and L Edd is the Eddington luminosity4 describing the limting luminosity needed to halt spherical accretion.Taken together (and assuming some fixed accretion efficiency), these assumptions imply high mass-accretion rates and thus arguably the most promising selection criteria for constraining the AGN triggering mechanism(s) in these systems.This sample has the same uniform color-selection criteria as M16.
Figure 1 shows the redshift distributions for our four samples.The 3CR radio-loud AGN fall between the V17 low-redshift (z ∼ 0.7) sample and the two high-redshift (z ∼ 2) samples (M16 and M19).The control samples have negligible redshift ranges and essentially sample distinct cosmic epochs close to, but generallly on either side, of the distributed 3CR cosmic ages.The V17 low-z sample has a mean cosmic age of ∼ 7.3 Gyr, while the high-z M16 and M19 samples have a mean cosmic age of ∼ 3.3 Gyr.The mean (and median) cosmic age of our 3CR radio-loud AGN sample is ∼ 4.6 Gyr, corresponding to z = 1.4.While the redshifts of our samples do not exactly match, the control groups effectively sample towards the beginning and right after the cosmic evolution of our 3CR sources.This gives us some leverage towards discriminating enhanced incidence of galaxy mergers among radio-loud sources while also considering potential redshift evolution in our interpretations.
In order to verify that the control samples are all radio-quiet, we cross-matched each against the FIRST, NVSS, and VLASS Very Large Array (VLA) all-sky surveys.These surveys have the required depth and sky coverage to assess the radio-loudness of our quasars, where we parametrize radio-loudness by the rest-frame ratio of 5 GHz radio to B-band optical flux density after K-corrections5 and use a cutoff of ten to delineate radioloud from radio-quiet objects (Kellermann et al. 1994).One of the M19 sources, one of the M16 sources, and five of the V17 sources were radio-loud.We removed these radio-loud sources from our control samples.
Figure 2 compares the black hole masses, bolometric luminosities, quasar colors, and Eddington ratios between the 3CR radio-loud AGN and the control radioquiet samples.For the 3CR, we compiled all black hole massess from different literature sources as given in Table 2. 3CR bolometric luminosities were determined from quasar PSF magnitudes obtained during our Galfit decompositions described further in section 3.2 and given in Table 2. Bolometric correction factors were obtained from Azadi et al. (2022) for radio-loud quasars at 1 < z < 2, after redshifting the F140W filter to the appropriate rest-frame values for each source.Black hole massess and bolometric luminosities are taken from the reference papers for the control samples; values are given in Table 3. Quasar U -V colors were determined in the rest frame, using nominal effective wavelengths of 245 6 nm and 580 nm for U and V band, respectively (U -V colors are given in Tables 2 and 3).In order to obtain reliable rest-frame color estimates, we used the appopriate observer-frame photometric band which most closesly redshifts to the rest frame U and V bands.To this end, we utilize SDSS DR14 ugriz photometry (Pâris et al. 2018), near-IR JHK photometry from the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007), and the HST observations described further in section 3. Once the appropriate photometric band is chosen, we K-corrected to the nominal rest-frame wavelengths using the near-IR-optical quasar template spectrum from Glikman et al. (2006).For the SDSS and HST data, our quasar colors were constructed using PSF magntiudes in order to avoid any host galaxy contamination (these are not available for UKIDSS data products).
In general, the distributions shown in Figure 2 cover similar regions of the overall parameter space.We also performed two-tailed Kolmogorov-Smirnov (KS) tests to assess whether the control samples and 3CR sample were consistent with being drawn from the same parent population (for all parameters of interest).The only sample comparisons in which we were able to reject the null hypothesis that the data came from the same population distributions at the 99% or greater significance level are the M19 vs 3CR black hole masses, M19 vs 3CR Eddington ratios, M16/V17 vs 3CR bolometric luminosities, and V17 vs 3CR quasar colors.The fact that the KS tests for the M19 vs 3CR black hole masses and Eddington ratios imply different parent distributions can be understood in terms of the M19 selection critera creating artifically narrow distributions.However, these M19 distributions are still fairly representative of the corresponding 3CR distributions.Similarly, while the M16 vs 3CR black hole masses and V17 vs 3CR quasar colors have different distributions, the range of values in the control samples is still fairly representative of the 3CR quasars.However, the V17 bolometric luminosities do seem to generally be lower than the 3CRs by ∼ an order of magnitude, even though there is substantial overlap.Along with the lower redshifts of the V17 sample, this is another factor to consider when comparing the V17 against 3CR merger fractions.The V17 sample was selected to have high bolometric luminosities (even if they do not approach the levels of the other samples), and were shown in V17 to be in excess of 10 45 erg s −1 where major mergers are expected to play a dominant role in AGN triggering based on theoretical analytical and simulaiton-based arguments.

HST OBSERVATIONS & DATA ANALYSIS
All of our quasar samples were observed with HST WFC3 near-IR observations.Table 1 gives the project code, exposure time, HST WFC3/IR filter, pivot wavelength, and surface brightness sensitivity limits for each sample.The sensitivity limits were estimated with the WFC3/IR online Exposure Time Calculator (ETC).For the ETC exposure calculations we used a 2 × 2 pixel extraction area, elliptical galaxy spectrum, median redshift for the sample, and neglected Milky Way extinction.The HST observations of the 3CR sample also included WFC3/UVIS F606W (pivot wavelength of 588.92 nm) observations in the same orbit as the IR channel exposures.We used UVIS aperture photometry when appropriate for determining quasar colors.Generally, there is fairly low host contamination in the UVIS channel because it probes rest-frame near-UV for the 3CR sample, thus allowing for robust quasar magnitude measurements from aperture photometry without PSFfitting (except for a few cases of star-forming regions, which we masked when measuring fluxes).
All samples were observed with ≥ 3 point dithers so as to remove hot pixels/cosmic rays and allow for improved resolution in the final sub-pixel imaging.All observations were non-destructively read out in a log-linear STEP sequence, intended to allow for high-dynamic range imaging.The exposure times for all samples allow for similar surface brightness detection thresholds, although the control samples are slightly deeper than the 3CR.We are testing the hypothesis that the 3CR radioloud AGN have a higher incidence of galaxy mergers, and greater sensitivity to low surface-brightness merger indicators in the control groups only serves to make our results more robust in the event of a higher 3CR merger fraction.
The 3CR sample was observed with the F140W wide filter (JH gap), with pivot wavelngth of 1392.3 nm.At the median redshift of the 3CR sample (z = 1.4), this corresponds to a rest-frame wavelength of ∼ 580 nm or V-band.The control-group samples were all observed with the wide H-band F160W filter, corresponding to a pivot wavelength of 1536.9 nm.For the V17 sample, this corresponds to a rest-frame wavelength of ∼ 904 nm or Z-band in the near-IR.For the high-redshift M16/M19 samples, this correpsonds to a rest-frame wavelength of ∼ 512 nm or V-band.The V17 observations were chosen to sample ∼ 1 µm rest-frame wavelengths in order to serve as a good proxy for stellar mass.The host galaxy morphologies should not appear substantially different between this wavelength and visible wavelengths, nor should the hosts be substantially contaminated by emission lines at these wavelengths, even if they are especially star-forming (Martins et al. 2013).As shown by Hilbert et al. (2016), the host galaxies of our 3CR sample should also not be substantially contaminated by emission lines, and the same is true for the z ∼ 2 M16/M19 samples which overlap in redshift with the 3CR.

Data Reduction
All flt image files were downloaded from the Mikulski Archive for Space Telescopes (MAST) after having first been processed through the calwf3 data reduction pipeline (Dressel 2022).We combined the dithered flt image files together with Astrodrizzle (Gonzaga et al. 2012), which also subtracts the sky background and makes corrections for cosmic ray artifacts and geometric distortion.Astrodrizzle relies on the "drizzle" algorithm, or variable-pixel linear reconstruction, in which pixels in dithered CCD images are mapped onto a sub-sampled output image after taking into account rotational and translational shifts between dithered input images and camera distortion (Fruchter & Hook 2002).We used a final_pixfrac parameter of 0.8 (describing the drop size of input pixels) and a final_scale of 0.06 ′′ for the final pixel scale of the output image.After experimentation, we found this combination of final_pixfrac and final_scale allowed for the best compromise between final image resolution, correlated noise properties, and sensitivity to low-surfacebrightness features.This 0.06 ′′ per pixel final plate scale of our output images oversamples the point spread function (PSF) by a factor of two for the F140W and F160W filters of HST's WFC3 camera.

Galfit Decompositions
After drizzling we performed 2D modeling of the brightness distribution within our images using the Galfit (Peng et al. 2010) software package.For simplicity and consistency, we used single Sérsic profiles to model the host galaxy light distrbution (Sérsic 1963;Sersic 1968), which have the following radial dependence: Here, I is the intensity as a function of radius (R), R e is the effective radius containing half the flux, I e is the intensity at R e , n is the Sérsic index, and κ is a parameter depending only on n.Thus, n solely determines how centrally concentrated the light profile is.Profiles with n = 1 correspond to exponential functions common to modeling galaxy disks (Freeman 1970), and n = 4 describes the de Vaucouleurs profile used to model classical bulges (de Vaucouleurs 1948).
While the Sérsic profile may not perfectly describe the light distributions of our galaxies, to the first order it is efficient in modeling the underlying brightness distribution of the non-disturbed host galaxy components without the worry of over-fitting.Since Galfit optimizes model parameters simultaneously when minimizing χ 2 through the non-linear "Levenberg-Marquardt" algorithm, this simple host light model helps keep us from oversubtracting the quasar PSF contribution described further below.
In addition to Sérsic profiles for the host galaxies, we modeled the sky background, and AGN point sources using synthetic PSFs constructed from the Tiny Tim software package (Krist et al. 2011).Our PSF model was constructed using the near-IR-optical quasar template spectrum from Glikman et al. (2006), redshifted appropriately for each source.We built PSF models separately for each dithered flt frame, using the WFC3 focus model from the following URL maintained by STScI: focustool.stsci.edu/cgi-bin/control.py.
Each PSF was constructed at the chip pixel location of peak flux measured in that flt frame, using a 7.5 ′′ diameter.After constructing each PSF model for the different dithered frames, we combined the PSFs by inserting them into blank flt frames and drizzling using the same parameters as our science images, except with a 2× finer plate scale of 0.03 ′′ per pixel.This yielded PSF models sampled to a ∼ 4× finer pixel scale than the PSF FWHM.We evaluated the quality of our PSFsubtraction technique for one of our quasars with an unresolved host galaxy in Appendix section B. We show all of the HST WFC3/IR images of our 3CR quasars before and after PSF-subtraction in Appendix section C. Similar use of drizzled Tiny Tim PSFs to subtract quasar light for host galaxy morphological analysis in other studies showed good results (Villforth et al. 2017;Zakamska et al. 2019).
We attempt fitting each source with both a Sérsic and PSF model component, allowing all relevant parameters to be free unless constraining some to physically motivated values helped improve the fit as determined by χ 2 ν (or allowed the fit to converge to physically meaningful values at all).The model parameters for the Sérsic profiles are the index n, effective radius R e , the axis ratio (semi-minor/semi-major axis), position angle, sky position, and integrated magnitude.For the PSF model, only the sky position and magnitude were free parameters, and the sky background was fit as a constant.We also used a hand-constructed bad-pixel mask for nearby sources and irregular host galaxy features so as to most accurately model the non-disturbed host and quasar components.Appendix section D compares this method of masking all pixels containing significant deviations from the quasar and its undisturbed quasar host galaxy component with an alternative method where we include a greater number of model components for very nearby or blended companion galaxies.The difference in magnitude between the two methods is less than the uncertainties found by Simmons & Urry (2008) for quasar host galaxies in simulations of HST ACS data for the GOODS survey (when comparing quasars with similar fluxes and flux ratios between the quasar and its host galaxy).For some cases of very bright quasars, we also masked the central few brightest pixels in order to account for non-linear pixel response where the PSF model would not match the saturated pixels.Tables 2 and 3 report the best-fit results of our Galfit decompositions.
In cases where the host galaxy appears to be marginally resolved or unresolved, we used a fixed Sérsic index of four and placed the Sérsic component at the same location as the PSF after optimizing a PSF-only fit.If the Sérsic effective radius was less than 3 pixels, equivalent to 0.18 ′′ or ∼50% beyond the PSF FWHM, then we considered the host unresolved and used only a PSF-fit in our Galfit galaxy model.These cases correspond to null values for m host in Tables 2 and 3.

CLASSIFICATION METHODOLOGY
For consistent comparison of our results to those of Chiaberge et al. (2015), who analyzed the host galaxies of the z > 1 type II 3CR radio galaxies, we used six human experts for our morphological merger classifications corresponding to six of the authors on this paper.This has the benefit of taking a more holistic view of galaxy mergers in comparison to quantitive, non-parametric measures of disturbance such as the Gini coefficient, second-order moment of the brightest 20% of the light (M 20 ), or asymmetry measure (Schade et al. 1995;Lotz et al. 2004).These quantitative estimates of disturbance can be effective at measuring post-galaxy-coalescence signatures (i.e., multiple nuclei, tidal features, and other large-scale asymmetries, see e.g., Lotz et al. 2011), but would miss many ongoing/incipient galaxy mergers or those with faint signatures that the human eye is more sensitive to.All measures are likely to be biased by the PSF-subtraction uncertainties, which are inherent to our classification images, as pointed out by previous WFC3/IR quasar merger studies (i.e., Glikman et al. 2015;Zakamska et al. 2019).
Each human classifier examined the entire image set before final source classifications were decided upon to ensure internal self-consistency.Each expert classified each quasar image presented in the form of Galfit panel decompositions.Figures 3 and 4 show examples.
These panels were constructed based upon our Galfit 2D modeling described in section 3.2, and show the original image data, our best-fit Galfit models, the PSFsubtracted images, and the residuals.As described in section 3.2, for unresolved objects for which we could not obtain a physically meaningful Sérsic fit, we only used a PSF model, and the residuals are the same as the PSFsubtracted image.Otherwise, the residuals are the data minus the best-fit PSF and Sérsic models.Including the residuals in the classifications helps to highlight the faint merger signatures (i.e., global asymmetries and tidal disturbances indicative of a recent or ongoing major galaxy merger), because the host galaxy components consistent with symmetric Sérsic profiles are removed.The images were stretched to highlight faint morphological features necessary for careful inspection and classification.Each panel was assigned a random number, and the experts had no identifying knowledge for each panel, so these are completely blind classifications.We classified images according to the following criteria: • A: Ongoing galaxy merger − these show two (or more) massive galaxies with clear signs of interaction such as tidal bridges, streams, shells, or tails.Cases with unresolved host galaxies that fit this category are still classified as ongoing mergers.
• B: Close companion(s)/incipient galaxy merger − these show a massive galaxy companion within the 25 kpc radius surounding the quasar but otherwise do not show any signs of interaction.Major/minor incipient galaxy mergers are further distinguished after classification by a flux ratio of 1:4 (i.e., minor incipient mergers are those where the less massive galaxy is ≥ 4× less massive).Minor incipient mergers were reclassified as non-mergers.
• C: Multiple nuclei − these show a common envelope around two or more nuclei, either stellar or unresolved quasar-like nucei.
• D: Post-merger Signatures − these do not have galaxy companions within 25 kpc but have host galaxies which show light profile asymmetries, tidal features, or other disturbances indicative of a post-coalescence major galaxy merger.
• E: Non-merger − these appear as smooth and symmetric host galaxies without any massive companions within 25 kpc.
• F: Host galaxy unresolved/non-detected or PSF-dominated residuals − for these, the clas-  sifier does not identify enough host galaxy component to classify it as any of the above letters.It may be some combination of very bright PSF, general PSF mismatch, compact host galaxy, undermassive host galaxy, high-redshift cosmological surface brightness dimming, or other factors that lead to this classification.
Figures 3 and 4 show examples firmly classified into each of the above categories.Each classifier's votes were consolidated into the categories of "merger" or "not merger", where categories of A−D counted towards "merger" and E/F were counted towards "not merger".This choice to reclassify unresolved sources as nonmergers implicitly assumes we would have seen evidence of a merger if it were there, regardless of incomplete host galaxy information.We discuss any potential bias resulting from this choice in Appendix section E. In the case of sources classified as F, or "unresolved", with detected companions within 25 kpc meeting the 1:4 flux ratio cutoff we use to discrimnate massive companions indicative of an incipient major galaxy merger, we reclassified them as "B" after voting and these would be considered mergers (see Appendix section F where we explain how this is done for cases without host galaxy fluxes determined from our Galfit decompositions).Finally, we further consolidated the votes into a consensus classification based upon the majority vote, where ties were broken in favor of "merger".

Merger Incidence
Table 4 gives the results of our consensus classifications after votes were consolidated into the categories "merger" and "not merger".In order to robustly quantify the 3CR merger fraction and compare to those of the control samples, we used a Bayesian framework as described below.The consolidated classification of each galaxy as "merger" or "not merger" turns each observed galaxy into a Bernoulli trial (i.e., "merger" would count as a "success", and "not merger" as "failure").Therefore, the likelihood of our data is given by the discrete binomial probability distribution: where Here, f m is the intrinsic merger fraction for a given sample (or probability of any successive Bernoulli trial being a merger), k is the number of mergers observed for that sample, and n is the number of Bernoulli trials observed from that sample.
We adopted an uninformative uniform proior on merger fraction for each sample, f m ∼ U [0, 1].Combining this prior with the binomial likelihood function given by equation 2, Bayes' theorem yields a posterior distribution for merger fraction of f m ∼ Beta(k+1, n−k+1).For these "shape parameters", the Beta distribution is given by: where The resulting posterior distributions for the merger fractions of our different samples based upon the consensus classification results are shown in Figure 5.A Bayesian two-sample proportion test based upon Monte-Carlo sampling the respective posterior distributions (as implemented in the R package BayesianFirstAid, Bååth 2014) shows that the 3CR sample has a greater intrinsic merger fraction than all control groups at the > 99.9% confidence level.

Monte Carlo Analysis Without Voter Consolidation
In order to strengthen the robustness for our result of an enhanced 3CR merger fraction in comparison to the control groups, we performed an additional analysis following a similar methdology to that of Cisternas et al. (2011) (see their Section 4.2.1),except within our existing Bayesian framework and using the same six human expert classifications described in section 4. In this case, we considered each of the classifier's "votes", or classification choices, individually, allowing each classifier to define their own posterior distribution of merger fraction for each sample.Then, for each control group, we performed a Monte Carlo sampling of the 3CR and control group merger fractions (denoted as fm,3CR and fm,CS , respectively) from their posterior distributions.This was done using 1,000 Monte Carlo samples of fm,3CR and fm,CS per classifier.We then computed the difference in merger fraction for each classifier from these Monte Carlo samples, scaled by each classifer's control sample merger fraction7 : fm,3CR− fm,CS fm,CS .This allows us to include input from all classifiers on how they perceive the 3CR sample relative to the control groups.Normalizing the difference in merger fraction, fm,3CR − fm,CS , by each classifier's control sample merger fractions allows for the combination of different classifiers results while adjusting for the difference in personal scale for each classifier.
Figure 6 shows the results of this analysis for all classifiers after co-adding their histograms.Positive values correspond to an incidence of higher 3CR merger fraction relative to the given control group being compared.Evidence for a statistically significant enhanced merger fraction of the 3CR sample relative to the control groups from this figure would correspond to cases where the confidence interval bounds fall above zero.As is evident, the 95% confidence interval bounds for the 3CR vs V17 and 3CR vs M16 sample comparisons falls above zero, as does the 68% bounds for the 3CR vs M19 sample.The heavy-tailed nature of the 3CR vs M16 histogram reflects a greater variance of control sample merger fractions among the classifers.This analysis supports the finding of an enhanced merger fraction of the 3CR sample relative to all of the control groups.However, the statistical significance of our results do appear stronger using the consensus classification methodology outlined in the previous section, rather than this method of combining the votes from all classifiers.This is consistent with our consensus classification removing some of the human bias inherent to our classifications by our consolidation schema.If this hypothesis is correct, the inclusion of more classifiers should only serve to bring the consensus classifications closer to the "truth" as is implicitly assumed in largescale voter-based morphological classifications in astronomy such as Galaxy Zoo (Lintott et al. 2008).

DISCUSSION
Our finding of high 3CR merger fraction and evidence for enhanced merger fraction of the 3CR sample relative to the radio-quiet AGN control groups, implies that major galaxy mergers are important to the triggering of radio-loud AGN and the launching of powerful radio jets in the cosmic bright ages (1 ≲ z ≲ 2).This result is consistent with that of Chiaberge et al. (2015), who found f m = 0.94 +0.06 −0.16 for the type II z > 1 3CR radio galaxies.The finding of high merger fraction in both the type I and type II z > 1 3CR radio-loud AGN is consistent with AGN unification models which predict AGN type based upon relative orientation with respect to the observer's line of sight (Lawrence & Elvis 1982;Antonucci & Miller 1985;Barthel 1989;Antonucci 1993).Given that not all major galaxy mergers host radio-loud AGN, galaxy mergers are consistent with being an important, although not sufficient, process for the launching of radio jets.
Our result of high merger fraction for the z > 1 3CR radio-loud quasars is also consistent with the work of Podigachoski et al. (2015), who find high star-formation rates for the same parent sample of z > 1 3CR radioloud AGN based upon Herschel far-infrared and Spitzer mid-infrared photometry with corresponding spectral energy distribution decomposition into AGN-heated and star-formation-heated dust components.Podigachoski et al. (2015) interpret these high star-formation rates as likely originating from major gas-rich galaxy mergers, although jet-induced star formation is another possibility and is hard to separate from merger-induced starformation due to the lack of required spatial resolution in the Herschel and Spitzer observations.Follow-up JWST mid-infrared observations would help to resolve the dustenshrouded star-forming regions in the z > 1 3CR sample and help better discern the origin of the prodigious, obscured star formation taking place in these systems.

Potential Impact of Quasar Environment and Jet-Host Interactions
In selecting our control groups, we focused exclusively on redshift and basic quasar properties.The purpose of this selection criteria was to make the presence of a powerful radio jet the major distinguishing variable between the radio-loud and radio-quiet AGN.However, we made no attempt to further control for quasar environment, as is the case in most other quasar merger studies.In principle, the large-scale density of group or cluster members, dark matter halo properties, intergalactic medium (IGM), and various host galaxy properties (e.g., stellar mass, gas/dust density and temparature, color, metallicity, size, etc.) are possible confounding variables to consider when devising a control group and trying to isolate radio-loudness as the most salient parameter influencing the importance of galaxy mergers in triggering AGN.In practice, it is not possible to control for all of these variables simultaneously.Thus, it is still possible radio-quiet quasars which reside in unique environments are preferentially triggered by galaxy mergers.However, it appears major galaxy mergers are an overwhelmingly essential ingredient to the triggering of AGN with powerful radio jets, and not to radio-quiet AGN at cosmic noon.
It is also important to consider the possibility of AGN jet feedback contaminating the merger classifications associated with asymmetric or disturbed host galaxy features.There are many examples in the literature of AGN jets strongly aligned with the semi-major axis of their host galaxies, which has been dubbed the "alignment effect" (e.g., McCarthy et al. 1987;Chambers et al. 1987;Dunlop & Peacock 1993;Best et al. 1998).One explanation for the alignment effect is jet-host interactions as the jet propagates through the ISM and triggers star formation (Rees 1989), thus indicating the possibility of jets altering the appearance of their host galaxies in our 3CR sample.Additionally, it is possibile electron or dust scattering of the AGN optical/UV continuum and the nebular continuum from an extended emission line region may also contribute to the alignment effect (Dickson et al. 1995;Tadhunter et al. 1989;De Young 1998;Tadhunter et al. 2002).We will explore the alignment effect further for the z > 1 3CR sample with our WFC3/UVIS F606W (rest-frame near-UV) HST observations in order to correlate star-formation and radio jet properties in a future study (Breiding et al., in prep.).But, given host asymmetry is only a minor factor in our comprehensive merger criteria and it is unclear jethost interactions have significantly affected the apparent asymmetry in our 3CR host galaxies, we do not believe jet-host interactions should substantially affect our results.Furthermore, jet-host interactions are not expected to affect ongoing mergers or induce tidal signatures which are tell-tale signs of galaxy mergers and pervasive in our 3CR host galaxies (the HST WFC3/IR images and PSF-subtracted images are shown in Appendix section C for our 3CR quasar sample).

A Possible Jet Formation Channel: the Merger of Binary SMBHs
One possible interpretation for our results, as first suggested by Wilson & Colbert (1995) (also see Hughes & Blandford 2003), is that major galaxy mergers can lead to radio jet formation in the so-called "spin paradigm" through the preceding merger of binary SMBHs and resultant spin-up of the remnant SMBH (also see Chiaberge & Marconi 2011;Chiaberge et al. 2015).In this scenario, major galaxy mergers result in the formation, and eventual coalescence, of binary SMBHs of similar mass ratio (Begelman et al. 1980).However, only with the right combination of initial spins and mass ratio will binary SMBH mergers spin-up the remnant black hole.In particular, mass ratios closer to unity and highly spinning black holes aligned with the orbital angular momentum of the binary will lead to the largest spin of the remnant black hole post-coalescence.In the case of gas-rich galaxy mergers, theoretical arguments support a scenario in which the spins of both black holes in the binary will tend to align with its orbital angular momentum through accretion torques mediated by the larger galactic-scale gas flow (Bogdanović et al. 2007).This mechanism is also supported by the observation of three nearby (z < 0.1), young (ages < 10 Kyr) compact symmetric objects showing jet axes normal to the galactic gas disks in near-IR/optical multiband imaging with HST (Perlman et al. 2001).
However, "dry mergers" may more naturally explain the galaxy cores (i.e., missing stellar light) more frequently observed in the centers of radio-loud AGN hosts (e.g., Capetti & Balmaverde 2006, 2007), possibly created through binary or recoiling SMBH scouring in our proposed model of jet formation (Begelman et al. 1980;Milosavljević & Merritt 2001;Nasim et al. 2021).Nevertheless, in the context of spin-powered jets such as those created through the "Blandford-Znajek" mechanism (Blandford & Znajek 1977), the increase of black hole spin following a binary SMBH merger (with the right properties, as described above) would allow for the formation of the radio jet by directly tapping the increased spin energy of the remnant black hole.Interestingly, observations seem to show that only high mass black holes (≳ 10 8 M ⊙ ) tend to launch powerful radio jets (e.g., Chiaberge & Marconi 2011), as supported by the 3CR black hole masses presented in this study (see Figure 2 and Table 2), possibly because the highest-mass black holes can achieve the highest, and most stable, spin magnitudes (Dotti et al. 2013).The "spin paradigm" is further supported by several recent numerical general relativstic magnetohydrodynamic (GRMHD) simulations favoring spin-powered over acretion-powered black hole jets in AGN (e.g., Tchekhovskoy et al. 2010Tchekhovskoy et al. , 2011;;Penna et al. 2013;Narayan et al. 2022).
Since this jet formation channel relies on a binary SMBH merger in order to spin-up the black hole and launch the jet, it makes several key predictions.The first is that the galaxy merger responsible for triggering the radio-loud AGN should be in some stage of postcoalescence, instead of an ongoing or incipient galaxy merger.This follows from the fact that SMBHs take some time to evolve towards the center of galaxies and eventually coaslesce themselves following a galaxy merger (binary evolution is described further in section 6.2.1).While a lot of the host galaxies from our 3CR sample are classified as post-coalescence systems, there are also some number of ongoing galaxy mergers.For these, we would predict a previous galaxy merger is responsible for launching the radio jet in the binary SMBH merger model of jet formation we consider.However, even though post-coalescence merger features can be observed for up to ∼ 1 Gyr (e.g., Lotz et al. 2008;Ji et al. 2014), this scenario is difficult to test since postcoalescence merger signatures are likely contaminated by stronger gravitational disturbances induced by the ongoing galaxy merger.The requirement of a recent previous galaxy merger in the 3CR host galaxies classified as ongoing major mergers can naturally be explained by our radio-loud AGN residing in environments where galaxy mergers are a frequent occurance.Indeed, observations have shown radio-loud AGN reside preferentially in galaxy overdensities (i.e., groups, clusters, and proto-clusters) (e.g., Venemans et al. 2007;Shen et al. 2009;Donoso et al. 2010;Wylezalek et al. 2013;Kotyla et al. 2016;Ghaffari et al. 2017;Retana-Montenegro & Röttgering 2017;Ghaffari et al. 2021), which are environments conducive to frequent galaxy mergers (e.g., Jian et al. 2012).

Binary SMBH Evolution & Recoiling SMBHs
Following the final coalescence of two galaxies in a merger, the SMBHs from each precursor system sink to the center of the merger remnant through dynamical friction (Chandrasekhar 1943), achieving a separation of ∼ 10 pc after a period of ∼ 100 Myr (e.g., Campanelli et al. 2007;Callegari et al. 2009).Subsequent evolution of the binary from ∼ 10 pc to sub-pc separations is highly uncertain, relying on three-body stellar hardening (e.g., Sesana et al. 2007), torqures from a circumbinary accretion disk (e.g., Escala et al. 2005), and possibly three-body interactions with another SMBH (Hoffman & Loeb 2007) to remove angular momentum until the binary enters the gravitational-wave-dominated regime.This evolutionary phase where the binary evolves from ∼ 10 pc to sub-pc separations is referred to as the "final pc problem" (Milosavljević & Merritt 2003), and can take anywhere from 10 Myr (e.g., Khan et al. 2015) to longer than the Hubble time (e.g., Yu 2002) to overcome.Once entering the gravitational-wave regime (separations of ≲ 0.01 − 0.1 pc), the binary rapidly evolves towwards its eventual coalsecence through the emission of low-frequency gravitational waves (this final stage takes ≲ tens of Myr, see e.g., Figure 1  After binary SMBH coalescence, the remnant black hole may experience a recoil velocity-kick due to the anisotropic emission of gravitational waves (Peres 1962;Bekenstein 1973).Depending on the initial spins and mass ratio of the binary, the velocity of the recoiling SMBH can be up to ∼ 5000 km s −1 , although kicks of ≲ a few hundred km s −1 are most common in simulations (Dotti et al. 2010;Lousto et al. 2012;Blecha et al. 2016).Mass ratios near unity, highly spinning black holes, and spins somewhat misaligned with the orbital angular momentum (i.e., "hangup kicks"), or anti-aligned and lying in the orbital plane (i.e., "super kicks") tend to produce the largest recoil velocities (e.g., González et al. 2007;Herrmann et al. 2007;Lousto & Zlochower 2011).
Thus, another potential observable prediction of our binary SMBH merger model is the formation of radio-loud AGN in recoiling SMBH systems.Recoiling SMBHs may manifest themselves observationally through AGN spatially offset from their host galaxy photocenters (e.g., Blecha et al. 2016) or the Doppler shifting of BLR emission lines (e.g., Eracleous et al. 2012, although this can also result from a binary SMBH) as the recoiling SMBH travels through the host galaxy.In fact, the source 3C 186 from our radio-loud AGN sample is possibly the best-evidenced recoiling SMBH in the literature, showing both a substantial 1.3 ′′ (or 11 kpc projected) spatial offset of the quasar from its host center (Chiaberge et al. 2017;Morishita et al. 2022) and ∼ 2000 km s −1 velocity offset between its broad emission lines and host galaxy rest frame (Chiaberge et al. 2017(Chiaberge et al. , 2018;;Castignani et al. 2022).The famous, and nearby (z = 0.004), radio galaxy M87 also shows evidence for a small ∼7 pc projected displacement between the AGN and host galaxy photocenter (Batcheldor et al. 2010), providing some further supportive evidence for this hypothesis.However, systematic searches (with high astrometric precision) for jetted AGN offset from their host galaxy photocenters are needed to more fully examine this hypothesis.
Our Galfit models showed evidence for only one other quasar-host-center offset among the 3CR sample, a 1.4 kpc projected offset in 3C 9 (where we describe the analysis of the 3C 9 offset in detail in Appendix section G).It is possible we do not find any quasarhost-center offsets among the other z > 1 3CR quasars due to some combination of lower recoil velocities 8 and not enough time having elapsed since the recoil event to produce a detectable astrometric offset.However, the very extreme spatial and velocity offset in 3C 186 suggests at least some cases of radio-loud AGN should lead to detectable recoiling SMBHs.It is possible the combination of compact high-z host galaxies and quasar PSF-contamination of the nuclear host galaxy light hide all but the most extreme quasar-host-center offsets 9 .One way to avoid this issue is to look for offset AGN in the type II 3CR radio galaxies where there is no bright quasar PSF cominating the host galaxy light.In this case, the SMBH location could be determined through high astrometric precision, very long baseline interferometry (VLBI) observations of the 3CR radio cores.We have a completed Very Long Baseline Array (VLBA) program to do just this (VLBA/23A-297, legacy ID BB446), where the results will be reported in a future publication.

The Case of the Radio-Loud Non-Merger: 3C 325
It is interesting to note the one radio-loud source with a clearly resolved and featureless host galaxy classified as a "non-merger" by our consensus classification: 3C 325.As shown in Figure 4, 3C 325's host galaxy has the appearance of a smooth and featureless elliptical, as supported by its relatively high Sérsic index of 3.4.One possible formation channel for elliptical galaxies is the merging of massive disk galaxies (Toomre & Toomre 1972), as supported through simulations (e.g., Farouki & Shapiro 1982).If the binary SMBH merger model of jet formation we consider is correct, this might be a system where 8 As described in section 6.2.1,only the right combination of binary spins and mass ratio will lead to large recoil velocities postbinary coalescence.Given that the most spin-up is experienced by remnant SMBHs where the binary spins are aligned with the orbital angular momentum, it is possible most recoil velocities in radio-loud AGN are quite small (see section 6.2 where we discuss favorable binary spins for producing large recoil velocities). 9As is the case in 3C 186.It is interesting to note that 3C 186 also has the largest size of any of our 3CR host galaxies.Perhaps the extended nature of the 3C 186 host galaxy makes a quasarhost-center offset easier to detect.
binary evolution took much longer than the rest of the 3CR sample and any morphological indicators of a past galaxy merger responsible for triggering the radio-loud AGN have long since dissapeared.Better understanding this system (and others like it among radio-loud AGN samples) which appears to be an outlier among the z > 1 3CR sources as having no indications of any ongoing or recent major galaxy mergers will help us better understand the jet-formation physics and AGN-triggering mechanisms operating in other radio-loud AGN.

Alternative Explanations for a High Galaxy Merger
Incidence among Radio-Loud AGN While our preceding discussion focuses on the merger of binary SMBHs as the main channel for the formation of powerful radio jets through the subsequent spinup of SMBHs in the spin-powered jet paradigm, there are other viable hypotheses not excluded by our results.First, it is possible galaxy mergers are still important for the production of powerful jets, but it is gas accretion resulting from a galaxy merger and not the merger of a binary SMBH which acts to spin up a SMBH leading to the formation of a jet through the Blandford-Znajek mechanism (Moderski & Sikora 1996a,b;Moderski et al. 1998;Sikora & Begelman 2013).One meaningful way to observationally discriminate accretion as the mechanism for SMBH spin up as opposed to the merger of binary SMBHs is to search for evidence of recently coalesced SMBHs among radio-loud AGN.As discussed in section 6.2.1, one avenue to explore in this regard is to search for evidence of gravitational-wave-recoiling SMBHs in radio-loud AGN.
Alternatively, it is possible galaxy mergers are not important triggers for radio-loud AGN but rather the findings of high rates of galaxy mergers among radioloud AGN samples are only a byproduct of the fact that radio-loud AGN preferentially reside in overdense group and cluster environments.This scenario presumes that the high radio luminosities of radio-loud AGN are linked to their environments in some way other than through galaxy mergers.One problem for this interpretation is the lack of correlations observed between cluster richness and radio luminosity in radio-loud AGN samples (see the results from the Clusters Around Radio-Loud AGN survey by Wylezalek et al. 2013 andHatch et al. 2014).Another challenge for this interpretation of our results is that while radio-loud AGN are found to statistically favor over-dense environments in comparison to radioquiet AGN, there are still a sizeable fraction of radioloud AGN found in relatively poor environments (Kotyla et al. 2016;Croston et al. 2019).For instance, 3C 298 lacks any evidence for a cluster (e.g., Ghaffari et al. 2021) but is a well known merger with disturbed morphology (Figure 14, and Barthel et al. 2018).This problem can be avoided in the scenario that galaxy mergers trigger radio-loud AGN, as dense environments would only be environments conducive to galaxy mergers but not necessary by themselves to the existence of the radio-loud AGN.One way to test the hypothesis that galaxy mergers are only a byproduct of dense environments and not important triggers for radio-loud AGN would be to control for the AGN environment in addition to the AGN and AGN host galaxy properties in future radio-loud AGN merger studies.

The Assembly, and co-Evolution, of SMBHs and their Host Galaxies
In Figure 7, we plot the host galaxy luminosities 10 for our quasars against their black hole masses (for all quasar host galaxies with reliable Galfit Sérsic model fits as described in section 3.2).The host galaxy luminosities were determined at R-band since it is the closest rest-frame band for most of our quasars (with an effective wavelength of 658 nm, and in νL ν ), using the magnitudes obtained from our Galfit decompositions and K-corrected with the elliptical host galaxy spectral template given in Mannucci et al. (2001) (assuming a simple linear correction redward of the λ4000Å break).We also plot the expected scalings for nearby inactive galaxy samples of "spheroid" 11 luminosity components against well-measured black hole masses from Graham (2007), Kormendy &Ho (2013), andMcConnell &Ma (2013).The Graham (2007) scaling was derived from the Marconi & Hunt (2003) K-band and McLure & Dunlop (2002) B-band relations (assuming R-band color corrections), using elliptical galaxies and the bulges of disk galaxies after correcting for some discrepencies and updating source measurements from the original papers (both of the scalings are essentially the same after the Graham 2007 corrections).The Kormendy & Ho (2013) and McConnell & Ma (2013) scalings are also based upon nearby samples of inactive elliptical galaxies and the bulges of disk galaxies (the latter decomposed from the total host galaxy light), color-corrected from K-band for Kormendy & Ho (2013) and V-band for McConnell 10 We plot luminosities instead of stellar masses since we have no color information for our host galaxies and thus can not obtain reliable mass-to-light ratios.While it is reasonable to substitute mass for host galaxy luminosity when interpreting the power-law scaling exponent, multi-band host color information and careful population synthesis modeling should be undertaken in order to obtain precise host galaxy masses. 11Here we refer to "spheroid" components as either pure elliptical galaxies or the bulge component of spiral or lenticular galaxies.).We distinguish galaxy "mergers" as filled yellow circles and "nonmergers" as filled purple circles as determined from our consensus results.We also mark the reddest quasars with red outlining circles, which have U-V colors in excess of 2. Middle: Same as top, but we distinguish sources based upon quasar sample.Black circles correspond to 3CR quasars, purple V17, and orange represent a combined M19 and M16 high-z sample.Regression lines (with 95% confidence bands) are shown for each sub-sample, with colors corresponding to the sub-sample plotted.Bottom: Same as top panel, but we plot host "bulge" luminosities after adjusting for B/T corrections following the procedure outlined in section 6.3.Linear regression was performed on these adjusted bulge luminosities.2 & 3.For Sérsic Index−log (MBH) correlation, the full sample includes those with measured Sérsic indices.b N is the sample size considered for the correlation.c We give the mean ρs and error bars defining a 68% confidence interval as determined from our bootstrapping procedure.The same is true for the corresponding regression coefficients.d The intrinsic scatter is just taken to be the standard deviation of residuals about the best-fit OLS regression line.& Ma (2013) to R-band again using the elliptical galaxy template of Mannucci et al. (2001).
We performed linear regression analyses on our host galaxy luminosities and black hole masses using an ordinary least squared (OLS) estimator (the best linear unbiased estimator) for our regression coefficients, with our regression equation given by: where α is the intercept, β is the regression slope, ϵ 0 describes the intrinsic scatter, log (L host ) is taken to be the dependant or "response" variable and log M BH /10 9 M ⊙ is taken to be the independant or "predictor" variable.We determined the strength of correlation using the non-parametric Spearman's rankorder correlation coefficient (Spearman 1904), ρ s , which is less sensitive to outliers than Pearson's correlation coefficient (e.g., Wilcox 2004).We also performed hypothesis tests for the log (L host )-log M BH /10 9 M ⊙ regression in order to determine significant postive correlations using a one-tailed Student's t-test (Zar 1972) (where we take p-values less than 0.05 to be cases where we can reject the null hypothesis that ρ s ≤ 0).In order to determine confidence intervals for our regression, we used bootstrapping with 10,000 trials per regression analysis.The results of our regression analyses are given in Table 5.We find a very significant, and tight, correlation for our full sample of quasars (combined from all subsamples considered in this work), with mean ρ s = 0.51 and median p-value of 0.00012.However, we find a much shallower slope of β = 0.25 in comparison to the nearlinear expectation of nearby elliptical galaxies and the bulges of disk galaxies.As show in Figure 7, we also find no obvious trends for galaxy mergers or extremely red quasars in the log (L host ) − log (M BH ) plane.In Figure 7 we also show separate regression analyses after splitting our quasars into the following subsamples: 3CR, V17, and a combined sample of M16 and M19 z ∼ 2 radioquiet sources.As is apparent, the regression slopes are fairly consistent across the samples, although the significance of the corresponding correlations is weakened after splitting the full sample into subsamples (with all but the V17 subsample exhibiting non-significant correlations when examining the median p-values).
A near-linear scaling between black hole mass and host spheroid luminosity (and roughly spheroid mass) can be produced from models in which the host galaxy central spheroid and black hole share a common gas supply regulated by gravitational torques (e.g., Anglés-Alcázar et al. 2017), and also from "merger avergaging" through the hierarchical build-up of both spheroid and black hole mass through galaxy mergers (e.g., Peng 2007;Jahnke & Macciò 2011).However, previous studies have also found evidence for shallower slopes in L spheroid − M BH scalings in the case of late-type galaxy samples ( M spheroid ∝ M 0.33−0.5 BH , e.g., Davis et al. 2019;Savorgnan et al. 2016) and wet (gas-rich) coreless galaxies possibly experiencing "quasar"-mode non-linear growth with respect to their cold gas content (L spheroid ∝ M 0.30−0.45BH , Graham & Scott 2013).Although high-z AGN samples seem to generally have slighly larger black-hole-to-bulge mass ratios (by a factor of ∼ 2 − 3, see Kormendy & Ho 2013), their slopes are broadly consistent with local L spheroid − M BH relations (e.g.Bennert et al. 2011;Schramm & Silverman 2013, although see Laor 2001 who present evidence for the non-linear M spheroid ∝ M 0.65 BH scaling for broad-lined quasars).One other possibility for our shallow slope (L host ∝ M 0.25 BH ) is negative quasar feedback.In this scenario high-z quasars evolve towards the steeper local relation by the suppression of star formation as they grow (effectively moving towards the bottom right in Figure 7).Finally, recent work has suggested that observational bias towards the most massive black holes and combined samples of ellipticals and massive bulges may have yielded near-linear correlations when separate offset, but parallel, non-linear relations (substantially shallower in L spheroid ∝ M BH ) might apply to merger-built elliptical galaxies and the bulges of spiral galaxies (Graham & Sahu 2023, where their Figure 9 best illustrates this potential bias).However, next we discuss the possibility of accounting for this discrepency in slope by our use of the total host galaxy luminosity rather than just the bulge component in the correlation.

MBH−Host Galaxy Correlations
The obvious discrepancy of shallow log (L host ) − log (M BH ) slope (β = 0.25) in comparison to the near-linear slopes found in samples of nearby, inactive spheroids regressed on black hole mass may also in part be explained by our use of total host galaxy luminosity instead of just the spheroid component.Given the high-z nature of our samples and PSF-contamination by our bright quasars, detailed decompositions into both disc and bulge host galaxy components is not possible for our sources.However, if we parametrize the bulge-to-total host galaxy luminosity, B/T, as a powerlaw with index of 3/4 (i.e., B/T ∝ (M BH ) 3/4 ), with B/T = 1 at log (M BH ) = 10, then we obtain the bottom panel in Figure 7 where we plot these adjusted "bulge" values instead of total host galaxy luminosity (we stress that these are not measured bulge values, just those inferred from the aforementioned parametrization).This parametrization yields a consistent linear regression slope with the near-linear slopes from our expected local spheroid relations.However, it predicts a rather extreme dependance of B/T on black hole mass.In our parametrization, B/T > 0.5 (i.e., bulgedominated) values are reached for log (M BH ) > 9.6 and B/T = 0.03 at log (M BH ) = 8.Thus, although the most extreme black hole masses from our samples would correspond to bulge-dominated host galaxies, a majority of our host galaxies would be disk-dominated.This result suggests a scenario where progressive major galaxy mergers lead to more bulge-dominated systems with more massive black holes, where elliptical galaxies would correspond to the most massive black holes and host galaxies with B/T = 1 (e.g., Hopkins et al. 2010).This scaling of B/T with black hole mass is also consistent with the positive correlations observed between Sérsic index12 and black hole mass (e.g., Graham & Driver 2007).We also regressed various other source properties against black hole mass following the same methodology as outlined above (except for these we assess significance using a two-tailed hypothesis test), with the results of these analyses given in Table 5.We found statistically significant correlations (as assessed from median p-values) only for Sérsic index and bolometric luminosity regressed on black hole mass (where the latter correlation has also been observed previously, see e.g., McLure & Dunlop 2002).In Figure 8 we plot Sérsic index against black hole mass, where there is evidently a very strong correlation consistent with the B/T scaling discussed above.We also considered the possibility that the host galaxy luminosities are inflated at low black hole masses due to less well-resolved host galaxies and host galaxy components absorbing some of the quasar flux in our Galfit decompositions.However, we would expect both smaller radii and flux ratios closer to unity for lower black hole masses if this effect were biasing our results, where no such correlations were observed.5.
In summary, we find the following plausible contributing factors for our shallow L host ∝ M 0.25 BH correlation: negative quasar feedback, B/T decreasing (substantially) towards low mass systems, gas-rich systems with highly non-linear black hole growth, or late-type galaxy morphologies which appear to follow a separate trend to strictly early-type galaxy samples.However, we do not know which, or if any, of these factors is responsible for our observed shallow correlation.Better host galaxy decompositions into bulge and disc components with higher-resolution, high-sensitivity near-IR observations with the JWST would help constrain which of these factors is most relevant in our high-z quasar samples, and thus help elucidate the galaxy formation physics and evolution taking place in these systems.Better understanding the physics and origins of different black hole mass scalings with host galaxy luminosity will aid in predictive modeling of the long-wavelength gravitationalwave background from nearby binary SMBHs, which is an imminent prospect in the coming years (Arzoumanian et al. 2020; Chen et al. 2021).

Possible Challenge to the "Blow Out" Paradigm
There have been several recent studies of dustreddened and obscured high-z quasar samples claiming high incidences of major galaxy mergers (Urrutia et al. 2008;Glikman et al. 2015;Fan et al. 2016, although see Zakamska et al. 2019 for a counter-example).It is interesting to note the possibility of radio-loud AGN contaminating these samples and consequently inflating the observed merger fractions, where the sample from We also mark various regions we would expect to find blue, typical colored, or extremely red quasars as discussed in section 6.4.2015) is radio-selected from the FIRST VLA all-sky survey (in conjunction with the 2MASS and UKIDSS infrared surveys).Nevertheless, in the context of some semi-analytic galaxy evolution models where major (gas-rich) galaxy mergers trigger quasar activity by funneling gas to the galactic center, the quasar starts out in an early highly obscured and dust-reddened phase coinciding with early starbursts (e.g., Sanders et al. 1988;Hopkins et al. 2005aHopkins et al. , 2006Hopkins et al. , 2008)).The highly obscured and reddened quasar then continues to grow in mass, and accretion luminosity, until a "blow out" of gas and dust by the quasar halts star formation and uncovers the (no-longer dust-reddened) blue central AGN.This "blow out" paradigm in galaxy and quasar co-evolution gained a lot of popularity as it also predicted the unambiguous association observed between ultraluminous infrared galaxies (ULIRGs) and major galaxy mergers, where in this picture ULIRGs would represent an early highly obscured quasar phase following a major galaxy merger (see the review by Sanders & Mirabel 1996).

Glikman et al. (
As shown in Figure 2, our radio-loud 3CR quasars exhibit similar rest-frame U-V colors to our radio-quiet quasar control samples, where the M16 and M19 samples were constructed from the uniform color-selection algorithm presented in Richards et al. (2002).However, in Figure 9 we also show rest-frame G 13 -I quasar colors (with effective wavelengths of 477 nm and 762.5 nm for G and I-bands, respectively) determined following the same methodology as our rest-frame U-V colors (out-lined in section 2).The sharp peak in the G-I bin [0.5, 1] is consistent with the findings from recent color studies of local quasars with the middle 50% of the SDSS g * −i * color distribution in the same 0.5 − 1 range (Fawcett et al. 2021).Our Glikman et al. (2006) near-IR optical quasar template has a G-I color of ∼ 0.2, but the authors acknowledge the small bias of a bluer continuum from this template than the SDSS DR1 quasar composite spectrum from Vanden Berk et al. (2001).Extremely red quasar samples have been defined in previous studies using observer-frame J-K colors in excess of 1.7 (i.e., Glikman et al. 2004Glikman et al. , 2007)), where observer-frame J-K colors roughly correspond to rest-frame G-I colors for a majority of our 3CR quasars.Only a handful of our 3CR quasars meet this extremely red quasar color selection criteria, with many 3CR quasars consistent with normal or blue quasar colors.However, our G-I distribution appears to have a prominent red skew, consistent with other studies finding more radio detections (many in young, compact radio jets) among red quasar samples (e.g., Klindt et al. 2019;Fawcett et al. 2020;Rosario et al. 2020Rosario et al. , 2021)).
Given that a lot of our 3CR radio-loud quasars are blue (or representative of typical quasar colors) and are observed in recent or ongoing major galaxy mergers, our results present a potential issue for models in which blue quasars are an older evolutionary phase following the "blow out" phase initiated by red quasar feedback.However, this tension can be resolved if the signatures of a major galaxy merger can last long enough for quasars to enter the blue, unobscured phase following the "blow out" event.Previous studies suggest major galaxy merger signatures can last up to ∼ 1 Gyr (e.g., Lotz et al. 2008;Ji et al. 2014), and typical quasar lifetime estimates range from ∼ 10 − 100 Myr (e.g., Martini & Weinberg 2001;Hopkins et al. 2005bHopkins et al. , 2006;;Kelly et al. 2010, although estimates can range up to ∼1 Gyr).Thus, it would appear blue quasars may still be viable phases of quasar evolution following a "blow out" period if no substantial lags occur between the onset of a major galaxy merger and quasar triggering.
As a further test of the "blow out" paradigm of galaxy evolutionary models, we plot our quasar bolometric luminosities against rest-frame U-V quasar colors in Figure 10, where the expectation from "blow out" evolutionary models is that quasar luminosities should reach a peak during the highly reddened blow-out phase (e.g., Hopkins et al. 2006).Performing a linear regression analysis following the same methodology as outlined in section 6.3 (results of which are given in Table 5), we find a signifcant positive correlation between quasar bolometric luminosity and quasar color (i.e., the red-U − V log L bol [erg s −1 ] q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −2 −1 0 Non−mergers or Unresolved Mergers q q Non−mergers or Unresolved Mergers Figure 10.Plot of quasar bolometric luminosity against rest-frame U-V quasar color.Red circles identify those systems classified as galaxy mergers by the consensus classifications (and black circles represent the "non-mergers").The gray line marks the best-fit OLS regression line, with dashed lines and the shaded region outlining the 95% confidence interval for the linear regression based upon our bootsrapped statistics.P-values and regression coefficients are reported in Table 5 for the correlation analysis.
dest quasars are the most luminous), albeit with a very large degree of intrinsic scatter and shallow slope.These findings are consistent with previous works showing red quasars tend to be more luminous than their blue counterparts (e.g., Glikman et al. 2012), and also appear consistent with the "blow out" picture of galaxy evolution.

SUMMARY & CONCLUDING REMARKS
Using human experts, we blindly classified HST WFC3/IR images of 28 type I z > 1 3CR radio-loud quasars, in combination with three separate radio-quiet quasar control samples (∼ 15−20 each), for the purpose of assessing the relevance of major galaxy mergers to AGN radio-loudness.In order to best examine the morphologies of the underlying host galaxies, we subtracted the bright quasar PSFs with Galfit.Our Galfit decompositions included single Sérsic models for the host galaxy components as to not oversubtract the quasar PSF, also allowing us robust host galaxy flux measurements.
Consistent with the results for the type II z > 1 3CR radio galaxies (Chiaberge et al. 2015), we find nearly all of the type I z > 1 3CR quasars are in recent or ongoing major galaxy mergers.Furthermore, we find evidence for a higher merger fraction in our radio-loud quasar sample compared to all three radio-quiet quasar control samples.Higher-resolution, high-sensitivity near-IR observations from the JWST may yield more robust host galaxy detections in these high-z quasar samples in order to better analyze the host galaxy morphologies and obtain more reliable classifications, and thus merger fractions, for all of our samples considered (indeed the JWST has recently detected the host galaxies of z > 6 quasars, see Ding et al. 2022).However, larger sample sizes of both radio-loud and radio-quiet AGN at all redshifts would help further constrain the role of galaxy mergers towards their initial triggering.The LOFAR Two-metre Sky Survey (LoTSS) of the northern sky at 144 MHz (with a sensitivity of 100 µJy beam −1 and angular resolution of ∼ 6 ′′ ) should have yielded the detection of lower-power jets at all redshifts compared to the 3CR sample of radio-loud AGN.This would allow for more leverage to examine any dependance of mergertriggering on jet power in radio-loud AGN (although see Chiaberge et al. 2015 who find no such dependance).To this end, the future Square Kilometer Array (SKA, Dewdney et al. 2009) will allow more sensitive and highresolution low-frequency radio-imaging of the southern sky.
These results support previous works finding a strong connection between radio-loud AGN and major galaxy mergers.We hypothesize that this connection results from a scenario where binary SMBH mergers allow for relativstic jet formation through the spin-up of the remnant SMBH.One prediction from this model is the formation of AGN jets in recoiling SMBHs, where magnitude of the recoil, and relative spin-up of the black hole, will depend on the configuration of binary spins prior to coalescence.These recoiling SMBHs can be observed as spatially offset AGN (as in the case of 3C 186, and possibly 3C 9), or as AGN with velocity-offset broad emission lines.Precise spatial offsets can be identified through follow-up VLBI observations of pc-scale radio cores tracing the SMBH positions (e.g., Breiding et al. 2021).Similarly, velocity-offset broad lines can be identified through careful modeling of their public SDSS spectra (e.g., Eracleous et al. 2012) or dedicated spectroscopic follow-up observations (e.g., Runnoe et al. 2017;Chiaberge et al. 2018).
Using our Galfit-measured host galaxy fluxes, we also examined the scaling between host galaxy luminosity and SMBH mass for our high-z quasar samples, finding a very tight and statistically significant correlation with a much softer slope of β = 0.25 than the nearlinear slopes found in nearby spheroid relations.This much softer relation may result from some combination of a non-linear "quasar"-mode type of SMBH growth, the negative quasar feedback suppression of star formation, predominantly late-type galaxy morphologies, or the strong dependance of bulge-to-total (B/T) galaxy luminosity on black hole mass for our quasar samples.Higher-resolution, high sensitivity near-IR observations with the JWST may help resolve the bulges in our samples, thus allowing for better assessment of the host galaxy morphologies and B/T decompositions.
Finally, major merger activity in our blue radio-loud quasars presents a possible tension with galaxy evolution models in which blue quasars are the end result of red quasar feedback in the "blow out" paradigm.However, this tension can be resolved if galaxy merger signatures outlive typical quasar lifetimes, as supported by recent numerical works.We find support for the notion that red quasars are more luminous than their blue counterparts, finding a statistically significant, but high-scatter, trend between quasar bolometric luminosity and restframe U-V quasar color.This correlation is consistent with the expected peak of quasar luminosity during the "blow out" phase when quasars are significantly dustreddened.
We thank the anonymous referee for the many useful and constructive comments which ultimately helped improve the quality of the manuscript.This research is based on observations made with the NASA/ESA Hubble Space

A. 3C 119
Figure 11 shows the strong point source contamination of the quasar 3C 119.While the contaminating point source is possibly a foreground star, it could not be matched to existing star catalogs beyond Gaia DR3.In fact, both 3C 119 and the unknown point source are Gaia DR3 sources exhibiting proper motions on the order of 1−3 mas yr −1 , which may be indcative of systematic astrometric error introduced by a dual AGN (Hwang et al. 2020).Additionally, both point sources show similar IR and UV fluxes from our WFC3 IR/UVIS imaging.Thus, we can not rule out the possibility of 3C 119 containing a double quasar without follow-up spectroscopy of its companion.Given such a close and bright point source contamination of 3C 119, and the fact that we can not distinguish between a foreground star and a dual quasar, we chose to drop 3C 119 from our merger analysis.

B. PSF-SUBTRACTION UNCERTAINTIES
In order to assess our PSF subtraction uncertainties, we PSF-subtracted a star in the field of the quasar 3C 190, separated from 3C 190 by ∼25 ′′ and with a magnitude of 18.24 (which is similar to many of our quasars in this study).We used the same PSF-creation methodology described in section 3.2.Figure 12 shows the corresponding Galfit decomposition, as well as the average surface brightness as a function of radius in 0.3 ′′ concentric annuli as determined by aperture photometry.There is a slight undersubtraction of the very few center pixels, but our PSF model very closely matches all annuli past 0.3 ′′ .The residual flux is ∼ 8 % of the original flux in the center 0.3 ′′ , and ∼ 20% past this region from 0.3 − 2.4 ′′ .Similarly, the RMS from the residual pixel flux distribution is ∼ 8% of the PSFsubtracted total within 0.3 ′′ , and < 1% from 0.3 − 2.4 ′′ .Thus, our PSF subtractions should allow for reliable recovery of extended host galaxy components beyond the most uncertain central core of the PSF.

D. MODELING TESTS FOR BLENDED GALAXY COMPONENTS
Here we examine the possibility of blended components other than the quasar and its host galaxy impacting our results, given we do not include models for blended components (e.g., very nearby companion galaxies) in our Galfit fitting procedure.We performed tests on two quasars from our 3CR sample, 3C 2 and 3C 245, which have host galaxies with very different types of blended light components.The 3C 2 WFC3/IR HST data reveal two very small components within the overall envelope of the quasar host galaxy, which may be small companion galaxies enveloped by the 3C 2 host galaxy, unrelated foreground galaxies, or possibly even sites of star formation.The 3C 245 WFC3/IR HST data reveal a massive companion galaxy participating in an ongoing major galaxy merger with the 3C 245 quasar host galaxy (as supported by our expert votes).
In this study, we used a single PSF and single Sérsic profile to model the quasar and its host galaxy in Galfit.We also constructed hand-made bad pixel masks for anything which did not correspond to the quasar PSF or its undisturbed host galaxy component (these bad pixel masks correspond to the pixels omitted from the Galfit fitting procedure).For our tests, we performed additional Galfit fits for 3C 2 and 3C 245, using additional Sérsic model components for the blended galaxy components.For the Galfit fits including the blended component models we used the same bad-pixel mask as our original fits, but with the blended components unmasked.We show the bad pixel masks we used for our original method and this blended model component method for 3C 2 and 3C 245 in Figures 15 and 16, respectively.In Figures 15 and 16, we refer to the models using our orginal method of masking the blended components as model 1.Similarly, we refer to the models containing additional Sérsic profiles for the blended components as model 2.
In Figures 15 and 16 we show Galfit decompositions for 3C 2 and 3C 245 using our original fitting procedure and the test fitting procedure based upon additional model components for the blended sources (models 1 and 2, respectively).We used a fixed Sérsic index of four for the 3C 245 host galaxy and its blended companion galaxy (in both model 1 and 2).We also used a fixed   In Section 4, we describe our classification criteria and our choice to reclassify those sources originally classified "unresolved" as "non-mergers" (corresponding to letters F and E, respectively, from the classification choices).Given no massive galaxy companions indicative of an ongoing or incipient major galaxy merger were identified in these cases, this reclassification choice could only bias our results in the event "post-merger" signatures would have been present if the host galaxy were more fully resolved and detected.
Galaxy mergers are a gravitationally violent process where stellar material can become widely distributed as the participating galaxies are coalescing, allowing for stellar material to persist beyond the HST angular resolution limits of our study before final coalescence and dynamical relaxation.Major galaxy mergers can also trigger galaxy-wide starbursts, significantly increasing the stellar mass and luminosity of the host galaxy as a result (e.g., Sanders et al. 1988;Elbaz & Cesarsky 2003;Barnes 2004).Furthermore, gas-rich galaxy mergers in which star formation would be the most dramatic are also shown to exhibit tidal features for longer periods of time (Lotz et al. 2010).Thus, we believe both ongoing and recent major galaxy mergers are preciesly the cases where the host galaxies of our quasars, and associated merger features, are more likely to be detected and resolved.The above arguments not withstanding, below we assess the potential influence of host galaxy non-detections on our results.
Using our human expert classifications, we define an "unresolved" source class to be cases where letter "F", corresponding to "unresolved" host galaxies, was voted for a majority of the time (ties broken towards F).In Table 6, we give our consensus results with the additional column "Unresolved" based upon this criteria.The following column labeled "Unresolved with Massive Companions" indicates the number of objects with massive galaxy companions reclassified as "B", standing for incipient major galaxy mergers.In Appendix section F, we describe how massive galaxy companions meeting the major-merger cutoff are determined for non-detected quasar host galaxies.
The z ∼ 2 high-z M16 and M19 quasar samples have the greatest number of unresolved sources.In part, this can likely be attributed to a combination of cosmological surface brightness dimming (Tolman 1930(Tolman , 1934) and more compact host galaxies (e.g., van Dokkum et al. 2008;Allen et al. 2017) at high redshift.The lower fraction of unresolved sources in the M16 sample compared to the M19 sample may be due to its high black hole mass selection criteria.The M16 sample has much greater black hole masses than the M19 sample.Thus, the well-known scaling between black hole mass and host galaxy stellar mass and luminosity would then predict less massive and less luminous host galaxies in the M19 sample (see Kormendy & Ho 2013 for a thorough review on this scaling relationship), as consistent with less host galaxy detections in the M19 vs M16 AGN samples.Given the greater number of unresolved sources in both of the high-z M16 and M19 control samples in comparison to the 3CR sample, below we examine the resulting potential bias.
As argued in the beginning of this section, we believe the sources classified as "unresolved" are most likely dominated by intrinsically non-merging galaxies.In support of this line of reasoning, we find all six of the M19 resolved sources not fitting the "unresolved" class discussed above are classified as galaxy mergers.However, we can not rule out a recent coalesced major galaxy merger that remains unresolved in these systems.Below we consider the possibility that some fraction of our unresolved sources in the M16 and M19 control samples are actually in these post-merger systems.The mean fraction of votes for "post-merger" (letter D) among merger votes for each galaxy merger from our consensus is ∼ 29%.Assuming all of the "unresolved" sources from Table 6 were actually intrinsically galaxy mergers and not "non-mergers", we should expect ∼ three post-mergers in both the M16 and M19 samples.Since the postcoalescence galaxy mergers are the only types of galaxy merger we should miss in systems with unresolved host galaxies, this assumption corresponds to increasing the number of galaxy mergers for our consensus classifications found in the M16 and M19 samples by ∼ three (and correspondingly three less non-mergers).This would correspond to a merger fraction of f m = 0.59 ± 0.11 for the M19 sample and f m = 0.45 ± 0.11 for the M16 sam-ple 14 .Both of these merger fractions are still well below that found for the 3CR sample, confirming our result of enhanced merger fraction for the radio-loud AGN in comparison to the radio-quiet AGN.As argued in the beginning of this section, recent major galaxy mergers representative of these post-merger systems should actually be more likely to have resolved and detected host galaxies, so we believe this potential source of bias to be insignificant and the merger fractions presented in Table 4 are closer to the true values.Furthermore, largescale tidal signatures indicative of post-merger systems should still be detectable well outside of the angular resolution limits of the HST observations presented in this study, and the most likely case is that the quasars with unresolved host galaxies from our study have not undergone a recent major galaxy merger.We used a 1:4 flux ratio throughout this work in order to distinguish between an ongoing or incipient major or minor galaxy merger.However, this requires using our Galfit-modeled fluxes of the quasar host galaxies for robust estimates.For cases in which we could not obtain reliable Galfit fluxes of the quasar host galaxy (as described in section 3.2), we estimate the likely quasar host galaxy flux using the empirical scaling relation shown in Figure 7 for the full set of quasars combining all subsamples.These fluxes are then compared to any potential companion galaxies to determine if they are incipient major or minor galaxy mergers (as described in section 4).
G. THE OFFSET QUASAR 3C 9 Figure 17 shows the HST WFC3/IR image of 3C 9 along with the PSF-subtracted image of its host galaxy.Our Galfit decompositions yielded a best-fit quasar PSF model offset from the the Sérsic component by ∼ 0.06 ′′ in the North-Eastern direction (with a 1σ uncertainty of only 6 mas reported by Galfit), or ∼ 0.5 kpc 14 Considering the standard error, we would require a ∼ 10σ deviation, or at least 80%, in post-merger fraction among the M19 "unresolved" post-coalescence galaxy mergers before the resulting merger fraction was more consistent with that obtained for the 3CR sample (again assuming all of the unresolved sources are actually in galaxy mergers).Given the mean fraction of postmerger "D" votes among merger votes for the M19 galaxy mergers is 0.27, and thus consistent with that found for the rest of our samples, this scenario is highly unlikely.However, even if all of the M16 "unresolved" sources were actually post-mergers, the resulting merger fraction would still be much lower than that found for the 3CR sample.at the redshift of 3C 9.In order to confirm this offset, we performed isophote fits of the PSF-subtracted host galaxy image using the photutils (Bradley et al. 2022) python package (the isophote-fitting alogithrm follows the methodology of Jedrzejewski 1987).We used 18 isophotes ranging from ∼ 1.3 − 1.8 ′′ in half-pixel increments of the semi-major axis in order to avoid the PSF-subtraction uncertainties which contaminate the center of the host galaxy.Our best-fit isophotes are shown in Figure 17, where we find an 0.17 ± 0.03 ′′ , or 1.4 ± 0.25 kpc projected, quasar-host-center offset also in the North-Eastarn direction.This roughly 5.6σ offset is suggestive of a possible recoiling SMBH system, but it is also possible the host asymmetry resulting from a recent major merger is enough to account for this offset.One way to test the recoil hypothesis in this source is to search for velocity-offset broad emission lines with follow-up spectroscopy.

Figure 1 .
Figure 1.Redshift histogram distributions for our samples.The shaded black histogram represents the 3CR sample, orange represents the M19, purple the V17, and cyan the M16 sample.See section 2 for a description of each of these control samples.

Figure 2 .
Figure 2. Histogram comparisons of the black hole mass, bolometric luminosity, Eddington ratio, and rest-frame U-V quasar color of the control samples to the radio-loud 3CR sample.The shaded black histograms represent the 3CR sample, orange represent the M19, purple the V17, and cyan the M16 sample.

Figure 3 .
Figure 3. Galfit decomposition panels used for classifications, where we show examples corresponding to each classification category as supported by the expert voting results.Images are constructed in the original detector frame, with varying right ascension/declination axis orientations.The classification is given in the upper left of each panel (as supported by the voting results), and source/sample are given in the upper right.All panel decompositions can be found as online-only figures.The green circles have radii corresponding to 25 kpc projected onto the plane of the sky for each quasar redshift.
20 † The error bars represent 95% credible intervals around the mean values.These are Bayesian merger fractions determined from the posterior distrbutions for each sample's consensus classification.

Figure 5 .
Figure 5. Posterior probability density function (PDF) of the merger fractions for each sample, based upon the consensus classification results given in Table4.Vertical dashed lines with semi-transparent bands show the means of these distributions with 68% credible intervals.

Figure 6 .
Figure 6.Co-added histograms of each classifer's Monte-Carlo sampled difference in merger fraction of the 3CR sample relative to each control group, scaled by the control sample merger fraction (again defined for each classifier).The 68% confidence interval is shown as dashed line, and the 95% confidence interval as the dot-dashed line.Positive values indicate an enhnced 3CR merger fraction relative to the given control group.

Figure 7 .
Figure 7. Top: Scatter plot of R-band host galaxy luminosities and black hole masses.Our best-fit mean linear regression line is given in dark gray, with 95% confidence band given in gray with outlining gray dashed lines.Expected correlation between the host spheroid component of nearby inactive galaxies with black hole mass is shown as solid colored lines (purple for McConnell & Ma 2013, cyan for Kormendy & Ho 2013, and green forGraham 2007).We distinguish galaxy "mergers" as filled yellow circles and "nonmergers" as filled purple circles as determined from our consensus results.We also mark the reddest quasars with red outlining circles, which have U-V colors in excess of 2. Middle: Same as top, but we distinguish sources based upon quasar sample.Black circles correspond to 3CR quasars, purple V17, and orange represent a combined M19 and M16 high-z sample.Regression lines (with 95% confidence bands) are shown for each sub-sample, with colors corresponding to the sub-sample plotted.Bottom: Same as top panel, but we plot host "bulge" luminosities after adjusting for B/T corrections following the procedure outlined in section 6.3.Linear regression was performed on these adjusted bulge luminosities.

e
Ratio of host to quasar flux density in the WFC3/IR HST filter used for each quasar.f We report median p-values from our bootstrapped statistics, with error bars corresponding to the 68% confidence intervals (confidence bands in Figures 7 & 8 are constructed using 95% intervals).P-values for log (L host ) − log (MBH) correlations are constructed from one-tailed hypothesis tests, and all others are two-tailed.g R eff are the effective radii in kpc.

Figure 8 .
Figure 8.Here we plot the Sérsic index against black hole mass as filled circles, with mean linear regression lines shown in dark gray and 95% confidence bands shown with outlined dashed gray lines (constructed from the bootstrapped statistics).A horizontal dashed black line is shown at the mean Sérsic index of 1.94.Regression results for this plot are given in Table5.

Figure 9 .
Figure9.Histogram of rest-frame G-I quasar colors for our 3CR radio-loud AGN sample.We also mark various regions we would expect to find blue, typical colored, or extremely red quasars as discussed in section 6.4.

Figure 11 .
Figure 11.HST WFC3/IR image of 3C 119, logarithmically scaled.The green circle centered on 3C 119 has a 25 kpc projected radius at the source redshift of z = 1.02.
Figures 13 and 14 show the HST WFC3/IR images and PSF-subtracted images of the z > 1 3CR quasar sample.

Figure 12 .
Figure12.Here we show an example PSF subtraction of a star within the field of 3C 190. Green circles mark concentric 0.3 ′′ spaced annuli used for the aperture photometry presented in the bottom plots.In the bottom, we plot average surface brightness and its logarithm for each aperture defined by the concentric annuli shown in the above figures.In black we plot the results for the image data, in gray our PSF model, and in green the PSF-subtraction results.The surface brightness is normalized to the value found in the central 0.3 ′′ circular aperture for the original image, corresponding to 9.0 × 10 −27 erg s −1 cm −2 Hz −1 arcsec −2 .

Figure 13 .
Figure 13.HST WFC3/IR images of the z > 1 broad-lined 3CR quasars.These images are the same as those given in the panel decompositions used for merger classifications.Green circles are centered on the quasars, with 25 kpc projected radii.Images are constructed in the original detector frame, with varying right ascension/declination axis orientations.

F
. CRITERIA FOR DETERMINING HOST-TO-COMPANION FLUX RATIOS IN UNRESOLVED SOURCES

Figure 17 .
Figure 17.Image panel illustrating the quasar-host galaxy spatial offset in 3C 9, with a logarithmic flux scale for all four images.(a) HST WFC3/IR image of 3C 9.The compass gives the orientation of both right ascension and declination axes.The horizontal arrow represents a 3 ′′ scale bar.The green cross marks the host center as determined from the isophote fits shown in panel c, and the cyan circle marks the quasar best-fit position determined by Galfit.(b) the same as panel a but for the PSF-subtracted image obtained from our Galfit decomposition.(c) HST WFC3/IR PSFsubtracted image of 3C 9. Isophotes used to determine the host galaxy center are shown in white.(d) same as panel c but instead of showing isophotes used during offset analysis we show the corresponding ellipsoid model subtracted from the host galaxy light.

Table 1 .
HST WFC3/IR Observations The exposures are given in the format of number of dithers × single frame exposure time.

Table 5 .
Results of Linear Regression Analyses Full sample for log (L host ) − log (MBH) correlations (and others correlations including host galaxy properties) is the combined sample of 3CR, V17, M16, and M19 quasars with host galaxy luminosities as determined from our Galfit decompositions, and measured black hole masses as given in Tables

Table 6 .
Consensus Results & Unresolved Sources