# A Close Quasar Pair in a Disk-Disk Galaxy Merger at $z=2.17$

Yu-Ching Chen1,2 Xin Liu1,2 Adi Foord3 Yue Shen1,2 Nianyi Chen4 Miguel
Holgado4 Tiziana Di Matteo4,5,6 Masamune Oguri7,8 Hsiang-Chih Hwang9 Nadia
Zakamska10

###### Abstract

Most local massive galaxies, if not all, are believed to harbor a supermassive
black hole (SMBH) at the center[1]. Galaxy mergers have long been thought to
drive strong gas inflows[2] and accretion onto either one or both central
SMBH, triggering single or dual quasars[3] as a natural stage of the
hierarchical galaxy and SMBH evolution. While many dual active galactic nuclei
– the low-luminosity counterparts of quasars – have been observed in low-
redshift mergers[4], no unambiguous dual quasar is known at cosmic noon
($z\gtrsim 2$) when both quasar activity and global star formation density
peaked[5, 6]. While a handful of dual quasar candidates were known at $z>1$[7,
8, 9, 10, 11, 12, 13, 14, 15], they still had competing explanations, either
as pre-merger projected quasar pairs at larger physical separations, or as
double images of single quasars produced by the strong gravitational lensing
by foreground galaxies. Here we report multi-wavelength observations of SDSS
J074922.97$+$225511.8 (SDSS J0749$+$2255 for short) as the first kiloparsec
(kpc)-scale dual quasar confirmed to be hosted by a galaxy merger at cosmic
noon. Hubble Space Telescope (HST) near-infrared (NIR) imaging reveals
extended host galaxies underlying the much brighter compact quasar nuclei
(with a separation of 0$.\\!\\!^{\prime\prime}$46 corresponding to 3.8 kpc at
$z=2.17$) and low surface brightness tidal features as conclusive evidence for
ongoing galactic interactions. The detection of both extended host galaxies
and tidal features unambiguously establishes SDSS J0749$+$2255 as a pair of
quasars hosted by an ongoing galaxy merger. We also present Keck NIR imaging
assisted by adaptive optics that better resolves the host galaxy structure,
spatially resolved HST optical and Gemini optical and NIR spectroscopy,
Chandra X-ray imaging spectroscopy, and NSF’s Karl Jansky Very Large Array
radio dual-band imaging in its most extended configuration, all lending
support to the dual quasar hypothesis. Unlike their low-redshift low-
luminosity counterparts, the high-redshift dual quasar is hosted by two
massive (with stellar masses of ${\sim}10^{11.5}M_{\odot}$) compact (with
effective radii of 1.8 kpc and 1.1 kpc) disk-dominated galaxies. The disk
dominance may be critical for efficient gas fueling onto the central SMBHs
producing two luminous quasars in the early-stage merger. The apparent lack of
stellar bulges and that SDSS J0749$+$2255 already follows the local SMBH
mass–host stellar mass relation are at odds with the canonical SMBH-host co-
evolution picture, suggesting that at least some SMBHs may have formed before
their host stellar bulges. While still at kpc-scale separations where the
host-galaxy gravitational potential dominates, the SMBHs may evolve into a
gravitationally bound binary system in $\sim$0.22 Gyr. The merger products at
low redshift are expected to be gravitational wave sources for pulsar-timing
arrays[16].

Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana,
IL 61801, USA

National Center for Supercomputing Applications, University of Illinois at
Urbana-Champaign, 605 East Springfield Avenue, Champaign, IL 61820, USA

Kavli Institute of Particle Astrophysics and Cosmology, Stanford University,
Stanford, CA 94305, USA

McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213, USA

NSF AI Planning Institute for Physics of the Future, Carnegie Mellon
University, Pittsburgh, PA 15213, USA

OzGrav-Melbourne, Australian Research Council Centre of Excellence for
Gravitational Wave Discovery, Melbourne, VIC 3122, Australia

Center for Frontier Science, Chiba University, Chiba 263-8522, Japan

Department of Physics, Graduate School of Science, Chiba University, Chiba
263-8522, Japan

School of Natural Sciences, Institute for Advanced Study, Princeton, 1
Einstein Drive, NJ 08540, USA

Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA

SDSS J0749$+$2255 is an optically selected type 1 (i.e., broad-line) quasar at
$z=2.17$ spectroscopically confirmed by the SDSS legacy survey[17]. It was
identified as a candidate to host a kpc-scale SMBH pair using a new
Varstrometry technique based on its large Gaia astrometric excess noise, which
may be induced by the asynchronous variability of an unresolved double
quasar[18]. The kpc-scale separation represents a physical regime close enough
for merger-induced effects to be important[3] yet wide enough to be directly
resolvable with current facilities. HST follow-up optical dual-band imaging
has revealed compact double nuclei with a projected nuclear angular separation
of 0$.\\!\\!^{\prime\prime}$46[12]. Close compact optical double nuclei could
represent a dual quasar, a projected quasar pair in two non-interacting
galaxies, a lensed quasar, or a star-quasar superposition. Previous Very Long
Baseline Array (VLBA) 9.0 GHz imaging detected strong radio emissions from
both nuclei, ruling out star-quasar superposition[12]. While morphological
differences in the VLBA 9.0 GHz image of the two radio nuclei favor a dual
quasar or a pre-merger projected quasar pair at larger physical separations,
lensing could not be conclusively ruled out[12].

To unambiguously confirm the nature of SDSS J0749$+$2255, we have conducted a
comprehensive multi-wavelength follow-up campaign. The new observations
include HST IR F160W ($H$ band) imaging, Keck adaptive optics (AO) assisted IR
($K_{p}$ band) imaging, HST/STIS optical spatially resolved slit spectroscopy,
Gemini spatially resolved optical (GMOS) and NIR (GNIRS) slit spectroscopy,
Chandra ACIS-S X-ray 0.5–8 keV imaging spectroscopy, and Very Large Array
(VLA) A-configuration (A-config) radio dual-band imaging.

First, after we subtract the bright central quasars modeled by two point-
spread functions (PSFs), the deep HST IR image uncovers two extended host
galaxies underlying both compact quasar nuclei (Figure 1). In addition, after
we further subtract the best axisymmetric models for the extended host
galaxies, the final residual image reveals low surface-brightness tidal
features indicative of ongoing galactic mergers interactions (Figure 1).
Detecting both the tidal features and the extended host galaxies, each of
which coincides with each compact quasar nucleus, unambiguously confirms SDSS
J0749$+$2255 as a dual quasar in an ongoing merger rather than a lensed quasar
or a pre-merger projected quasar pair at larger physical separations.
Complementary higher angular-resolution Keck AO-assisted IR imaging
independently supports the detection of the extended host galaxies and the
non-detection of a central foreground lens galaxy. The alternative scenario of
a foreground lens galaxy as a single extended source located in between two
compact quasar nuclei is quantitatively ruled out by lensing mass model tests
for the Keck AO image. However, the tidal features seen in the HST image are
too faint to be detected with the shallower Keck AO image (Methods).

Second, Chandra ACIS-S detects both quasar nuclei in the hard X-rays (2–8 keV)
as compact X-ray point sources whose spatial profiles are consistent with the
PSF (Figure 2). The X-ray positional offsets between the two nuclei are
consistent with those of the HST optical nuclei. The optically fainter NE
nucleus is brighter in the X-rays than the SW nucleus. This contradicts naive
expectations from a lensed quasar, since gravitational lensing is in general
achromatic: the deflection angle of a light ray (and the lensing
magnification) does not depend on its wavelength so that consistent flux
ratios are expected at different wavelengths from different lensed images
given the same underlying source spectrum. While the wavelength-dependent
geometry of the different emission regions may result in chromatic effects,
the differences in the flux ratios between different images at different
wavelengths should be minor in a lensed quasar. However, the 2–8 keV X-ray
flux ratio between the two nuclei is 45${}^{+41}_{-33}$, which is
significantly different from that observed by HST in the optical
(0.229$\pm$0.001 in F475W and 0.284$\pm$0.002 in F814W), disfavoring lensing.
The X-ray-to-optical luminosity ratio typically characterized by the spectral
slopes $\alpha_{{\rm OX}}$ are 1.46$\pm$0.01 for the NE nucleus and
2.32$\pm$0.13 for the SW nucleus. The estimated unabsorbed hard X-ray (2–8
keV) luminosities (assuming a redshift of z = 2.17 for both X-ray point
sources) are $\sim$10${}^{45.29^{+0.05}_{-0.03}}$ erg s-1 for the NE nucleus
and $\sim$10${}^{43.62^{+0.32}_{-0.58}}$ erg s-1 for the SW nucleus, which
greatly exceed the most X-ray luminous starburst galaxies[19].

Third, both quasar nuclei are detected as compact radio sources by VLA
A-config continuum imaging at 6 GHz and 15 GHz (Figure 3). The radio point
source centers are consistent with those of the HST optical quasar nuclei. The
radio flux ratios between the two nuclei are 0.379$\pm$0.001 at 6 GHz and
0.351$\pm$0.003 at 15 GHz, which are different (by $>5\sigma$) than those
observed by HST in the optical (0.229$\pm$0.001 in F475W and 0.284$\pm$0.002
in F814W), disfavoring lensing. The radio spectral indices are
$-$0.397$\pm$0.004 for the NE nucleus and $-$0.314$\pm$0.002 for the SW
nucleus. The radio loudness parameter $R_{6\,cm/2500\,\AA}$, commonly defined
as the flux density ratio at the rest-frame 6 cm and that at 2500 Å, is
1079$\pm$74 for the NE nucleus and 659$\pm$16 for the SW nucleus, placing both
quasars in the radio-loud (i.e., $R_{6\,cm/2500\,\AA}>$10) population
according to the canonical definition[20].

Finally, HST spatially resolved optical spectroscopy shows significantly
different rest-frame UV continuum spectral slopes between the two nuclei
($\alpha_{\nu}$ = $-$0.70$\pm$0.24 for the NE nucleus and 0.08$\pm$0.06 for
the SW nucleus), independently confirmed by Gemini/GMOS spatially resolved
optical slit spectroscopy (Figure 4) despite some discrepancy between the
STIS-based and the GMOS-based estimates since GMOS only marginally separated
the two nuclei. The different spectral slopes disfavor a lensed quasar,
although gravitational lenses may also produce different spectral slopes due
to varying reddening along different light paths through the foreground lens
galaxy[21]. There are noticeable differences in the broad emission lines, but
the low spectroscopic data quality cannot conclusively rule out lensing based
on the line shapes alone. The best-fit systemic redshifts estimated from the
optical-IR spectra of the two quasar nuclei suggest a line-of-sight velocity
offset of $\sim$310$\pm$50 km s-1 (Extended Data Table 1), consistent with
expectations of typical orbital velocities in galaxy mergers[22] but
inconsistent with lensing (which would produce identical redshifts of the two
quasars).

In summary, our new multi-wavelength follow-up observations unambiguously
establish SDSS J0749$+$2255 as a dual quasar hosted by an ongoing galaxy
merger at $z=2.17$. By modeling the PSF-subtracted host-galaxy-only two-band
(HST F160W and Keck AO $K_{p}$) photometry of SDSS J0749$+$2255, we estimate
the stellar masses of the two host galaxies as
${\sim}10^{11.46^{+0.02}_{-0.26}}M_{\odot}$ and
${\sim}10^{11.50^{+0.04}_{-0.27}}M_{\odot}$ (with the errors accounting for
the 1$\sigma$ statistical and partial systematic uncertainties from surface
brightness profile fitting but ignoring additional systematics from population
synthesis modeling) for the NE nucleus and the SW nucleus. The massive hosts
appear to be compact and disk-dominated (with best-fit Sérsic indices close to
unity; Methods), typical of galaxies[23] and hosts of single quasars[24] at
$z>2$. The disk dominance of the hosts in SDSS J0749$+$2255 is dramatically
different from the more bulge-dominated characteristics of low-redshift low-
luminosity dual AGN hosts[25]. The latter may be more subject to the
stabilization effect of stellar bulges preventing more efficient SMBH fueling
in the early merging stages, resulting in low-luminosity AGNs rather than the
luminous quasars observed in SDSS J0749$+$2255.

By jointly modeling the optical-NIR spectrum (Methods), we estimate the SMBH
masses as ${\sim}10^{9.1\pm 0.4}M_{\odot}$ and ${\sim}10^{9.2\pm
0.4}M_{\odot}$ (with the errors representing 1$\sigma$ total uncertainties)
for the NE and the SW nuclei. Bolometric luminosities estimated from the
optical continua suggest Eddington ratios of $\sim$0.12$\pm$0.03 and
$\sim$0.25$\pm$0.07 for the NE and SW nuclei, respectively. The SMBH masses
and host total stellar masses are consistent with the empirical scaling
relations observed in both single quasars at $z\sim 2$ and local inactive
galaxies within uncertainties (Methods).

SDSS J0749$+$2255 represents the first robust example of a kpc-scale pair of
merging SMBHs in a galaxy merger at cosmic noon ($z\gtrsim 2$), when both the
space density of quasars and the global star formation rate density of
galaxies peaked[5, 6]. It fills a long-standing gap in the expected population
of dual quasars hosted by high-redshift ongoing galaxy mergers. Such systems
are expected to be common in hierarchical merger models, which have been
otherwise successful in explaining much of the quasar phenomenology[3]. In
comparison, the previous redshift record holder of kpc-scale dual quasars
hosted by confirmed galaxy mergers was LBQS 0103$-$2753 at $z=0.86$, which was
serendipitously discovered[26].

A handful of small scale (sub-10 kpc) dual quasar candidates were known at
$z>1$[7, 8, 9, 10, 11, 12, 13, 14, 15]. However, none of them has been
confirmed to be hosted by ongoing galaxy mergers – as expected in bona fide
dual quasars (Methods). The establishment of SDSS J0749$+$2255 as a dual
quasar hosted in a galaxy merger conclusively demonstrates that at least some
high-redshift kpc-scale dual quasar candidates may indeed be dual quasars
rather than gravitational lenses or pre-merger projected quasar pairs at large
separations. On larger (i.e., $>10$ kpc) scales, it has long been suggested
that most wide separation (3′′–10′′) double quasars are quasar pairs rather
than gravitational lenses based on statistical arguments[27]. There are
$\sim$tens of known “binary” quasars at $z>1$ with projected separations of
$>10$ kpc[28] and even a quasar “quartet”[29]. However, the separations of
most of these systems are still too large ($\gtrsim$ tens of kpc) for the host
galaxies of the quasar pairs to be in direct galactic interactions. Unlike
SDSS J0749$+$2255, none of them has been observed with tidal features as
direct evidence for ongoing interactions.

Despite significant merit and intense theoretical interest, direct evidence
for compact binary SMBHs – SMBHs that are gravitationally bound with each
other – is still lacking[4]. Their separations (typically less than a few
parsecs) are too small to resolve beyond the local universe with current
facilities. High-redshift kpc-scale dual quasars represent promising
progenitors of low-redshift compact binary SMBHs. Using a stellar dynamical
friction argument, we estimate that the 3.8 kpc SMBH pair in SDSS J0749$+$2255
may form a compact bound SMBH binary in $\lesssim$0.22${}^{+0.50}_{-0.16}$ Gyr
(Methods). Future infrared integral field unit spectroscopy with the James
Webb Space Telescope will map the stellar and gas kinematics, enabling better
predictions on the subsequent SMBH pair orbital evolution and the impact of
dual quasar on their host galaxies.

The confirmation of SDSS J0749$+$2255 as a dual quasar provides a proof-of-
concept for combining Gaia and HST to discover high-redshift kpc-scale dual
quasars[12]. The discovery underscores the need for high-resolution deep NIR
imaging, which is crucial for detecting the merging galaxy hosts and for
robustly distinguishing high-redshift dual quasars from lensed quasars and
pre-merger projected quasar pairs at larger physical separations. Future high-
resolution wide-area deep surveys[30] will uncover a larger population of
high-redshift dual quasars at the peak epoch of galaxy and SMBH assembly.

##  0.1 Target Selection

SDSS J0749$+$2255 is one of the dual quasar candidates selected from
spectroscopically confirmed quasars in the SDSS DR7 and DR14 quasar
catalogs[17, 31, 32] using the new astrometric technique dubbed
“Varstrometry”[18, 12]. Varstrometry utilizes the centroid jitters due to the
asynchronous variability of an unresolved dual quasar. Dozens of high-redshift
kpc-scale dual quasar candidates were selected by applying Varstrometry to
Gaia DR2 data and follow-up HST optical dual-bang imaging[12, 33]. SDSS
J0749$+$2255 is a high-piority target for detailed follow-up observations
given the high redshift of $z=2.17$ and the small separation of
0$.\\!\\!^{\prime\prime}$46\. Extended Data Table 1 lists its basic
photometric and spectroscopic properties. Previous high angular-resolution
VLBA observations showed two compact radio cores, confirming both sources as
quasars[12]. However, the alternative scenarios of a lensed quasar or a pre-
merger projected quasar pair at larger physical separations could not be
completely ruled out based on the previous HST dual-band optical and VLBA
single-band radio images.

##  0.2 Comparison with Known Dual Quasar Candidates

Extended Data Figure 1 shows SDSS J0749$+$2255 in comparison with all the
previously known $1<z\lesssim 3$ small-scale (sub-10 kpc) candidate dual
quasars from the literature. Four of the candidates[9, 10, 11] are nearly
identical quasars (NIQs; shown in grey) – quasars with nearly identical
optical spectra, which are generally thought to be lensed quasars lacking a
definitive detection of foreground lenses due to observational
limitations[11]. Three have more noticeable spectral differences[7, 8, 13]
(shown in magenta), but without the host galaxy mergers unambiguously
observed, they are still debatable as either pre-merger projected quasar pairs
at larger physical separations or gravitational lenses. The dual quasar
candidate[12] (shown in blue) exhibits differences in the broad emission line
strengths[34] but is still controversial as a quasar lens or a pre-merger
projected quasar pair at larger physical separations due to the lack of high-
resolution deep NIR imaging. Finally, a $z=2.4$ quadruply lensed quasar (shown
in green) has been proposed to contain a lensed kpc-scale dual quasar[15], but
the interpretation is subject to lense modeling uncertainties.

##  0.3 HST/WFC3 NIR Imaging

High-resolution NIR imaging of SDSS J0749$+$2255 was conducted on 2 February
2022 UT (Program GO 16892; PI: X. Liu) with the HST Wide Field Camera 3 (WFC3)
in the F160W filter in a single orbit. The F160W filter operates in the wide
$H$ band with a pivot wavelength of 1536.9 nm and a width of 268.3 nm[35]. The
total net exposure time was 2055 s. The IR detector has a pixel size of
0$.\\!\\!^{\prime\prime}$13, which undersamples the PSF with a full width at
half maximum (FWHM) of 0$.\\!\\!^{\prime\prime}$15\. We dithered the
observations to properly sample the PSFs while rejecting cosmic rays and bad
pixels. We reduced the HST data following standard procedures[36]. Geometric
distortion and pixel area effects were corrected. We combined dithered frames
while rejecting cosmic rays and hot pixels. The final image has a pixel scale
of 0$.\\!\\!^{\prime\prime}$06 for the $H$-band image to be Nyquist-sampled.
The photometric zeropoint is 25.936 (28.177) in the AB (ST) system, which was
derived based on observations of four white dwarf standard stars[35]. Figure 1
shows the HST F160W image of SDSS J0749$+$2255.

##  0.4 Keck AO-Assisted NIR Imaging

Keck AO-assisted NIR imaging was conducted on 21 December 2021 UT (Program
N072; PI: Y. Shen) with the near-infrared camera NIRC2 in the $K_{p}$ filter.
The $K_{p}$ filter has a central wavelength of 2.124 um and a bandpass width
of 0.351 um. The observations were conducted in the Laser Guide Star (LGS)
mode. The artificial LGS was used as a reference object to correct for
atmospheric distortion. We used the star N85C000573 (located
52$.\\!\\!^{\prime\prime}$62 from SDSS J0749$+$2255) as the reference tip-tilt
star. The images were dithered using a 4${}^{\prime\prime}\times$4′′
9-position box. The total net exposure time was 1350 s. The individual images
were reduced following standard procedures and retrieved through the Keck
Observatory Archive. Individual exposures are then co-added together after
correcting for the dither offsets. The final image has a pixel scale of
0$.\\!\\!^{\prime\prime}$01\. The zeropoint of the Keck $K_{p}$-band image is
24.80 (Vega), which was derived based on measurements taken on AO engineering
stars during laser checkouts. Extended Data Figure 2 shows the Keck image of
SDSS J0749$+$2255.

##  0.5 Host Galaxy Morphology, Structural Decomposition, and Lensing Mass
Model Test

We study host galaxy morphology and perform PSF-host structure decomposition
for the HST F160W (H-band) and Keck AO $K_{p}$ images. We adopt GALFIT[37] to
model the 2D surface brightness profiles for the galaxy merger. A field star
is adopted to build the PSF model for the HST image. For the Keck image, a
standard star, S852-C, which was observed during the same night closest to
SDSS J0749$+$2255, is adopted to build the PSF model. Since the tip-tilt
reference star is much fainter than the standard star, the strehl ratio of
SDSS J0749$+$2255’s Keck image is much lower. Therefore, the standard star is
convolved with a wider Gaussian to match the PSF size of SDSS J0749$+$2255’s
observation (with a FWHM of 0$.\\!\\!^{\prime\prime}$05) to correct for this
mismatch.

The baseline model adopts two PSFs (for the two unresolved quasars), two
Sérsic models (for the two extended host galaxies), and a constant background
to fit the image. The Sérsic model is given by

$\Sigma(r)=\Sigma_{e}\,{\rm
exp}\bigg{[}-\kappa\bigg{(}\Big{(}\frac{r}{r_{e}}\Big{)}^{1/n}-1\bigg{)}\bigg{]},$
(1)

where $\Sigma(r)$ is the pixel surface brightness at radial distance $r$,
$\Sigma_{e}$ is the pixel surface brightness at the effective radius $r_{e}$,
and $\kappa$ is a parameter related to the Sérsic index $n$. $n=1$ for an
exponential profile, whereas $n=4$ for a de Vaucouleurs profile. Bulge-
dominated galaxies have high $n$ values ($n>2$), while disk-dominated galaxies
have $n$ close to unity. We try different models until reaching the minimum
reduced $\chi^{2}$ (i.e., $\chi_{\nu}^{2}$) for the smallest number of
parameters. In all models considered, the centers of the Sérsic profiles were
left as free parameters, but the best-fit Sérsic profile centers in all the 2
PSF + 2 Sérsic models are consistent with the best-fit PSF centers (Figure 1
and Extended Data Figure 5), suggesting that the best-fit models are
reasonable. The only exception was the 2 PSF + 2 Sérsic model for the Keck
image, where the center of the Sérsic profile for the NE host was fixed at the
NE nucleus PSF model center to help break model degeneracy given the relative
shallowness of the AO imaging.

Figure 1 shows the best-fit GALFIT model for the HST F160W image. The extended
host galaxies are detected after the PSFs are subtracted (i.e., panel “Data-
PSFs”). After we further subtract the extended hosts (i.e., panel “Residual”),
low surface-brightness tidal features are revealed in the residual image,
providing unambiguous evidence for an ongoing galaxy merger.

Extended Data Table 5 lists the best-fit parameters of the GALFIT models for
the HST F160W image. In addition to the baseline model (with 2 PSFs + 2
Sérsic), we also consider a model for the lensing scenario (Model 2; with 2
PSFs for the quasars plus 1 Sérsic for the foreground lens galaxy) and a more
complicated model to attempt the bulge-disk decomposition (Model 3; with 2
PSFs for the quasars plus 4 Sérsic for the two hosts, where the Sérsic indices
are fixed at $n=1$ for the disk and $n=4$ for the bulge for each host to help
break model degeneracy). The baseline model is statistically preferred over
the lensing model based on the $\chi^{2}_{\nu}$ value comparison.

To more quantitatively test the lensing scenario, we have tried mass modeling
using the glafic software[38], assuming positional errors of
0$.\\!\\!^{\prime\prime}$01 for PSFs and 0$.\\!\\!^{\prime\prime}$05 for the
lensing galaxy. Assuming a simple singular isothermal ellipsoid (SIE) lens
model, we were able to obtain a mass model that fit the positions and flux
ratio of the double optical nuclei, so the lensing hypothesis cannot be
rejected based on the Model 2 best-fit positions for the HST image (Extended
Data Table 5) alone.

While the more complicated model (Model 3) yields a slightly smaller
$\chi^{2}_{\nu}$ than the baseline model, the improvement is too small to
warrant the larger number of parameters. For the HST image, the best-fit
Sérsic index in the baseline model for the host of the SW nucleus is $n=5.2$,
suggesting a bulge-dominated host (but see caveats below) with an effective
radius of 1.4 kpc, whereas that of the host of the NE nucleus is $n=0.9$,
suggesting a disk-dominated host with an effective radius of 4.3 kpc.

Extended Data Figure 2 shows the best-fit GALFIT model for the Keck AO $K_{p}$
image. While the sensitivity of the Keck image is much lower than that of the
HST image, the Keck image offers a higher angular resolution of
$\sim$0$.\\!\\!^{\prime\prime}$05 with much better PSF sampling
(0$.\\!\\!^{\prime\prime}$01 pixel-1), which is useful for better
disentangling the quasars and the hosts, characterizing the host structures,
and revealing any foreground lens galaxy in the lensing scenario. After we
subtract the two unresolved PSFs, the extended host galaxies of the two
quasars are also detected in the Keck AO image.

Extended Data Table 6 lists the best-fit parameters of GALFIT models for the
Keck AO $K_{p}$ image. The lensing scenario (i.e., with 2 PSFs + 1 Sérsic) is
statistically ruled out based on the $\chi^{2}_{\nu}$ value comparison. A more
complicated model is not warranted to attempt the bulge-disk decomposition for
the Keck image, given its much lower sensitivity, which is not enough to
detect any additional extended disks. The tidal features seen in the HST image
(in the GALFIT model residual) are too faint to detect in the Keck AO image.

We have also tried mass modeling using glafic to more quantitatively test the
lensing scenario, assuming positional errors of 0$.\\!\\!^{\prime\prime}$01
for PSFs and 0$.\\!\\!^{\prime\prime}$05 for the lensing galaxy. Assuming a
simple SIE lens model, we were unable to obtain a mass model that could fit
the positions and flux ratio of the double optical nuclei. With a more
complicated SIE plus shear model, we got a reasonably good fit, but it would
require very large external shear ($\sim$0.6) and therefore quite unrealistic.
In conclusion, our mass modeling test suggests that the Model 2 best-fit
positions for the Keck AO image (Extended Data Table 6) are inconsistent with
the lensing hypothesis.

For the Keck AO image, the best-fit Sérsic index is $n=0.2$ for the host of
the SW nucleus (with an effective radius of 1.1 kpc) and $n=0.6$ for the host
of the NE nucleus (with an effective radius of 1.8 kpc), suggesting disk-
dominated galaxies for both hosts. The Sérsic index estimate for the NE
nucleus host is broadly consistent between the HST and Keck images, although
the effective radius estimate is larger for the HST image than that for the
Keck AO image, likely due to the differences in the image resolution (i.e.,
the more compact disk component is better resolved with Keck AO) and
sensitivity (i.e., the HST image is deep enough to reveal any additional
extended disk component). On the other hand, the SW nucleus host’s
morphological type is more controversial. The discrepancy may be caused by the
different image resolutions and PSF samplings, i.e., the smaller SW host is
more challenging to robustly resolve with HST imaging given its lower angular
resolution and poorer PSF sampling. We conclude that both hosts are likely to
be disk-dominated galaxies based on the Keck-based measurements, given the
higher angular resolution (PSF FWHM of 0$.\\!\\!^{\prime\prime}$05 for Keck
compared to 0$.\\!\\!^{\prime\prime}$15 for HST) and finer PSF sampling
(0$.\\!\\!^{\prime\prime}$01 pixel-1 for Keck compared to
0$.\\!\\!^{\prime\prime}$06 pixel-1 for HST) for more robustly disentangling
the quasars and their host galaxies and measuring the host galaxy structures,
even though the HST PSF is more stable.

Simulations have shown that galaxy structure plays a dominant role in
regulating activity in major mergers: inflows in bulgeless galaxies occur
earlier in galaxy mergers whereas stellar bulges stabilize galaxies against
bar modes and inflow until the galaxies merger[39]. The disk dominance of the
hosts in SDSS J0749$+$2255 may be critical for efficient SMBH fueling to
produce two luminous quasars in the early-stage merger without the
stabilization effect of stellar bulges which are often present in low-redshift
low-luminosity dual AGNs[25].

##  0.6 HST/STIS Optical Slit Spectroscopy

Spatially resolved optical slit spectroscopy was conducted for SDSS
J0749$+$2255 on 16 February 2021 UT (Program GO 16210; PI: X. Liu) with the
HST/Space Telescope Imaging Spectrograph (STIS) using the G750L grating. A
slit width of 0$.\\!\\!^{\prime\prime}$2 was used. The slit direction was
oriented at P.A. = 236.4∘ (E of N) to simultaneously cover the two nuclei. The
total net exposure time was 2136 s. The data were reduced using the
$stis\\_cti$ package, which removes trails and other artifacts caused by
Charge Transfer Inefficiency effects[40]. The exposures were flat-fielded and
combined to reject cosmic rays. The final 1-dimensional (1D) spectra for each
nucleus were extracted using an aperture size of 7 pixels, corresponding to
0$.\\!\\!^{\prime\prime}$35\. Extended Data Figure 3 shows the HST/STIS
2-dimensional (2D) spectrum and demonstrates the extraction of the 1D spectra
for each nucleus. The two nuclei are well separated by HST/STIS. Figure 4
shows the final STIS spectra for the double nucleus.

##  0.7 Gemini/GMOS Optical Slit Spectroscopy

Gemini spatially resolved optical slit spectroscopy was conducted on April 20
2022 UT (Program GN-2022A-Q-139; PI: X. Liu) with the Multi-Object
Spectrographs (GMOS) on Gemini-North. The R150 grating and a slit width of
0$.\\!\\!^{\prime\prime}$5 were used, yielding a spectral resolution of $R\sim
631$. The slit direction was oriented at P.A.=236.1∘ (E of N) to
simultaneously cover the two nuclei. The total net exposure time was 1876 s.
The average seeing during the observations was
$\sim$0$.\\!\\!^{\prime\prime}$49 in the optical wavelengths. The spectra
cover the wavelength range of 400-1200 nm. The data were reduced using the
Gemini PyRAF/IRAF pipeline following standard procedures. The reduction steps
included bias subtraction, flat-field correction, cosmic-ray rejection,
wavelength calibration, and flux calibration. The CuAr spectra with the same
instrumental configurations were used for the wavelength calibration. EG131
was observed as the standard star for spectrophotometric calibration.

Extended Data Figure 4 shows the Gemini/GMOS 2D spectrum and demonstrates the
source decomposition of the 1D spectra for each nucleus. Based on the 1D
spatial profile of the wavelength-averaged spectrum, the two nuclei are well
modeled by two Gaussian components whose separation is consistent with that
observed in the HST optical images. We have decomposed the two nuclei using
two Gaussians using a fixed separation of 0$.\\!\\!^{\prime\prime}$46 across
all wavelengths to help with the decomposition. The FWHM of the Gaussian
components is a free parameter during the fitting process to account for the
variable seeing as a function of wavelength. Figure 4 shows the final
decomposed GMOS spectra for the double nucleus.

##  0.8 Gemini/GNIRS NIR Slit Spectroscopy

Gemini spatially resolved NIR slit spectroscopy was conducted on 22 December
2020 UT (Program GN-2020B-FT-113; PI: X. Liu) with the Near Infra-Red
Spectrograph (GNIRS) on Gemini-North. We used the cross-dispersion (XD) 32/I
mode with the blue camera and a slit width of 0$.\\!\\!^{\prime\prime}$45,
yielding a spectral solution of $R\sim 1000$. The total net exposure time was
2400 s. The average seeing during the observations was
$\sim$0$.\\!\\!^{\prime\prime}$34 in the NIR wavelengths. The detector
controller artifacts in the raw images were removed using the cleanir.py
script. Then, the data were reduced with the Gemini IRAF pipeline following
standard procedures. The reduction steps included flat-field correction,
S-distortion correction, wavelength calibration, telluric correction, and flux
calibration. The pinhole images were used to correct the S-distortion in the
images and then the Argon spectra were used for the wavelength calibration.
HIP26000 was observed as the standard star for telluric correction and
spectrophotometric calibration.

Extended Data Figure 5 shows the Gemini/GNIRS 2D spectrum and demonstrates the
source decomposition of the 1D spectra for each nucleus. Based on the 1D
spatial profile of the wavelength-averaged spectrum, the two nuclei are
marginally resolved. The 1D profile is consistent with two Gaussian components
with separation of 0$.\\!\\!^{\prime\prime}$46, although a third faint
extended component is seen near the SW nucleus, which is likely due to the
extended host galaxy emission in the merger. We adopted the local minimum
between two Gaussian as the central aperture boundaries. The final 1D spectra
for each nucleus were extracted using a boxcar aperture with a size of 21
pixels (corresponding to 1$.\\!\\!^{\prime\prime}$05) for the SW nucleus and 7
pixels (corresponding to 0$.\\!\\!^{\prime\prime}$35) for the NE nucleus. The
aperture for the NE nucleus does not fully cover the tail of the emission
($>$91 pixels) because it reaches the edge of the slit. Figure 4 shows the
final decomposed GNIRS NIR spectra for the double nucleus.

##  0.9 Joint Spectral Modeling of the HST/STIS and Gemini/GNIRS Data

We jointly fit the optical and NIR spectra of each nucleus to estimate the
systemic redshift and measure emission-line properties for each quasar. We use
the PyQSOFit fitting code[12] following the procedures as described in detail
in[41]. Here we briefly summarize the fitting procedures and model
assumptions. Each combined spectrum has been de-reddened for Galactic
extinction using the Milky Way reddening law[42] and $E(B-V)$ derived from the
dust map[43]. The spectrum is modeled as a linear combination of a pseudo-
continuum (consisting of a power-law plus polynomial quasar continuum and Fe
II emission), broad, and narrow emission lines. A Balmer continuum component
has also been considered in the models but made negligible contribution in the
best-fit results. First, to estimate the systemic redshift, we perform a
global fitting to the joint optical-NIR spectrum. Given the moderate S/N of
the spectra, the velocities of all the narrow and broad emission lines have
been tied to be the same to help break model degeneracy. The best fit systemic
redshift is $2.1679\pm 0.0003$ for the SW nucleus and $2.1665\pm 0.0010$ for
the NE nucleus which we take as our baseline estimates. Alternatively,
allowing for a nonzero velocity offset between the narrow and broad emission
lines yields a best fit systemic redshift of $2.1686\pm 0.0007$ for the SW
nucleus and $2.1653\pm 0.0021$ for the NE nucleus based on the narrow emission
lines, which are consistent with the baseline estimates within uncertainties.

Second, to make emission-line measurements and derive physical quasar
properties, we perform local fittings for each emission-line complexes. Given
the moderate S/N of the existing spectra, each narrow and broad emission line
component is modeled by a single Gaussian, where an upper limit of FWHM = 1200
km s-1 has been imposed on the narrow lines. The velocities of all the narrow
emission lines have been fixed to be the same as the best-fit systemic
redshift of each quasar, whereas non-zero velocity offsets are allowed for
each broad emission lines (i.e., with respect to the systemic redshift) to
better account for possible velocity offsets in different broad lines due to
ionization stratification, accretion disk winds, and/or IGM absorption. The
widths of all the narrow emission lines are tied to be the same to help break
model degeneracy, but the widths of the broad emission lines are allowed to
vary among different lines. The flux ratios of the [O III] and [N II] doublets
are fixed to be their theoretical values to help break model degeneracy and
reduce ambiguities in decomposing the emission line complexes. Extended Data
Figure 6 shows the best-fit spectral models for both the global and the local
fits. Extended Data Table 2 lists the emission-line measurements from the
local best-fit models.

##  0.10 Chandra ACIS-S X-ray 0.5–8 keV Imaging Spectroscopy, X-ray Data
Reduction, and Data Analysis Using BAYMAX

We observed SDSS J0749$+$2255 with the ACIS-S on board the Chandra X-ray
Observatory on 23 Dec 2021 UT (Program 23700377, ObsID=25710 and 26245; PI: X.
Liu). SDSS J0749$+$2255 was observed on axis on the S3 chip with a total
exposure time of 30 ks. We follow a standard data reduction procedure for each
observation using _Chandra_ Interactive Analysis of Observations software
(CIAO) v4.13[44]. We correct for astrometry by cross-matching _Chandra_
-detected point-like sources with available optical catalogs. We require a
minimum of 3 matches between the X-ray observations and optical catalogs, and
we required each matched pair to be separated by less than 2′′. For each
Chandra observation, we find background flaring to be negligible, with no time
interval containing a background rate 3$\sigma$ above the mean level. We
restrict our analysis to photons with energies between 0.5–8 keV. We analyze
the photons that are contained within a 20${}^{\prime\prime}\times$20′′ box
centered on the nominal X-ray coordinates of SDSS J0749$+$2255\. Extended Data
Figure 7 shows the binned Chandra data, with the soft X-rays (0.5–2 keV) shown
in red on the left (total of 114 counts) and the hard X-rays (2–8 keV) shown
in purple on the right (total of 126 counts).

We analyze the Chandra observations for the presence of two X-ray point
sources using python tool BAYMAX (Bayesian AnalYsis of Multiple AGN in
X-rays)[45, 46, 47]. BAYMAX allows for a quantitative analysis of whether a
source in a given _Chandra_ observation is more likely composed of one or two
point sources. The main component of BAYMAX is the calculation of the Bayes
factor $\mathcal{B}$, which represents the ratio of the plausibility of
observed data, given each model. BAYMAX takes calibrated _Chandra_ events and
compares them to the expected distribution of counts for a single point source
($M_{1}$) versus a dual point source ($M_{2}$) model.

When calculating the Bayes Factor, BAYMAX uses python nested sampling package
nestle, where both probability and prior densities need to be defined. The
probability densities of $M_{1}$ and $M_{2}$ are estimated by comparing the
Chandra sky coordinates ($x$, $y$) and energies ($E$) of each detected X-ray
event to simulations based on single and multiple point source models[45]. The
prior vectors for $M_{1}$ and $M_{2}$ are $\theta_{1}=\\{\mu_{1},f_{{\rm
bkg}}\\}$ and $\theta_{2}=\\{\mu_{1},\mu_{2},f,f_{{\rm bkg}}\\}$,
respectively. Here, the sky-x, sky-y coordinates of each point source
(encompassed by $\mu$) are uniform distributions across the entire
20${}^{\prime\prime}\times$20′′ field of view. The parameter $f_{{\rm bkg}}$
represents the number of counts associated with a spatially uniform background
component versus a point-source component; the prior density is represented in
log-space with log-normal distribution with a mean of $-2$ and a large
standard deviation ($\sigma$=0.5) to allow flexible movement in prior space.
Lastly, the parameter $f$ represents the count ratio between the secondary and
primary point-source; the prior density if represented in log-space with a
log-normal distribution set between $-4$ and $0$. For sources with more than
one Chandra observations, as in the case of SDSS J0749$+$2255, the code also
fits for astrometric shifts between each observation, although the point
spread function model used in the probability densities are unique to each
observation.

These are the standard “non-informative” priors used by BAYMAX when analyzing
X-ray point sources[46, 47]. Running BAYMAX on the Chandra observations of
SDSS we find log $\mathcal{B}=3.3\pm 0.7$. The error bars are returned by
nestle and represent the statistical sampling error of the nested sampling
procedure. Previous analyses have found this error to be consistent with the
$\sim 68\%$ confidence interval when running BAYMAX 100 times on the same
source[46].

In addition to calculating a Bayes factor, BAYMAX estimates optimal values for
parameters such as the separation ($r$) and count ratio ($f$) using python
MCMC package PyMC3. We find best-fit values of
$r=1\mbox{$.\\!\\!^{\prime\prime}$}54_{-0\mbox{$.\\!\\!^{\prime\prime}$}36}^{+0\mbox{$.\\!\\!^{\prime\prime}$}39}$
and log$f=-1.43_{-0.30}^{+0.30}$, at the 95% confidence interval. The
overestimation of $r$ as compared to the resolved separation of the two
optical cores (along with the large error-bars) are likely a result of the
secondary X-ray point source having, on average, only 6 counts. The low number
of counts associated with a putative secondary X-ray point source ($\sim$6)
affects the accuracy of estimations on $r$. However, the estimated position
angle between the two X-ray point sources ($350_{-24}^{+37}$ deg.) is
consistent with the two resolved cores in the _HST_ imaging (330 deg.). Deeper
_Chandra_ observations can help better discern whether the X-ray emission is
spatially consistent with the _HST_ imaging.

Sampling from the posterior distributions, the code will probabilistically
assign counts to different model components (which, for $M_{2}$ is the primary
point source, the secondary point source, and a background component), from
which a spectra can be created and analyzed. This allows for the individual
spectral fitting of each X-ray point source in a putative AGN pair. Iterating
this process results in a set of 100 spectral realizations for each model
component. We fit and model each spectral realization, for each X-ray point
source, using XSPEC[48]. This allows us to quantify the average values of the
spectral shape and X-ray luminosity, across all spectral fits.

Extended Data Figure 8 shows the spectral fits to the primary (left) and
secondary (right) X-ray point source detected by BAYMAX. We find both sets of
spectral realizations are best-fit with a simple absorbed power law (XSPEC
model components: phabs $\times$ zphabs $\times$ zpow). Given the low number
of counts for both components (average of $\sim 191.5$ and $\sim 6$ counts for
the primary and secondary, respectively), we fix the power law index,
$\Gamma$, to a value of 1.8 for all spectral fits. Extended Data Table 3 lists
the median values and the confidence interval of the best-fit spectral
parameter distributions.

##  0.11 X-ray-to-optical Luminosity Ratio as Characterized by the Spectral
Slope $\alpha_{{\rm OX}}$

The slope of a hypothetical power law from 2500 Å to 2 keV is usually defined
by $\alpha_{{\rm OX}}=0.3838\log(l_{2500\,\AA}/l_{2\,keV})$, where
$l_{2500\,\AA}$ is the rest-frame monochromatic optical luminosity (in units
of erg s-1 Hz-1) at 2500 Å measured from the dereddened spectrum for each
nucleus and $l_{2\,keV}$ is the rest-frame 2 keV luminosity density (in units
of erg s-1 Hz-1) [49]. The X-ray-to-optical spectral indices $\alpha_{{\rm
OX}}$ are calculated as 1.46$\pm$0.01 for the NE nucleus and 2.32$\pm$0.13 for
the SW nucleus. As an intrinsic probe of the accretion processes,
$\alpha_{{\rm OX}}$ characterizes the relation between the optical/UV thermal
blackbody radiation from the accretion disk and the X-ray emission from the
hot corona. $\alpha_{{\rm OX}}$ of optically selected quasars has been
observed to strongly correlate with optical luminosity with a possible
residual evolution with redshift[49]. Compared with the typical optically
selected quasars at similar luminosities[49], $\alpha_{{\rm OX}}$ of the NE
nucleus is consistent the general population. On the other hand, $\alpha_{{\rm
OX}}$ of the SW nucleus (which also has a higher estimated Eddington ratio) is
much larger, which may be evidence for the SMBH experiencing some unusual
accretion state induced by the galaxy merger. Alternatively, the large
$\alpha_{{\rm OX}}$ of the SW nucleus could be evidence of high circumnuclear
absorption as expected in galaxy mergers, although the absorption material
would have to be at the smallest scales that only affect the X-ray emitting
region but do not obscure the optical emitting region given the Type 1 nature
of the optical nucleus.

##  0.12 VLA Observations, Data Reduction, and Measurements

Dual-band radio imaging was conducted on 11 and 12 December 2020 UT (Program
20B-242; PI: X. Liu) with the VLA in its most extended A configuration. C-band
(6 GHz) and Ku-band (15 GHz) were adopted with an on-source time of 0.5 hours
for each band. The images were calibrated and reduced through the standard VLA
Calibration Pipeline version 5.6.2 in the Common Astronomy Software
Applications package (CASA)[50]. Sources 3C147 and 3C286 were used as the flux
and bandpass calibrators. Source J0738+1742 was used as the gain calibrator as
well as the phase reference source in the Ku-band. We fit 2D Gaussian
components to each nucleus to obtain the flux densities. Two nuclei are well
resolved in the Ku-band image and are consistent with two point sources. The
two nuclei are marginally resolved in the C-band image, therefore, we fit the
C-band image using two Gaussian components with a fixed shape, same as the
synthesized beam, to obtain accurate flux estimations. Figure 3 shows the VLA
A-config two-band radio images for SDSS J0749$+$2255\. Extended Data Table 4
lists the radio measurements.

##  0.13 Estimation of Virial Black Hole Masses, Bolometric Luminosities, and
Eddington Ratios

We estimate SMBH properties for each nucleus using the combined GMOS+GNIRS
spectra. These include SMBH masses, bolometric luminosities, and Eddington
ratios. The GMOS spectrum is adopted in our baseline estimates because of the
higher S/N than that of the STIS spectrum.

We estimate SMBH masses using the single-epoch virial mass estimator[41]. The
virial SMBH mass is estimated as

$\text{log}_{10}\bigg{(}\frac{M_{\text{BH}}}{M_{\odot}}\bigg{)}=a+b\text{
log}_{10}\bigg{(}\frac{\lambda L_{\lambda}}{10^{44}\,\text{erg
s}^{-1}}\bigg{)}+c\text{ log}_{10}\bigg{(}\frac{\text{FWHM}}{\text{km
s}^{-1}}\bigg{)},$ (2)

where $L_{\lambda}$ is the monochromatic continuum luminosity at wavelength
$\lambda=$5100Å, FWHM is the full width at half maximum of the H$\alpha$ line,
and $a=1.390$, $b=0.555$, $c=1.873$ are the coefficients from the calibration
of [41]. The H$\alpha$-based virial BH masses are ${\sim}10^{9.1\pm
0.4}M_{\odot}$ for the NE nucleus and ${\sim}10^{9.2\pm 0.4}M_{\odot}$ for the
SW nucleus, whereas the errors indicate total 1$\sigma$ uncertainties
including systematic errors. We have also estimated the virial BH masses using
the broad Mg II lines as an independent double check. The Mg II-based BH
masses are estimated as ${\sim}10^{8.9\pm 0.4}M_{\odot}$ for the NE nucleus
and ${\sim}10^{9.3\pm 0.4}M_{\odot}$ for the SW nucleus, which are consistent
with the H$\alpha$-based BH masses within the systematic uncertainty of
$\sim$0.4 dex[51].

We estimate the bolometric luminosity using the rest-frame optical flux
density at 5100 Å and the relevant Bolometric Correction (BC)[5]. The
bolometric luminosity is estimated as

$L_{\rm bol}=\lambda L_{\lambda}\times BC,$ (3)

where $L_{\lambda}$ is the monochromatic continuum luminosity at rest-frame
wavelength $\lambda$=5100 Å and BC = 10.33 is the bolometric correction at
5100 Å[5]. The bolometric luminosities are estimated as
$\sim$(2.0$\pm$0.4)$\times$1046 erg/s for the NE nucleus and
$\sim$(5.0$\pm$1.0)$\times$1046 erg/s for the SW nucleus.

To infer the Eddington ratio, we calculate the Eddington luminosity from the
BH mass as,

$L_{\rm Edd}\sim 1.26\times 10^{38}\bigg{(}\frac{M_{\rm
BH}}{M_{\odot}}\bigg{)}\ {\rm ergs}^{-1}.$ (4)

The Eddington ratios $L_{bol}/L_{Edd}$ are estimated as $\sim$0.12$\pm$0.03
for the NE nucleus and $\sim$0.25$\pm$0.07 for the SW nucleus.

In addition to the inferences based on the optical measurements, we also
estimate the bolometric luminosities and Eddington ratios based on the X-ray
measurements as an independent double check. We use the X-ray luminosity
$L_{2-8\,{\rm keV}}$ and apply the X-ray bolometric correction for radio-loud
quasars (Equation 14 in [52]) . The estimated bolometric luminosities are
$\sim$10${}^{46.59^{+0.03}_{-0.02}}$ erg/s for the NE nucleus and
$\sim$10${}^{45.72^{+0.17}_{-0.30}}$ erg/s for the SW nucleus. The inferred
Eddington ratios are $\sim$0.23${}^{+0.04}_{-0.03}$ for the NE nucleus and
$\sim$0.03${}^{+0.01}_{-0.02}$ for the SW nucleus. The bolometric luminosities
and Eddington ratios are in broad agreement with the optically based estimates
considering substantial scatters in the X-ray bolometric correction[52].

##  0.14 X-ray-to-Radio Multi-wavelength Spectral Energy Distributions

Extended Data Figure 9 shows the multi-wavelength spectral energy
distributions (SEDs) for both nuclei in SDSS J0749$+$2255\. The data include
photometry in the radio, infrared, optical and X-ray wavelengths, as well as
the optical and NIR spectra. For context, also shown are various optically-
selected quasar SED templates from the literature[5, 53]. Overall, the SEDs of
both nuclei are consistent with the SED templates of typical optically
selected (single) quasars. Both nuclei are classified as radio-loud quasars.

##  0.15 Estimation of Host-galaxy Stellar Masses

We estimate the total stellar mass for the host galaxy of each nucleus using
the HST F160W and Keck $K_{p}$-band PSF-subtracted host-only magnitudes. We
first obtain the HST F160W magnitudes and Keck $K_{p}$-band magnitudes of the
host galaxies by decomposing the 2D surface brightness profiles using GALFIT.
We then run the stellar population synthesis code GALAXEV[54] to create
composite stellar population (CSP) models to fit the observed two-band host-
only photometry. While only two-band data are available for the hosts, they
allow a reasonable estimate of the mass-to-light ratio (M/L) to estimate the
stellar mass. The BaSeL spectral library, a Chabrier initial mass function
(IMF)[55], and an exponentially decaying star formation history (SFH) are
assumed to build the CSP SED. The evolution tracks for three different
metallicities ($Z$=0.02, 0.008, and 0.004)[56, 57] are adopted to model any
CSP change due to different metallicities. Dust attenuation assuming a fixed
optical depth is adopted. We create a grid of CSPs from assuming different
stellar population parameters. The fitted free parameters include the total
stellar mass, age, metallicity, the exponential decay star formation
timescale, and the dust optical depth. The CSP models are fit to the PSF-
subtracted HST F160W and Keck $K_{p}$-band magnitudes to estimate the best-fit
stellar population parameters of the host galaxies.

Extended Data Figure 10 shows the best-fit CSP models and the PSF-subtracted
host-only two-band photometry. The total stellar masses are estimated as
$\sim$10${}^{11.46^{+0.02}_{-0.26}}$ $M_{\odot}$ for the NE nucleus and
10${}^{11.50^{+0.04}_{-0.27}}$ $M_{\odot}$ for the SW nucleus, whereas the
errors include 1$\sigma$ statistical uncertainties and the partial systematics
from surface brightness profile fitting. This partial systematic uncertainty
of total stellar masses is estimated to be $\sim$0.2 dex from assuming
different models during the surface brightness decomposition (e.g., 2 Sérsic
versus 1 Sérsic for each host) and should be taken as a lower limit since it
does not include the additional systematics from adopting different
assumptions about the IMF, the SFH, or the unknown dust configuration, which
may result in an additional ${\sim}0.3$ dex uncertainty at least[58].

##  0.16 Black Hole Mass versus Stellar Mass Scaling Relations

To put SDSS J0749$+$2255 into context, we study the relation between BH mass
and host-galaxy total stellar mass by comparing SDSS J0749$+$2255 with
observations of the host galaxies of single quasars and AGNs, as well as
inactive galaxies in the literature. Extended Data Figure 11 shows the
relation between BH mass and host total stellar mass for SDSS J0749$+$2255, in
comparison with various literature samples. Included for comparison are the
local samples of both inactive and active galaxies[59, 60, 61] and samples of
AGNs and quasars at $z>0.2$ selected from various methods[62, 63, 64, 65, 66].
Overall, the two nuclei of SDSS J0749$+$2255 are broadly consistent with both
the local scaling relations for inactive galaxies[59] and HST-imaged AGNs[60]
and the high-redshift quasar samples, considering the large systematic errors
of our measurements. The hosts of SDSS J0749$+$2255 already follow the local
scaling relations without prominent stellar bulges. This suggests that at
least some SMBHs may have formed before their host stellar bulges, in contrast
to naive expectations from the canonical SMBH-bulge co-evolution scenario
where SMBHs and their host galaxy stellar bulges form coevally from merger-
induced starbursts and quasar fueling[67, 68].

Given the high redshift and high luminosities of SDSS J0749$+$2255, it is
challenging to construct a large tailored sample of control quasars for a
fairer comparison with similar selection methods, properties, and data
qualities. The high angular resolution images in the rest-frame optical are
expensive to obtain for even larger samples. Previous works have studied
optical selected quasars and estimated the total stellar mass from decomposing
HST NIR images[64, 66]. However, some of the studied quasars[66] are at
significantly lower redshifts ($0.2<z<0.6$) than that of SDSS J0749$+$2255\.
Others are focused on quasars with high SMBH masses
(10${}^{9}-$10${}^{10}M_{\odot}$) only[64]. While some work also studied
optical selected quasars[63], the stellar mass was estimated from fitting the
broad-band SED rather than spatially decomposing the PSF and the host. Some
works were based on X-ray selected AGNs in the COSMOS field[62, 65], which are
likely to trace a different SMBH population than optically selected quasars
given that X-ray selected AGNs in small deep extragalactic fields tend to have
lower bolometric luminosities than typical optically selected bright quasars
from wide-area surveys.

##  0.17 Estimation of Merging Timescales

Chandrasekhar dynamical friction timescale arguments are adopted to estimate
characteristic timescales for the galaxy-galaxy merger and the inspiraling
timescales of the two SMBHs. Assuming singular isothermal spheres for the host
galaxy density profiles and circular orbits, the dynamical friction timescale
for a satellite galaxy with velocity dispersion $\sigma_{{\rm S}}$ inspiraling
from radius $r$ in a host galaxy with velocity dispersion $\sigma_{{\rm H}}$
is estimated as[69]

$t^{{\rm gal}}_{{\rm fric}}=\frac{2.7}{\ln\Lambda}\frac{r}{30\,{\rm
kpc}}\bigg{(}\frac{\sigma_{{\rm H}}}{200\,{\rm
km\,s^{-1}}}\bigg{)}^{2}\bigg{(}\frac{100\,{\rm km\,s^{-1}}}{\sigma_{{\rm
S}}}\bigg{)}^{3}~{}{\rm Gyr},$ (5)

where $\Lambda$ is the Coulomb logarithm with typical values of
3${\lesssim}\ln\Lambda{\lesssim}$30[70]. We assume $\ln\Lambda{\sim}$2, which
is appropriate for equal mass mergers[71]. Assuming the SW galaxy as the host
(primary) and the NE one as the satellite (secondary), the two galaxies would
merge in $t_{{\rm fric}}{\sim}$7${}^{+64}_{-5}$ Myr, corrected for projection
effects assuming random orientation (i.e., $r=r_{p}/\sin{<I>}$, where
$r_{p}=3.8$ kpc is the projected separation between the two galaxies and the
median inclination angle is $<I>=60^{\circ}$ for random orientations).

After merging together with their individual hosts, the two SMBHs would
inspiral to the center of the merged galaxy under dynamical friction with the
stellar background. Assuming a singular isothermal sphere for the density
distribution of the merged host galaxy, the dynamical friction timescale of a
BH of mass $M_{\rm BH}$ on a circular orbit of radius $r$ is estimated as

$t^{{\rm BH}}_{{\rm fric}}=\frac{19}{\ln\Lambda}\bigg{(}\frac{r}{5\,{\rm
kpc}}\bigg{)}^{2}\frac{\sigma_{{\rm H}}}{200\,{\rm
km\,s^{-1}}}\frac{10^{8}\,M_{\odot}}{M_{\rm BH}}~{}{\rm Gyr},$ (6)

where $\ln\Lambda{\sim}$6 for typical values[70]. We estimate $t_{{\rm
fric}}{\sim}$0.22${}^{+0.50}_{-0.16}$ Gyr for the SMBH in the NE nucleus
reaching the center of the merger products and forming a gravitationally bound
binary SMBH system with the SW nucleus, assuming a radius of $r=4$ kpc at the
start of the orbital decay.

The effects of friction by gas have been neglected in the estimation, which
will accelerate the merger processes. On the other hand, the effects of tidal
stripping of stars in the satellites have been neglected, which will delay the
merger. The effect of BH accretion have also been neglected in the simple
estimation, which may change the masses of the BHs and the resulting merger
timescales. The actual merger involving the stellar and BH components with gas
friction and gas accretion onto two SMBHs is likely much more complicated.
Numerical simulations with model parameters tailored to those of SDSS
J0749$+$2255 are needed to make more realistic predictions of the subsequent
merger of the system[72]. Regardless of the actual sequence of the mergers,
the SMBHs will form a gravitationally interacting bound binary system.

##  0.18 Implications on the Abundance of $z\sim 2$ Kpc-Scale Dual Quasars

Observational constraints on the occurrence rate of kpc-scale dual quasars at
cosmic noon are still accumulating[73, 74]. A recent study using Gaia DR3
resolved pairs around luminous quasars at $z\sim 2$ reported a dual quasar
fraction of $\sim 6\times 10^{-4}$ over $\sim 3-30\,{\rm kpc}$ pair
separations[74]. About one third of the sample have pair separations less than
$\sim 10\,{\rm kpc}$. In this measurement both members of the pair are
optically-unobscured quasars with similar physical properties as SDSS
J0749+2255. Some of the double quasars in the sample[74] are gravitationally
lensed quasars, which require further confirmation in follow-up observations;
therefore the quoted dual quasar fraction is an upper limit. Confirming SDSS
J0749$+$2255 as a genuine kpc-scale dual quasar is important in the sense that
it implies a substantial fraction of the sub-arcsec double quasars at $z\sim
2$ are physical pairs.

##  0.19 Comparison with Cosmological Hydrodynamic Simulations

To compare SDSS J0749$+$2255 with theoretical predictions for similar systems,
we utilize the very recent large-volume, high-resolution (gravitational
softening of $1.5\,{\rm kpc}/h$; dark matter mass resolution of $9.6\times
10^{6}\,M_{\odot}$) cosmological hydrodynamic simulation ASTRID[75, 76, 77].
ASTRID includes subgrid models for star formation, SMBH accretion, and the
associated supernova and AGN feedback, recently updated with an SMBH seed
population between $3\times 10^{4}\,M_{\odot}/h$ and $3\times
10^{5}\,M_{\odot}/h$ and a sub-grid dynamical friction model to follow the
SMBH dynamics down to kpc scales. With a comoving volume of $(250\,{\rm
Mpc}/h)^{3}$, ASTRID is the largest galaxy formation simulation up to date
that covers the epoch of the cosmic noon, and can provide a statistical sample
of the rare quasar pair population down to galactic scales.

By searching through a redshift range of $z=2.0$–$2.3$ in ASTRID, we have
identified 8 quasar pairs with $L_{\rm bol}>10^{46}~{}{\rm erg/s}$ for both
quasars in the pair and with 3D spatial separations of $\Delta r<30\,{\rm
kpc}$. Out of these 8 close quasar pairs, 6 only become both luminous ($L_{\rm
bol}>10^{46}~{}{\rm erg/s}$) after the host galaxy merger and when the pair
separation is below $\sim 2\,{\rm kpc}$. Only two pairs are in the early
stages of the galaxy merger and hence have separate host galaxies, similar to
SDSS J0749$+$2255\. For both pairs (PairA and PairB) we find that the quasars
each contain SMBHs with masses of $M_{{\rm BH1,2}}\sim[0.8-4]\times
10^{9}\,M_{\odot}$, and are embedded in galaxies with $M_{\rm gal}\sim$ a few
$10^{11}\,M_{\odot}$.

Extended Data Figure 12 shows the $M_{\rm BH}-M_{\rm gal}$ relation for the
eight bright quasar pairs. The colored lines mark PairA and PairB with
apparently distinct galaxy bulges, although the kinematic bulge-disk
classification adopted by simulations cannot be reliably used in interacting
systems, whose morphology is typically disturbed and the disks are in the
process of being disrupted in the merger process. PairB is the closest match
with SDSS J0749$+$2255 in both the SMBH masses and the galaxy stellar masses.

In order for the pairs to be the closest match with SDSS J0749$+$2255 (i.e.
$L_{\rm bol1,2}>10^{46}~{}{\rm erg/s}$ and in separate galaxies), they are
likely in an early phase of the galaxy merger, and in particular at their
first few pericentric passages. This is a stage where both the quasar
luminosities and the host galaxy star formation are enhanced from the merger-
induced gas inflow. Extended Data Figure 13 shows the orbits of the two
quasars in one of such pairs along with the host galaxy merger. This pair has
a very eccentric orbit with an eccentricity of $0.94$ before the merger, and
significant host galaxy tidal disruption is seen, after which the secondary
quasar becomes offset from the galaxy merger remnant. The shown pair (PairB)
merges within $\sim 80\,{\rm Myrs}$ after its properties are closest to those
of SDSS J0749$+$2255\. PairA is at a later stage of the quasar pairing than
pairB, and merges in the simulation within $\sim 10\,{\rm Myrs}$. Both merging
timescales are broadly consistent with the estimate of the SDSS J0749$+$2255
dynamical friction time (Equations 5 and 6).

##  0.20 Implications for Gravitational Waves

After dynamical friction, SDSS J0749$+$2255 is expected to undergo hardening
via loss-cone scattering of nearby stars, then GW-dominated evolution in the
pulsar-timing array (PTA) band. Given that the virial SMBH mass estimates for
the two quasars are comparable, dynamical friction should proceed without the
binary stalling at $\sim$10–100 pc[78]. Regardless of the binary mass ratio,
tidal interactions between the nuclear star clusters may accelerate the
binary’s orbital decay[79]. Depending on how gas-rich the resulting nucleus
is, the sub-pc SMBH binary’s orbital evolution may be affected by
circumbinary-disk accretion at sub-pc separations. Depending on how long this
phase lasts and how massive the gas supply is, the binary may also accrete an
appreciable amount of mass, which in turn decreases the GW coalescence
timescale[80].

Various PTAs have observed a common-spectrum stochastic process in their data
that may be consistent with a stochastic GW background[81, 82, 83]. Hellings-
Downs correlations[84], however, have not yet been observed, which are needed
in order to determine that this signal is of GW origin. If the common-process
signal that current PTAs are observing is assumed to be due to a stochastic GW
background and if the cosmic SMBH binary population at nanohertz GW
frequencies is assumed to be the sole source of this background, then the
amplitude of this signal may be used to constrain the cosmic merger rate of
SMBH binaries. This can then be compared to the observed occurrence rate of
dual AGN from current and future wide-area deep surveys. The growing number of
observations of dual AGNs/quasars across cosmic time as well as PTAs improving
their characterization of a potential stochastic GW background will continue
to provide vital tests for astrophysical models of SMBH binary formation.

##  0.21 Code availability

The code used to model the HST and Keck data is publicly available at
https://users.obs.carnegiescience.edu/peng/work/galfit/galfit.html. The code
used to perform the optical-NIR spectroscopic analysis is publicly available
at https://github.com/legolason/PyQSOFit. The BAYMAX code used to model the
Chandra data is available upon reasonable request.

##  0.22 Data availability

The SDSS spectrum (Plate = 1203, FiberID = 576, MJD = 52669) is publicly
available at https://www.sdss.org/. The HST, Chandra, VLA, Gemini, and Keck
data are all available through their separate public data archives (HST
program numbers GO-16210, and GO-16892, Chandra GO-23700377, VLA 20B-242,
Gemini 2020B-FT-113 and 2022A-Q-139).

## References

  * [1] Kormendy, J. & Richstone, D.  Inward Bound—The Search For Supermassive Black Holes In Galactic Nuclei.  _ARA &A_ 33, 581–+ (1995).
  * [2] Toomre, A. & Toomre, J.  Galactic Bridges and Tails.  _ApJ_ 178, 623–666 (1972).
  * [3] Hopkins, P. F., Hernquist, L., Cox, T. J. & Kereš, D.  A Cosmological Framework for the Co-Evolution of Quasars, Supermassive Black Holes, and Elliptical Galaxies. I. Galaxy Mergers and Quasar Activity.  _ApJS_ 175, 356–389 (2008).  0706.1243.
  * [4] Bogdanovic, T., Miller, M. C. & Blecha, L.  Electromagnetic Counterparts to Massive Black Hole Mergers.  _Living Reviews in Relativity_ (2021).  2109.03262.
  * [5] Richards, G. T. _et al._ The Sloan Digital Sky Survey Quasar Survey: Quasar Luminosity Function from Data Release 3.  _AJ_ 131, 2766–2787 (2006).  arXiv:astro-ph/0601434.
  * [6] Madau, P. & Dickinson, M.  Cosmic Star-Formation History.  _Annu. Rev. Astro. Astrophys._ 52, 415–486 (2014).
  * [7] Inada, N. _et al._ The Sloan Digital Sky Survey Quasar Lens Search. II. Statistical Lens Sample from the Third Data Release.  _AJ_ 135, 496–511 (2008).  0708.0828.
  * [8] Schechter, P. L. _et al._ First Lensed Quasar Systems from the VST-ATLAS Survey: One Quad, Two Doubles, and Two Pairs of Lensless Twins.  _AJ_ 153, 219 (2017).  1607.07476.
  * [9] Anguita, T. _et al._ The STRong lensing Insights into the Dark Energy Survey (STRIDES) 2016 follow-up campaign - II. New quasar lenses from double component fitting.  _MNRAS_ 480, 5017–5028 (2018).  1805.12151.
  * [10] Lemon, C. A., Auger, M. W., McMahon, R. G. & Ostrovski, F.  Gravitationally lensed quasars in Gaia - II. Discovery of 24 lensed quasars.  _MNRAS_ 479, 5060–5074 (2018).  1803.07601.
  * [11] Lemon, C. _et al._ The STRong lensing Insights into the Dark Energy Survey (STRIDES) 2017/2018 follow-up campaign: discovery of 10 lensed quasars and 10 quasar pairs.  _MNRAS_ 494, 3491–3511 (2020).  1912.09133.
  * [12] Shen, Y. _et al._ A hidden population of high-redshift double quasars unveiled by astrometry.  _Nature Astronomy_ 5, 569–574 (2021).  2105.03298.
  * [13] Tang, S. _et al._ Optical Spectroscopy of Dual Quasar Candidates from the Subaru HSC-SSP program.  _ApJ_ 922, 83 (2021).  2105.10163.
  * [14] Yue, M., Fan, X., Yang, J. & Wang, F.  A Candidate Kiloparsec-scale Quasar Pair at z = 5.66.  _ApJ_ 921, L27 (2021).  2110.12315.
  * [15] Lemon, C. _et al._ J1721+8842: a gravitationally lensed binary quasar with a proximate damped Lyman-$\alpha$ absorber.  _A &A_ 657, A113 (2022).  2109.01144.
  * [16] Arzoumanian, Z. _et al._ The NANOGrav 11 yr Data Set: Limits on Supermassive Black Hole Binaries in Galaxies within 500 Mpc.  _ApJ_ 914, 121 (2021).  2101.02716.
  * [17] Schneider, D. P. _et al._ The Sloan Digital Sky Survey Quasar Catalog. V. Seventh Data Release.  _AJ_ 139, 2360–2373 (2010).  1004.1167.
  * [18] Hwang, H.-C., Shen, Y., Zakamska, N. & Liu, X.  Varstrometry for Off-nucleus and Dual Subkiloparsec AGN (VODKA): Methodology and Initial Results with Gaia DR2.  _ApJ_ 888, 73 (2020).  1908.02292.
  * [19] Zezas, A., Alonso-Herrero, A. & Ward, M. J.  Searching for X-Ray Luminous Starburst Galaxies.  _Ap &SS_ 276, 601–607 (2001).  arXiv:astro-ph/9907147.
  * [20] Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B. & Green, R.  VLA Observations of Objects in the Palomar Bright Quasar Survey.  _AJ_ 98, 1195 (1989).
  * [21] Sluse, D., Hutsemékers, D., Courbin, F., Meylan, G. & Wambsganss, J.  Microlensing of the broad line region in 17 lensed quasars.  _A &A_ 544, A62 (2012).  1206.0731.
  * [22] Contreras-Santos, A. _et al._ Galaxy pairs in The Three Hundred simulations: a study on the performance of observational pair-finding techniques.  _MNRAS_ (2022).  2207.13451.
  * [23] Shapley, A. E.  Physical Properties of Galaxies from z = 2 to 4.  _Annual Review of Astronomy and Astrophysics_ 49, 525–580 (2011).
  * [24] Ding, X., Silverman, J. D. & Onoue, M.  Opening the era of quasar host studies at high redshift with JWST.  _arXiv e-prints_ arXiv:2209.03359 (2022).  2209.03359.
  * [25] Shangguan, J. _et al._ Chandra X-Ray and Hubble Space Telescope Imaging of Optically Selected Kiloparsec-scale Binary Active Galactic Nuclei. II. Host Galaxy Morphology and AGN Activity.  _ApJ_ 823, 50 (2016).  1603.07925.
  * [26] Shields, G. A. _et al._ LBQS 0103-2753: A Binary Quasar in a Major Merger.  _ApJ_ 744, 151 (2012).  1109.1524.
  * [27] Kochanek, C. S., Falco, E. E. & Muñoz, J. A.  Why Quasar Pairs are Binary Quasars and Not Gravitational Lenses.  _ApJ_ 510, 590–596 (1999).
  * [28] Hennawi, J. F. _et al._ Binary Quasars at High Redshift. I. 24 New Quasar Pairs at z ~ 3-4.  _ApJ_ 719, 1672–1692 (2010).  0908.3907.
  * [29] Hennawi, J. F., Prochaska, J. X., Cantalupo, S. & Arrigoni-Battaia, F.  Quasar quartet embedded in giant nebula reveals rare massive structure in distant universe.  _Science_ 348, 779–783 (2015).  1505.03786.
  * [30] Guy, L. P. _et al._ Rubin-Euclid Derived Data Products: Initial Recommendations.  _arXiv e-prints_ arXiv:2201.03862 (2022).  2201.03862.
  * [31] Shen, Y. _et al._ A Catalog of Quasar Properties from Sloan Digital Sky Survey Data Release 7.  _ApJS_ 194, 45–+ (2011).  1006.5178.
  * [32] Pâris, I. _et al._ The Sloan Digital Sky Survey Quasar Catalog: Fourteenth data release.  _A &A_ 613, A51 (2018).  1712.05029.
  * [33] Chen, Y.-C. _et al._ Varstrometry for Off-nucleus and Dual Subkiloparsec AGN (VODKA): Hubble Space Telescope Discovers Double Quasars.  _ApJ_ 925, 162 (2022).  2108.01672.
  * [34] Mannucci, F. _et al._ Unveiling the population of dual- and lensed- AGNs at sub-arcsec separations with Gaia.  _arXiv e-prints_ arXiv:2203.11234 (2022).  2203.11234.
  * [35] Dressel, L., 2022. “Wide Field Camera 3 Instrument Handbook, Version 14.0” (Baltimore: STScI).
  * [36] Sahu, K.C., et al., 2021, “WFC3 Data Handbook”, Version 5.0, (Baltimore: STScI).
  * [37] Peng, C. Y., Ho, L. C., Impey, C. D. & Rix, H.-W.  Detailed Decomposition of Galaxy Images. II. Beyond Axisymmetric Models.  _AJ_ 139, 2097–2129 (2010).  0912.0731.
  * [38] Oguri, M.  The Mass Distribution of SDSS J1004+4112 Revisited.  _PASJ_ 62, 1017 (2010).  1005.3103.
  * [39] Mihos, J. C. & Hernquist, L.  Gasdynamics and Starbursts in Major Mergers.  _ApJ_ 464, 641–+ (1996).  arXiv:astro-ph/9512099.
  * [40] Anderson, J. & Bedin, L. R.  An Empirical Pixel-Based Correction for Imperfect CTE. I. HST’s Advanced Camera for Surveys.  _PASP_ 122, 1035 (2010).  1007.3987.
  * [41] Shen, Y. & Liu, X.  Comparing Single-epoch Virial Black Hole Mass Estimators for Luminous Quasars.  _ApJ_ 753, 125 (2012).  1203.0601.
  * [42] Cardelli, J. A., Clayton, G. C. & Mathis, J. S.  The relationship between infrared, optical, and ultraviolet extinction.  _ApJ_ 345, 245–256 (1989).
  * [43] Schlegel, D. J., Finkbeiner, D. P. & Davis, M.  Maps of Dust Infrared Emission for Use in Estimation of Reddening and Cosmic Microwave Background Radiation Foregrounds.  _ApJ_ 500, 525–+ (1998).  arXiv:astro-ph/9710327.
  * [44] Fruscione, A. _et al._ CIAO: Chandra’s data analysis system.  In Silva, D. R. & Doxsey, R. E. (eds.) _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , vol. 6270 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 62701V (2006).
  * [45] Foord, A. _et al._ A Bayesian Analysis of SDSS J0914+0853, a Low-mass Dual AGN Candidate.  _ApJ_ 877, 17 (2019).  1904.06363.
  * [46] Foord, A. _et al._ A Second Look at 12 Candidate Dual AGNs Using BAYMAX.  _ApJ_ 892, 29 (2020).  2002.01033.
  * [47] Foord, A., Gültekin, K., Runnoe, J. C. & Koss, M. J.  AGN Triality of Triple Mergers: Detection of Faint X-Ray Point Sources.  _ApJ_ 907, 71 (2021).  2012.00761.
  * [48] Arnaud, K. A.  XSPEC: The First Ten Years.  In Jacoby, G. H. & Barnes, J. (eds.) _Astronomical Data Analysis Software and Systems V_ , vol. 101 of _Astronomical Society of the Pacific Conference Series_ , 17 (1996).
  * [49] Green, P. J. _et al._ A Full Year’s Chandra Exposure on Sloan Digital Sky Survey Quasars from the Chandra Multiwavelength Project.  _ApJ_ 690, 644–669 (2009).  0809.1058.
  * [50] McMullin, J. P., Waters, B., Schiebel, D., Young, W. & Golap, K.  CASA Architecture and Applications.  In Shaw, R. A., Hill, F. & Bell, D. J. (eds.) _Astronomical Data Analysis Software and Systems XVI_ , vol. 376 of _Astronomical Society of the Pacific Conference Series_ , 127 (2007).
  * [51] Shen, Y.  The mass of quasars.  _Bulletin of the Astronomical Society of India_ 41, 61–115 (2013).  1302.2643.
  * [52] Runnoe, J. C., Brotherton, M. S. & Shang, Z.  Updating quasar bolometric luminosity corrections.  _MNRAS_ 422, 478–493 (2012).  1201.5155.
  * [53] Shang, Z. _et al._ The Next Generation Atlas of Quasar Spectral Energy Distributions from Radio to X-Rays.  _ApJS_ 196, 2 (2011).  1107.1855.
  * [54] Bruzual, G. & Charlot, S.  Stellar population synthesis at the resolution of 2003.  _MNRAS_ 344, 1000–1028 (2003).  astro-ph/0309134.
  * [55] Chabrier, G.  Galactic Stellar and Substellar Initial Mass Function.  _PASP_ 115, 763–795 (2003).  astro-ph/0304382.
  * [56] Fagotto, F., Bressan, A., Bertelli, G. & Chiosi, C.  Evolutionary sequences of stellar models with new radiative opacities. IV. Z=0.004 and Z=0.008.  _A &AS_ 105, 29–38 (1994).
  * [57] Bressan, A., Fagotto, F., Bertelli, G. & Chiosi, C.  Evolutionary Sequences of Stellar Models with New Radiative Opacities. II. Z = 0.02.  _A &AS_ 100, 647 (1993).
  * [58] Conroy, C.  Modeling the Panchromatic Spectral Energy Distributions of Galaxies.  _Annu. Rev. Astro. Astrophys._ 51, 393–455 (2013).
  * [59] Kormendy, J. & Ho, L. C.  Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies.  _ARA &A_ 51, 511–653 (2013).  1304.7762.
  * [60] Bennert, V. N. _et al._ A Local Baseline of the Black Hole Mass Scaling Relations for Active Galaxies. IV. Correlations Between MBH and Host Galaxy $\sigma$, Stellar Mass, and Luminosity.  _ApJ_ 921, 36 (2021).  2101.10355.
  * [61] Reines, A. E. & Volonteri, M.  Relations between Central Black Hole Mass and Total Galaxy Stellar Mass in the Local Universe.  _ApJ_ 813, 82 (2015).  1508.06274.
  * [62] Jahnke, K. _et al._ Massive Galaxies in COSMOS: Evolution of Black Hole Versus Bulge Mass but not Versus Total Stellar Mass Over the Last 9 Gyr?  _ApJ_ 706, L215–L220 (2009).  0907.5199.
  * [63] Dong, X. Y. & Wu, X.-B.  Herschel Observed Stripe 82 Quasars and Their Host Galaxies: Connections between AGN Activity and host Galaxy Star Formation.  _ApJ_ 824, 70 (2016).  1604.06189.
  * [64] Mechtley, M. _et al._ Do the Most Massive Black Holes at z = 2 Grow via Major Mergers?  _ApJ_ 830, 156 (2016).  1510.08461.
  * [65] Suh, H. _et al._ No Significant Evolution of Relations between Black Hole Mass and Galaxy Total Stellar Mass Up to z $\sim$ 2.5.  _ApJ_ 889, 32 (2020).  1912.02824.
  * [66] Li, J. I. H. _et al._ The Sloan Digital Sky Survey Reverberation Mapping Project: The MBH-Host Relations at 0.2 $\lesssim$ z $\lesssim$ 0.6 from Reverberation Mapping and Hubble Space Telescope Imaging.  _ApJ_ 906, 103 (2021).  2006.02522.
  * [67] Di Matteo, T., Springel, V. & Hernquist, L.  Energy input from quasars regulates the growth and activity of black holes and their host galaxies.  _Nature_ 433, 604–607 (2005).  arXiv:astro-ph/0502199.
  * [68] Di Matteo, T., Colberg, J., Springel, V., Hernquist, L. & Sijacki, D.  Direct Cosmological Simulations of the Growth of Black Holes and Galaxies.  _ApJ_ 676, 33–53 (2008).  0705.2269.
  * [69] Chandrasekhar, S.  Dynamical Friction. I. General Considerations: the Coefficient of Dynamical Friction.  _ApJ_ 97, 255 (1943).
  * [70] Binney, J. & Tremaine, S.  _Galactic dynamics_ (1987).
  * [71] Dubinski, J., Mihos, J. C. & Hernquist, L.  Constraining Dark Halo Potentials with Tidal Tails.  _ApJ_ 526, 607–622 (1999).  arXiv:astro-ph/9902217.
  * [72] Renaud, F., Appleton, P. N. & Xu, C. K.  N-body Simulation of the Stephan’s Quintet.  _ApJ_ 724, 80–91 (2010).  1009.2740.
  * [73] Silverman, J. D. _et al._ Dual Supermassive Black Holes at Close Separation Revealed by the Hyper Suprime-Cam Subaru Strategic Program.  _ApJ_ 899, 154 (2020).  2007.05581.
  * [74] Shen, Y. _et al._ Statistics of Galactic-Scale Quasar Pairs at Cosmic Noon.  _arXiv e-prints_ arXiv:2208.04979 (2022).  2208.04979.
  * [75] Bird, S. _et al._ The ASTRID simulation: galaxy formation and reionization.  _MNRAS_ 512, 3703–3716 (2022).  2111.01160.
  * [76] Ni, Y. _et al._ The ASTRID simulation: the evolution of supermassive black holes.  _MNRAS_ 513, 670–692 (2022).  2110.14154.
  * [77] Chen, N. _et al._ Massive black hole mergers with orbital information: predictions from the ASTRID simulation.  _MNRAS_ (2022).  2112.08555.
  * [78] Dosopoulou, F. & Antonini, F.  Dynamical Friction and the Evolution of Supermassive Black Hole Binaries: The Final Hundred-parsec Problem.  _ApJ_ 840, 31 (2017).
  * [79] Ogiya, G., Hahn, O., Mingarelli, C. M. F. & Volonteri, M.  Accelerated orbital decay of supermassive black hole binaries in merging nuclear star clusters.  _Monthly Notices of the Royal Astronomical Society_ 493, 3676–3689 (2020).
  * [80] Siwek, M. S., Kelley, L. Z. & Hernquist, L.  The effect of differential accretion on the gravitational wave background and the present-day MBH binary population.  _Monthly Notices of the Royal Astronomical Society_ 498, 537–547 (2020).
  * [81] Arzoumanian, Z. _et al._ The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background.  _ApJL_ 905, L34 (2020).
  * [82] Chen, S. _et al._ Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search.  _Monthly Notices of the Royal Astronomical Society_ 508, 4970–4993 (2021).
  * [83] Goncharov, B. _et al._ On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array.  _ApJL_ 917, L19 (2021).  Publisher: American Astronomical Society.
  * [84] Hellings, R. W. & Downs, G. S.  Upper limits on the isotropic gravitational radiation background from pulsar timing analysis.  _ApJ_ 265, L39–L42 (1983).

We thank M. Leveille, A. Vick, R. Campbell, R. McGurk, J. Cortes, T. R.
Geballe, S. Leggett, A. Nitta, T. Seccull, and H. Medlin for their help with
our HST, Keck, Gemini, and VLA observations. Y.-C.C. thanks J. Li and H. Guo
for their help with the GALFIT and PyQSOFit codes. M.O. acknowledges support
from JSPS KAKENHI grants JP20H00181, JP20H05856, JP22H01260. This work is
supported by the Heising-Simons Foundation and Research Corporation for
Science Advancement and NSF grants AST-2108162 and AST-2206499. This research
was supported in part by the National Science Foundation under PHY-1748958.
Support for Program number 23700237 was provided by NASA through Chandra Award
Number GO2-23099X issued by the Chandra X-ray Observatory Center, which is
operated by the Smithsonian Astrophysical Observatory for and on behalf of
NASA under contract NAS 8-03060. Support for Program number HST-GO-15900 (PI:
H. Hwang), HST-GO-16210, and HST-GO-16892 was provided by NASA through grants
from the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Incorporated, under
NASA contract NAS5-26555. The National Radio Astronomy Observatory is a
facility of the National Science Foundation operated under cooperative
agreement by Associated Universities, Inc. Based in part on data obtained at
the W. M. Keck Observatory, which is operated as a scientific partnership
among the California Institute of Technology, the University of California and
the National Aeronautics and Space Administration. The Observatory was made
possible by the generous financial support of the W. M. Keck Foundation. This
research has made use of the Keck Observatory Archive (KOA), which is operated
by the W. M. Keck Observatory and the NASA Exoplanet Science Institute
(NExScI), under contract with the National Aeronautics and Space
Administration. The authors wish to recognize and acknowledge the very
significant cultural role and reverence that the summit of Maunakea has always
had within the indigenous Hawaiian community. We are most fortunate to have
the opportunity to conduct observations from this mountain. Based in part on
observations obtained at the international Gemini Observatory (Program IDs
2020B-FT-113 and GN-2022A-Q-139; PI: X. Liu), a program of NSF’s NOIRLab,
which is managed by the Association of Universities for Research in Astronomy
(AURA) under a cooperative agreement with the National Science Foundation. on
behalf of the Gemini Observatory partnership: the National Science Foundation
(United States), National Research Council (Canada), Agencia Nacional de
Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e
Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e
Comunicações (Brazil), and Korea Astronomy and Space Science Institute
(Republic of Korea). This work was enabled by observations made from the
Gemini North telescope, located within the Maunakea Science Reserve and
adjacent to the summit of Maunakea. We are grateful for the privilege of
observing the Universe from a place that is unique in both its astronomical
quality and its cultural significance. This work makes use of SDSS-I/II and
SDSS-III/IV data (http://www.sdss.org/ and http://www.sdss3.org/).

Y.-C.C. conducted the VLA and Keck observations, reduced and analyzed the HST,
Gemini, and VLA data, performed the GALFIT analysis and the optical and NIR
spectroscopic modeling, and produced Figures 1, 3, and 4, Extended Data
Figures 1–6 and 9–11, and Extended Data Tables 1, 2, and 4–6. X.L. led the
study, was PI of the Chandra, Gemini, HST, and VLA observing programs,
conducted the Keck observations, and reduced the Keck data. A.F. analyzed the
Chandra data, performed the BAYMAX analysis and produced Figure 2, Extended
Data Figures 7 and 8, and Extended Data Table 3. Y.S. was PI of the Keck
observing program, conducted the Keck observations, and wrote the optical and
NIR spectroscopic modeling pipeline. N.C. produced Extended Data Figures 12
and 13. M.O. performed the lensing mass model tests. X.L., Y.-C.C., and A.F.
co-wrote the manuscript with the help of N.C., Y.S., and M.H.. All authors
contributed to the results and commented on the manuscript.

Reprints and permissions information is available at www.nature.com/reprints.
The authors declare no competing financial interests. Publisher’s note:
Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations. Correspondence and requests for
materials should be addressed to X.L. (xinliuxl@illinois.edu).

|  | F475W | F814W | F160WPSF | $K_{p,{\rm PSF}}$ | F160Whost | $K_{p,{\rm host}}$ | log$M_{\ast}$ | log$M_{\rm BH}$ | log$M_{{\rm BH,\,vir}}$
---|---|---|---|---|---|---|---|---|---|---
Nucleus | Redshift | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | ($M_{\odot}$) | ($M_{\odot}$) | ($M_{\odot}$)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
SDSS J0749$+$2255 NE | 2.1665 $\pm$ 0.0010 | 21.08 | 20.56 | 20.61 | 20.85 | 21.35 | 20.75 | 11.46 | 9.07 | 9.12
SDSS J0749$+$2255 SW | 2.1679 $\pm$ 0.0003 | 19.48 | 19.19 | 19.34 | 19.62 | 21.05 | 20.51 | 11.50 | 9.12 | 9.20
Extended Data Table 1: Basic photometric and spectroscopic properties of SDSS J0749$+$2255\. Col. 1: NE: northeast nucleus. SW: southwest nucleus; Col. 2: Spectroscopic redshift measured from the best-fit model of the joint optical-NIR spectral fitting (Extended Data Figure 6); Cols. 3 & 4: HST optical dual-band AB magnitude from [33]; Cols. 5 & 6: HST and Keck NIR PSF AB magnitude; Cols. 7 & 8: HST and Keck NIR model AB magnitude for the extended host galaxies. Col. 9: Total host galaxy stellar mass estimate derived from $K_{p,\,{\rm host}}$ and the F160Whost \- $K_{p,\,{\rm host}}$ color (Extended Data Figure 10). Col. 10: Black hole mass estimate inferred from host galaxy total stellar mass assuming the $M_{\rm BH}$–$M_{\ast}$ relation observed in local HST-imaged broad-line AGNs with virial BH mass estimates[60]. Col. 11: Virial black hole mass estimate inferred from broad H$\alpha$ assuming the calibration of [41].  Emission-line Measurements | SDSS J0749$+$2255 NE | SDSS J0749$+$2255 SW
---|---|---
C III] Flux (10-17 erg s-1 cm-2) (1) | 246$\pm$13 | 736$\pm$16
Mg II Flux (10-17 erg s-1 cm-2) (2) | 110$\pm$4 | 466$\pm$8
H$\beta$ Flux (10-17 erg s-1 cm-2) (3) | 51$\pm$14 | 233$\pm$27
[O III] $\lambda$5007 Flux (10-17 erg s-1 cm-2) (4) | 8$\pm$4 | 22$\pm$9
Broad H$\alpha$ Flux (10-17 erg s-1 cm-2) (5) | 151$\pm$10 | 309$\pm$21
Narrow H$\alpha$ Flux (10-17 erg s-1 cm-2) (6) | 8$\pm$3 | 21$\pm$10
C III] EW (Å) (7) | 29$\pm$2 | 22$\pm$1
Mg II EW (Å) (8) | 28$\pm$1 | 32$\pm$1
H$\beta$ EW (Å) (9) | 45$\pm$14 | 81$\pm$13
[O III] $\lambda$5007 EW (Å) (10) | 7$\pm$4 | 8$\pm$3
Broad H$\alpha$ EW (Å) (11) | 181$\pm$16 | 150$\pm$14
Narrow H$\alpha$ EW (Å) (12) | 9$\pm$4 | 10$\pm$5
C III] FWHM (km s-1) (13) | 8480$\pm$1200 | 7900$\pm$290
Mg II FWHM (km s-1) (14) | 3980$\pm$160 | 4450$\pm$100
H$\beta$ FWHM (km s-1) (15) | 7600$\pm$1800 | 7400$\pm$1000
[O III] $\lambda$5007 FWHM (km s-1) (16) | 910$\pm$410 | 680$\pm$200
Broad H$\alpha$ FWHM (km s-1) (17) | 5550$\pm$400 | 4680$\pm$430
Extended Data Table 2: Emission-line measurements from jointly analyzing the Gemini/GMOS and Gemini/GNIRS spectra. Lines (1)–(6): emission-line flux intensity. For Mg II , only the broad component is listed here. Lines (7)–(12): rest-frame equivalent width. Lines (13)–(17): Full width at half maximum. All errors quoted are 1-$\sigma$ uncertainties estimated from Monte Carlo simulations.  |  | flux${}_{0.5-2\,{\rm keV}}$ | flux${}_{2-8\,{\rm keV}}$ |  |  | NH | log $L_{0.5-2\,{\rm keV}}$ | log $L_{2-8\,{\rm keV}}$
---|---|---|---|---|---|---|---|---
Nucleus | Counts | (10-6 photons cm-2 s-1) | HR | $\Gamma_{X}$ | (1022 cm-2) | (erg s-1) | (erg s-1)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
SDSS J0749$+$2255 NE | 191.5 | 12.26${}^{+3.71}_{-1.06}$ | 10.52${}^{+0.24}_{-0.48}$ | 0.37${}^{+0.01}_{-0.02}$ | 1.8 | 4.71${}^{+1.3}_{-2.6}$ | 44.92${}^{+0.02}_{-0.03}$ | 45.29${}^{+0.05}_{-0.03}$
SDSS J0749$+$2255 SW | 6 | 0.56${}^{+0.51}_{-0.41}$ | 0.23${}^{+0.21}_{-0.17}$ |  $-$0.34${}^{+0.48}_{-0.60}$ | 1.8 |  $<$10-2 | 43.17${}^{+0.32}_{-0.60}$ | 43.62${}^{+0.32}_{-0.58}$
Extended Data Table 3: X-ray measurements from Chandra ASIS-S observations. Col. 2: Total net counts in 0.5–8 keV band. Cols. 3 & 4: Observed photon flux. Col. 5: Hardness ratio HR$\equiv(H-S)/(H+S)$ where $H$ and $S$ are the number of counts in the hard and soft X-ray bands, respectively. Col. 6: Photon index value (fixed), assuming a power-law model where $n(E)\propto E^{-\Gamma_{X}}$. Col. 7: Best-fit intrinsic column density assuming a power-law model. Cols. 8 & 9: Unabsorbed luminosity. All values and errors represent the median value and 99.7% confidence interval when analyzing 100 spectral realizations of each X-ray point source.  | S${}^{{\rm Peak}}_{{\rm 6.0\,GHz}}$ | S${}^{{\rm Int}}_{{\rm 6.0\,GHz}}$ | log $\nu$L${}_{{\rm\nu,6.0\,GHz}}$ | S${}^{{\rm Peak}}_{{\rm 15.0\,GHz}}$ | S${}^{{\rm Int}}_{{\rm 15.0\,GHz}}$ | log $\nu$L${}_{{\rm\nu,15.0\,GHz}}$
---|---|---|---|---|---|---
Nucleus | (mJy/beam) | (mJy) | (erg s-1) | (mJy/beam) | (mJy) | (erg s-1 )
(1) | (2) | (3) | (4) | (5) | (6) | (7)
SDSS J0749$+$2255 NE | 22.91$\pm$0.05 | 22.94$\pm$0.05 | 43.382$\pm$0.002 | 15.64$\pm$0.08 | 15.95$\pm$0.14 | 43.622$\pm$0.004
SDSS J0749$+$2255 SW | 60.43$\pm$0.05 | 60.51$\pm$0.05 | 43.763$\pm$0.001 | 44.66$\pm$0.08 | 45.38$\pm$0.14 | 44.036$\pm$0.002
Extended Data Table 4: Radio measurements from VLA A-config dual-band imaging. Col. 2: 6.0 GHz (C-band) peak flux density. Col. 3: 6.0 GHz integrated flux density from decomposing the two marginally resolved nuclei. Col. 4: Rest-frame 6.0 GHz luminosity density. Col. 5: VLA 15.0 GHz (Ku band) peak flux density. Col. 6: VLA 15.0 GHz integrated flux density from decomposing the two well resolved nuclei. Col. 7: Rest-frame 15.0 GHz luminosity density. All errors quoted represent 1$\sigma$ statistical uncertainties.  HST (F160W) GALFIT Parameters | SDSS J0749$+$2255 NE | SDSS J0749$+$2255 SW
---|---|---
Model 1 (2 PSFs + 2 Sérsics): $\chi^{2}_{\nu}=$ 1.331 |  |
mPSF (AB) (1) | 20.605$\pm$0.002 | 19.344$\pm$0.019
m${}_{\rm S\acute{e}rsic}$ (AB) (2) | 21.349$\pm$ 0.034 | 21.053$\pm$0.095
$R_{e}$ (kpc)(3) | 4.29$\pm$0.07 | 1.44$\pm$0.28
n(4) | 0.91$\pm$0.03 | 5.20$\pm$1.32
Model 2 (2 PSFs + 1 Sérsic): $\chi^{2}_{\nu}=$1.354 |  |
mPSF (AB) (5) | 20.613$\pm$0.003 | 19.260$\pm$0.001
CPSF (hh:mm:ss.ss+dd:mm:ss.s)(6) | 07:49:23.00+22:55:12.0 | 07:49:22.97+22:55:11.8
m${}_{\rm S\acute{e}rsic}$ (AB) (7) | 20.734$\pm$0.008
C${}_{\rm S\acute{e}rsic}$ (hh:mm:ss.ss+dd:mm:ss.s)(8) | 07:49:22.98+22:55:11.9
$R_{e}$ (kpc)(9) | 3.47$\pm$0.04
n(10) | 1.38$\pm$0.04
Model 3 (2 PSFs + 4 Sérsics): $\chi^{2}_{\nu}=$1.327 |  |
mPSF (AB) (11) | 20.613$\pm$0.003 | 19.373$\pm$0.014
m${}_{\rm S\acute{e}rsic,\,1}$ (AB) (12) | 21.476$\pm$0.038 | 21.257$\pm$0.057
m${}_{\rm S\acute{e}rsic,\,2}$ (AB) (13) | 23.794$\pm$0.876 | 22.186$\pm$0.157
$R_{e,\,1}$ (kpc)(14) | 3.89$\pm$0.09 | 1.94$\pm$0.29
$R_{e,\,2}$ (kpc)(15) | 9.94$\pm$13.91 | 0.81$\pm$0.14
n1(16) | 1 | 4
B/T (17) | 0.11$\pm$0.08 | 0.70$\pm$0.05
Extended Data Table 5: GALFIT structural decomposition results for the HST/WFC3 F160W image. Lines (1) & (2): Best-fit component magnitude. Lines (3) & (4): Best-fit Sérsic component effective radius and index. Lines (5)–(10): Similar to above, but for a single Sérsic for the foreground lens. Lines (6) & (8): J2000 coordinates for the best-fit centers. Lines (11)–(17): Similar to above, but for two Sérsic models for each host to attempt bulge-disk decomposition. The Sérsic index is fixed at $n=1$ ($n=4$) for the disk (bulge) in each host to help break model degeneracy. Line (17): Bulge-to-total luminosity ratio.  Keck ($K_{p}$) GALFIT Parameters | SDSS J0749$+$2255 NE | SDSS J0749$+$2255 SW
---|---|---
Model 1 (2 PSFs + 2 Sérsics): $\chi^{2}_{\nu}=$1.416 |  |
mPSF (AB) (1) | 20.847$\pm$0.006 | 19.620$\pm$0.002
m${}_{\rm S\acute{e}rsic}$ (AB) (2) | 20.745$\pm$0.017 | 20.514$\pm$0.008
$R_{e}$ (kpc)(3) | 1.76$\pm$0.04 | 1.07$\pm$0.01
n(4) | 0.56$\pm$0.05 | 0.16$\pm$0.02
Model 2 (2 PSFs + 1 Sérsic): $\chi^{2}_{\nu}=$1.521 |  |
mPSF (AB) (5) | 20.706$\pm$0.004 | 19.573$\pm$0.002
CPSF (hh:mm:ss.ss+dd:mm:ss.s)(6) | 07:49:23.00+22:55:12.0 | 07:49:22.97+22:55:11.8
m${}_{\rm S\acute{e}rsic}$ (AB) (7) | 19.925$\pm$0.011
C${}_{\rm S\acute{e}rsic}$ (hh:mm:ss.ss+dd:mm:ss.s)(8) | 07:49:22.97+22:55:11.7
$R_{e}$ (kpc)(9) | 1.79$\pm$0.04
n(10) | 1.405$\pm$0.04
Extended Data Table 6: GALFIT structural decomposition results for the
Keck/NIRC2 AO $K_{p}$ image. Lines (1) & (2): Best-fit component magnitude.
Lines (3) & (4): Best-fit Sérsic component effective radius and index. Lines
(5)–(10): Similar to above, but for a single Sérsic model for the foreground
lens. Lines (6) & (8): J2000 coordinates for the best-fit centers. The
absolute astrometry of the Keck image has been registered to that of the HST
image according to the PSF centers.

Figure 1: HST/WFC3 F160W (H band) image and 2D structural decomposition
results using GALFIT[37] modeling for SDSS J0749$+$2255. From left to right:
the raw HST image, our best-fit model consisting of two PSFs for the central
quasars and two Sérsic models for the two host galaxies, the two-PSF-
subtracted image highlighting the detection of two host galaxies underlying
both quasars in the merger, the residual image after subtracting two PSFs and
two Sérsic models finally revealing the detection of low surface-brightness
tidal tails indicative of ongoing galactic interactions, and the 1D radial
profiles of the data and the best-fit model and the individual components. The
detection of both the host galaxies and the tidal tails unambiguously confirms
the system as a dual quasar rather than a lensed quasar or a pre-merger
projected quasar pair at larger physical separations. In addition to the two
Sérsic models for the two quasar host galaxies, we have also considered the
alternative scenario of a single Sérsic model for a central foreground lens
galaxy which is disfavored by the data based on the $\chi^{2}_{\nu}$ values.
Figure 2: Chandra ASIS-S X-ray image and best-fit BAYMAX results for SDSS
J0749$+$2255\. From left to right: Binned Chandra data (0.5$-$8 keV), centered
on the nominal X-ray coordinates of quasar with best-fit location of primary
X-ray point source (red “x”) and secondary X-ray point source (red “+”) with
95% confidence interval shown in red contour; joint-posterior distribution of
the separation (_r_ , arcsecond) and count ratio (log $f$; ratio between the
number of counts associated with secondary versus the primary). Blue contour
levels represent 68%, and 95% confidence intervals. We find the best-fit
separation
$r=1\mbox{$.\\!\\!^{\prime\prime}$}54_{-0\mbox{$.\\!\\!^{\prime\prime}$}36}^{+0\mbox{$.\\!\\!^{\prime\prime}$}39}$
and log$f=-1.43_{-0.30}^{+0.30}$, at the 95% confidence interval. Although the
estimated separation is statistically larger than estimated from _HST_ imaging
(0$.\\!\\!^{\prime\prime}$46), the estimated position angle between the two
X-ray point sources ($350_{-24}^{+37}$ deg.) is consistent with the position
angle of the two resolved cores resolved in the _HST_ imaging ($330$ deg.).
Figure 3:  VLA A-config Ku-band (Left) and C-band (Right) continuum images for
SDSS J0749$+$2255\. The contours represent 10, 30, 60, 120, 200, 300, 400,
600, 1000 $\times$ 1$\sigma$ noise level. The green crosses mark the centroids
of the fitted Gaussian components. The synthesized beam sizes are shown in
grey in the bottom left corner.  Figure 4:  Spatially resolved optical and NIR
spectroscopy from HST/STIS, Gemini/GMOS, and Gemini/GNIRS for SDSS
J0749$+$2255\. The orange and cyan curves represent Gemini spectra for the
northeastern (NE) and southwestern (SW) nuclei, respectively. The green and
magenta curves show the HST spectra for the NE and SW nuclei. Labeled are the
detected major strong emission lines. Bottom panel shows the flux ratio
between the two nuclei, where the HST-based flux ratio is shown in red and the
Gemini-based flux ratio (for both GMOS and GNIRS) is shown in black. The
circled-cross symbols exhibit the telluric absorption region. The two nuclei
have significantly different rest-frame UV continuum spectral slopes and
strengths of the broad and narrow emission lines. The discrepancy between the
STIS- and GMOS-based flux ratios is likely cause by differences in the slit
coverage and/or observing time. We adopt the STIS spectra to estimate the
continuum power-law spectral indices because STIS cleanly separated the two
nuclei.  Extended Data Figure 1: SDSS J0749$+$2255 in comparison to the
previously known candidate dual quasars in the literature at $z>1$ and with
projected separations $<$10 kpc. Nearly identical quasars (NIQs)[9, 10, 11],
shown with grey symbols, are marked with ”(NIQ)” in the legend. They are
generally thought to be lensed quasars lacking a definitive detection of
foreground lenses due to observational limitations. Three other dual quasar
candidates have more noticeable spectral differences[7, 8, 13] (shown in
magenta), but they do not show unambiguous host merger features, and are
debatable as project pairs at much larger separations or lenses. The dual
quasar candidate of [12] (shown in blue) shows differences in the broad
emission line strengths[34] but is still controversial as a quasar lens due to
the lack of high-resolution deep NIR imaging. The quadruply lensed quasar[15]
(shown in green) has been suggested to contain a lensed kpc-scale dual quasar,
but the interpretation is subject to lense modeling.  Extended Data Figure 2:
Similar to Figure 1, but for Keck AO-assisted $K_{p}$ imaging of SDSS
J0749$+$2255 and 2D structural decomposition results from GALFIT analysis. The
best-fit GALFIT model contains two PSFs for the two unresolved quasars and two
Sérsic models for the two host galaxies with the $\chi^{2}_{\nu}$ value
labeled on the residual image. After we subtract the two PSFs, the extended
host galaxies are detected (i.e., Data-PSFs). The tidal features detected in
the HST image (in the GALFIT model residuals) are too faint to detect in the
shallower Keck image. The lensing hypotheis is disfavored based on both the
$\chi^{2}_{\nu}$ values and lensing mass model tests.

Extended Data Figure 3: 2D HST/STIS spectrum for SDSS J0749$+$2255 and
illustration of the slit position and 1D aperture extraction for each nucleus.
Top Left: Illustration of the STIS slit position and orientation. The
background shows the HST optical two-band (F475W and F814W) color image[12].
The width of the white rectangle denotes the slit width. Top Right: Projected
1D spatial profile of the wavelength-averaged STIS spectrum. The blue filled
circles denote the total spectrum. The black vertical lines mark the
extraction apertures for the two nuclei. Bottom: 2D STIS spectrum around the C
III] line. The white horizontal lines mark the extraction apertures for the
two nuclei.

Extended Data Figure 4: 2D Gemini/GMOS spectrum for SDSS J0749$+$2255 and
illustration of the slit position and source decomposition for each nucleus.
Top Left: Illustration of the GMOS slit position and orientation. The
background shows the HST optical two-band (F475W and F814W) color image[12].
The width of the white rectangle denotes the slit width. Top Right: Projected
1D spatial profile of the wavelength-averaged GMOS spectrum. The blue filled
circles denote the total spectrum. The green (purple) open circles represent
the decomposed spatial profile of the SW (NE) nucleus. Bottom: 2D GMOS
spectrum around the C III] line.

Extended Data Figure 5: 2D Gemini/GNIRS spectrum for SDSS J0749$+$2255 and
illustration of the slit position and Top: The GNIRS slit position and
orientation. The background shows the HST F160W PSF-subtracted (top left) and
residual (top right) images. The white rectangles marked with “NE” and “SW”
show the extraction apertures of both nuclei. The three dashed lines mark the
peak positions of three components seen in the projected 1D profile. Bottom
Left: Projected 1D spatial profile of the wavelength-average GNIRS spectrum.
The blue filled circles show the total spectrum. The open circles represent
the SW nucleus (shown in green), the NE nucleus (purple), and the extended
host (orange). The top horizontal black lines mark the extraction apertures
for both nuclei. Bottom Right: 2D GNIRS spectrum around the H$\alpha$ line.
The horizontal black lines denote the extraction apertures for both nuclei.

Extended Data Figure 6: Global and local spectral fitting analysis for the
joint optical-NIR spectrum of the two nuclei in SDSS J0749$+$2255\. The top
and bottom two rows show results for the SW and NE nuclei, respectively. For
each nucleus, the upper panels show the global fit of the pseudo-continuum (in
cyan), which consists of a power-law plus polynomial (in orange), and a fit
with the Fe II template. The horizontal grey bars denote the wavelength
windows adopted in the pseudo-continuum fit. The lower panels show the local
fit for individual emission lines (with the $\chi^{2}_{\nu}$ values marked)
including C IV through H$\alpha$ using the pseudo-continuum-subtracted
spectra. Red curves show models for the broad emission lines, whereas the
green and blue curves denote models for the narrow and total emission lines,
respectively.  Extended Data Figure 7: Binned Chandra data, centered on the
nominal X-ray coordinates of SDSS J0749$+$2255\. Soft photons (0.5–2 keV) are
shown in red on the left (total of 114 counts) and hard photons (2–8 keV)
shown in purple on the right (total of 126 counts). In general, we find that
the hardness ratio of the primary point source is consistent with an AGN
($HR=0.37^{+0.01}_{-0.02}$, while the secondary point source has a softer
spectrum (albeit, with larger error bars; $HR=-0.34^{+0.48}_{-0.60}$)

Extended Data Figure 8: Spectral fits to the primary (left) and secondary
(right) X-ray point source detected by BAYMAX. We fit and model 100 spectral
realizations for each X-ray point source (blue) using XSPEC to quantify the
average values of the spectral shape and X-ray luminosity. We show one of
these spectral realizations in black, for each X-ray point source. We find
both sets of spectral realizations are best-fit with a simple absorbed power
law (XSPEC model components: phabs $\times$ zphabs $\times$ zpow). Given the
low number of counts for both components (average of $\sim 191.5$ and $\sim 6$
counts for the primary and secondary, respectively), we fix the power law
index, $\Gamma$, to a value of 1.8 for all spectral fits. We list the median
values and the confidence interval of the best-fit spectral parameter
distributions in Extended Data Table 3.  Extended Data Figure 9: Multi-
wavelength Spectral Energy Distributions (SEDs) of the SW (top) and NE nuclei
(bottom) in SDSS J0749$+$2255\. For comparison, the SEDs of typical optically
selected SDSS quasars[5] are shown in grey, with the mean values shown in
solid (“All”), and the optically luminous (“OL”) and dim (“OD”) sub-
populations shown in dotted curves. Also shown are the SEDs of optically
bright, non-blazar quasars from [53], including both radio-loud (“RL”; shown
in dash-dotted curves) and radio-quiet (“RQ”; shown in dashed curves)
populations. All comparison SEDs are normalized to match the luminosities of
SDSS J0749$+$2255 in HST F814W.  Extended Data Figure 10: Top: HST F160W and
Keck AO $K_{p}$ PSF-subtracted host-only photometry and the best-fit composite
stellar population (CSP) models for the NE (in cyan) and SW (in orange) nuclei
of SDSS J0749$+$2255\. The thick curves represent the best-fit CSP models
whereas the thin ones are 100 random draws from 500k MCMC realizations (with
1k walkers and 1k steps with 500 burn-ins) for illustration. Bottom:
Transmission curves for the HST F160W (dotted) and Keck $K_{p}$ (dashed)
filters, which roughly sample the rest-frame $g$ and $R$ bands for SDSS
J0749$+$2255. Extended Data Figure 11: Black hole mass versus host-galaxy
total stellar mass. The two nuclei in SDSS J0749$+$2255 are plotted as blue
(SW) and orange (NE) stars with the error bars representing 1$\sigma$
uncertainties. Three local scaling relations are shown with a red solid
line[59], a green dashed line[60], and a yellow dash-dotted line[61]. The
shaded regions denote the 1$\sigma$ scatters for these relations. Optical
selected quasars at z$>$0.2 with HST IR images are shown as dark gray
hexagons[64] and purple diamonds[66]. X-ray selected AGN are shown as light
grey squares[65] and light pink pentagons[62]. Quasars whose stellar masses
are estimated from broad-band SED fitting (i.e., rather than spatially
decomposing the PSF and the host) are shown as light brown circles[63].
Extended Data Figure 12: $M_{\rm BH}-M_{\rm gal}$ relation for the eight
simulated bright quasar pairs in Astrid. The colored lines mark the two pairs
(PairA and PairB) with apparently distinct galaxy bulges, although the
kinematic bulge-disk classification adopted by simulations cannot be reliably
used in interacting systems whose morphology is typically disturbed and the
disks are in the process of being disrupted in the merger process. PairB is
the closest match with SDSS J0749$+$2255 in both the SMBH masses and the host
galaxy stellar masses.  Extended Data Figure 13: A simulated quasar pair
(shown in red and cyan crosses) at $z\sim 2$ closely matching the properties
of SDSS J0749$+$2255, shown before and after the host galaxy merger. $\Delta
t$ refers to the time difference between the snapshot and the SMBH merger. The
second frame is similar to SDSS J0749$+$2255 at the time of observation, with
slightly lower SMBH and galaxy masses. The two SMBHs merge in the simulation
$83\,{\rm Myrs}$ after the “observation”, and the merger remnant has $M_{\rm
BH}\sim 3\times 10^{9}\,M_{\odot}$.