# Varstrometry for Off-nucleus and Dual sub-Kpc AGN (VODKA): Investigating the
Nature of SDSS J0823$+$2418 at $z=1.81$: a Likely Lensed Quasar

Arran C. Gross Department of Astronomy, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA  Yu-Ching Chen Department of Astronomy,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Adi Foord
Kavli Institute of Particle Astrophysics and Cosmology, Stanford University,
Stanford, CA 94305, USA  Xin Liu Department of Astronomy, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA  National Center for
Supercomputing Applications, University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA  Yue Shen Department of Astronomy, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA  National Center for
Supercomputing Applications, University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA  Masamune Oguri Center for Frontier Science, Chiba
University, Chiba 263-8522, Japan  Department of Physics, Graduate School of
Science, Chiba University, Chiba 263-8522, Japan  Andy Goulding Department of
Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
Hsiang-Chih Hwang School of Natural Sciences, Institute for Advanced Study,
Princeton, 1 Einstein Drive, NJ 08540, USA  Nadia L. Zakamska Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Yilun Ma (马逸伦) Department of Astrophysical Sciences, Princeton University,
Princeton, NJ 08544, USA  Liam Nolan Department of Astronomy, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

Dual quasars at small physical separations are an important precursor phase of
galaxy mergers, ultimately leading to the coalescence of the two supermassive
black holes. Starting from a sample of dual/lensed quasar candidates
discovered using astrometric jitter in Gaia data, we present a pilot case
study of one of the most promising yet puzzling candidate dual quasars at
cosmic noon ($z\sim 1.8$). Using multi-wavelength imaging and spectroscopy
from X-ray to radio, we test whether the SDSS J0823$+$2418 system is two
individual quasars in a bound pair at separation$\sim$0.64″, or instead a
single quasar being gravitationally lensed by a foreground galaxy. We find
consistent flux ratios ($\sim$1.25$-$1.45) between the two sources in optical,
NIR, UV, and radio, and thus similar spectral energy distributions, suggesting
a strong lensing scenario. However, differences in the radio spectral index,
as well as changing X-ray fluxes, hint at either a dual quasar with otherwise
nearly identical properties, or perhaps lensing-based time lag of $\sim$3 days
paired with intrinsic variability. We find with lens mass modeling that the
relative NIR positions and magnitudes of the two quasars and a marginally
detected central galaxy are consistent with strong lensing. Archival SDSS
spectra likewise suggest a foreground absorber via Mg ii absorption lines. We
conclude that SDSS J0823$+$2418 is likely a lensed quasar, and therefore that
the VODKA sample contains a population of these lensed systems (perhaps as
high as 50%) as well as dual quasars.

black hole physics — galaxies: active — quasars: general — surveys

##  1 Introduction

Our current understanding of galaxy evolution indicates that the massive
galaxies observed in the present-day universe have been built up in a series
of major mergers of less massive galaxies (Di Matteo et al., 2005). It is
believed that nearly all massive galaxies contain a central supermassive black
hole (SMBH) (Kormendy & Richstone, 1995), which grows primarily through
accretion of mass (Yu & Tremaine, 2002), and whose evolution is linked to the
growth of its host galaxy’s stellar bulge (yielding the $M_{\rm
BH}-\sigma_{\star}$ relation, McConnell & Ma, 2013). Galaxy mergers are often
invoked to explain the growth necessary to produce observed massive SMBHs
($e.g.,$ Zhang et al., 2021). Such phases of accretion via gas inflows to the
SMBH also trigger periods of immense activity (e.g., Hernquist & Mihos, 1995;
Wassenhove et al., 2012; Capelo & Dotti, 2016), seen observationally as an
active galactic nucleus (AGN), the most extreme of which are classified as
quasars ($L_{\rm bol}\gtrsim 10^{45}$erg s-1).

It is predicted that at small physical separations ($\lesssim$10 kpc) in the
later stages of a galactic merger (e.g, Wassenhove et al., 2012; Rosas-Guevara
et al., 2019), both SMBHs in the constituent galaxies have a significant
chance of being active simultaneously as dual quasars. The pair will
eventually form a more compact binary SMBH, ending in a BH merger event
(Begelman et al., 1980). Identifying dual quasar systems at small separations
is thus essential for constraining the rate of binary black hole mergers,
which in turn is needed to inform predictions of gravitational waves emitted
as a consequence of those mergers (e.g., Peters, 1964; Dotti et al., 2012;
Abbott et al., 2016; Holgado & Ricker, 2019; Goulding et al., 2019). The
activity of galaxy mergers and quasars is thought to peak around the epoch of
cosmic noon ($z\sim 2$) (Richards et al., 2006a). However, despite decades of
searches, both serendipitous and systematic discoveries of dual quasars have
remained fairly rare (see, e.g., Satyapal et al., 2017; De Rosa et al., 2019;
Gross et al., 2023, and references therein), yielding even fewer confirmed
dual AGN at the high redshift and small separations that are most relevant for
probing cosmic growth (see Chen et al., 2022, for a recent review). The lack
of observed dual quasars is partially due to their apparent intrinsic rarity.
Using Gaia EDR3-resolved pairs, Shen et al. (2023) estimate that the fraction
of double quasars (3 kpc$<$sep$<$30 kpc) among luminous unobscured quasars at
1.5$<z<$3.5 is $\sim 6.2\pm 0.5\times 10^{-4}$. Limitations of current
observatories to probe such small spatial scales in wide area surveys also
hamper efforts to identify promising candidates.

The confirmation of dual quasar candidates then requires multi-wavelength
follow-up observations, where systematic issues can lead to spurious
classifications. For example, X-ray observations can unambiguously classify
bright sources as AGN, but the low spatial resolution of X-ray telescopes
($\gtrsim$0.5″ for on-axis observations with $Chandra$) precludes identifying
dual sources with small separations. While IR surveys using the Wide-field
Infrared Survey Explorer (WISE; Wright et al., 2010) have had much success in
systematically recovering the obscured AGN population (Barrows et al., 2021),
resolution limitations also hamper robust detection of dual quasars with
separations on the kiloparsec scale (e.g., the $W2$ band has a point spread
function (PSF) full width at half maximum (FWHM) of 6″, corresponding to a
separation of $\sim$50 kpc at $z=2$).

A contaminant population for the dual quasar search is that of small-
separation gravitationally lensed single quasars. At high redshift, lensing
galaxies can be difficult to detect due to limited spatial resolution, as well
as their inherent faintness compared to the more luminous and potentially
blended quasar images. Indeed, many systematic searches for dual quasars have
turned up lensed quasars incidentally (see Shajib et al., 2022, for a recent
review). These systems are intriguing in their own right; cases of strong
gravitational lensing have been used in recent decades to address cosmological
questions including the value of the Hubble constant (Kochanek et al., 2006).
In particular, lensed quasars are novel probes of dark matter, and high-
resolution imaging of lensed quasar radio jets have been used to constrain the
so-called “fuzzy” dark matter model (Hui et al., 2017; Powell et al., 2023,
and references therein). Accurately determining the nature of a given dual
quasar candidate as a bona fide dual AGN or a single lensed quasar is
therefore critically important to its astrophysical implications.

Throughout this series of papers, we have endeavored to probe the hitherto
poorly explored regime of high-$z$ dual quasar candidates at sub-arcsec
separations through Varstrometry for Off-nucleus and Dual sub-Kpc AGN (VODKA).
The technique of varstrometry was first described by Shen et al. (2019) and
Hwang et al. (2020), hereafter Paper I and Paper II, respectively. Briefly,
varstrometry relies on variable jitters in the astrometric centroid position
of an unresolved source over a series of observations. Such a signature can
indicate the presence of two separate sources with intrinsic photometric
variability, making it possible to compile samples of candidate dual quasars
at kiloparsec-scale physical separations.

In Paper II we compiled a sample of candidate dual quasars by primarily
applying the varstrometry selection to the Gaia Data Release 2 (DR2) catalog
(Gaia Collaboration et al., 2018), as well as searching for resolved pairs.
Using a series of selection criteria based on astrometric noise and cross-
matches to SDSS, WISE, and Pan-STARRS, a catalog of 150 targets was assembled
in Chen et al. (2022), hereafter Paper III. From this catalog 84 targets were
subsequently observed in Paper III as part of follow-up two-band Hubble Space
Telescope ($HST$) imaging (SNAP-15900, PI: Hwang) to obtain optical
photometric colors, yielding a sample of VODKA targets composed of 45 resolved
pairs in Gaia with separations 0.4″$<$sep$<$0.7″ (limited by the minimal
separation in Gaia DR2; Arenou et al., 2018). Recently, we conducted a pilot
study in Chen et al. (2023), hereafter Paper IV, focused on one of these dual
quasar candidates, SDSS J0749+2255, first reported in Shen et al. (2021). A
robust combination of multi-wavelength imaging and spectroscopy confirmed this
system to be a dual quasar. The system is composed of two $\sim
10^{9}M_{\odot}$ SMBHs hosted by a disk-disk galaxy merger with tidal features
and stellar masses of $\sim 10^{11.5}M_{\odot}$, making it the first discovery
of such a system at nuclear separation of $\lesssim$4 kpc at cosmic noon.

In this paper, we build on the previous work by performing a case study on
another promising VODKA target, SDSS J082341.08$+$241805.0, hereafter SDSS
J0823$+$2418\. Our main goal in this work is to determine whether the SDSS
J0823$+$2418 system is a dual quasar or a lensed single quasar, and to
characterize its properties. SDSS J0823$+$2418 is spatially resolved by Gaia,
and received preliminary follow-up imaging in Paper III using $HST$ WFC3 bands
F475W and F814W, shown in Figure 1. The difference in photometric colors
suggest that both the North (N) and South (S) components are quasars, and not
foreground stars (see Paper III for selection criteria). We also consider SDSS
J0823$+$2418 a promising candidate dual quasar because it exhibits two
distinct detections in preliminary radio imaging, indicative of two radio AGN
with different spectral indices. If confirmed as a dual quasar, SDSS
J0823$+$2418 would be part of a select population of high redshift ($z\sim
1.81$) systems at an advanced merger stage (sep$\sim 5.4$ kpc).

Figure 1: HST/WFC3 color composite image from Paper III, where F475W is coded
in blue, F814W in red, and the average of F475W and F814W bands in green). The
Gaia DR2 source positions are marked by the green crosses. The angular and
physical separations between the double sources are given in the bottom right
corner.

We organize the paper as follows. In §2 we detail the multi-wavelength
observations and their associated data reduction procedures. In §3 we perform
various analyses and give the general interpretation case by case. We combine
all of the results as lines of evidence in §4 to weave together a coherent
picture of SDSS J0823$+$2418\. We also compare it to similar systems from the
literature. In §5 we summarize our findings and offer concluding remarks for
future campaigns with similar strategies. Throughout this work, we assume a
flat $\Lambda$CDM cosmology with values of $\Omega_{\rm\Lambda}$=0.7,
$\Omega_{\rm m}$=0.3, and $H_{0}$=70 km s-1 Mpc-1.

##  2 Observations and Analysis

To assess whether the SDSS J0823$+$2418 system is a dual quasar, we employ
multi-wavelength observations ranging from X-rays to radio, using a mix of
photometry and spectroscopic observations. In some wavebands, we have obtained
complementary observations to account for various differences in observational
performance, i.e., higher angular resolution versus deeper observations for
faint structure.

In addition to the observations detailed below, we have also acquired slit
spectroscopy observations using the Magellan Folded-port Infrared Echellette
(FIRE) spectrograph. However, the redshift of SDSS J0823$+$2418 places the
prominent rest-frame optical emission lines (H$\alpha$, [O iii], and H$\beta$)
in windows that are heavily affected by telluric absorption, even after
correction. The traces of the 2D spectra are only marginally separated,
complicating the extraction. While the broad H$\gamma$ emission line does show
some differences between the N and S sources, there might be some
contamination from Fe ii emission. We do not discuss this observation further
owing to these complications.

###  2.1 HST Imaging

Infrared (IR) imaging in the observed frame has the potential to uncover low-
surface brightness features critical for the interpretation of the SDSS
J0823$+$2418 system. For example, detection of faint tidal features, such as
trailing galactic arms, would be highly suggestive of an ongoing galactic
merger between the two nuclei. Conversely, detection of a faint central galaxy
between the two nuclei could be evidence of gravitational lensing. Both of
these scenarios require careful subtraction of the bright point source quasars
to uncover the fainter residual features.

We obtained high-resolution imaging of SDSS J0823$+$2418 using the HST Wide
Field Camera 3 (WFC3) on 8 February 2022 (Program GO-16892; PI: X. Liu). We
use the NIR range F160W filter to capture the faint underlying hosts. This
filter operates in the wide (width 168.3 nm) $H$-band ($\lambda_{\rm pivot}$ =
1536.9 nm), corresponding to a wavelength of $\lambda\sim 546$ nm in the rest-
frame of the host galaxy. We obtain a net exposure time of 2055 s. The
detector’s PSF (FWHM = 0.15″) is undersampled by the pixel scale (0.13″)
necessitating dithering during the observations. The dithered frames are
stacked during data reduction following standard routines for WFC3 (Sahu,
2021), which clean the frames for cosmic ray trails and pixel area effects. We
show the processed F160W image in the leftmost panel of Figure 5. The reduced
image has a pixel scale of 0.06″ and a photometric zeropoint of 25.936
(28.177) in the AB(ST) system. We reach a surface brightness limit of
$\sim$25.65 mag/arcsec2.

###  2.2 Keck Imaging

We pursue a similar approach as with $HST$ using NIR imaging from the Keck
Telescope under Program N072 (PI: Y. Shen). The adaptive optics (AO)-assisted
NIRC2 camera is able to achieve a finer resolution than $HST$, with a PSF FWHM
of 0.05″ sampled at 0.01″/pix. This enables us to better model the point
sources and extended host components so that a central lens galaxy can be more
accurately constrained.

Observations of SDSS J0823$+$2418 were conducted using the Keck NIRC2 camera
on 18 December 2021 UT. We observed the target using the NIR $K_{p}$ filter,
which has a central wavelength of 2.124 $\mu$m and width of 0.351 $\mu$m. The
nearby star N910000396, 41.68″ from SDSS J0823$+$2418, was used as a tip-tilt
reference star. Dithering was done using a 4″$\times$4″ 9-position box. After
standard reduction procedures, we co-add the exposures for a total exposure
time of 1260 s. The Laser Star Guide mode was used for the adaptive optics
corrections for atmospheric distortions during observations, and we
subsequently use the guide stars to compute the photometric zeropoint of the
$K_{p}$ filter to be 26.56 (AB). We reach a surface brightness limit of
$\sim$24.94 mag/arcsec2. We show the reduced $K_{p}$ image in the leftmost
panel of Figure 6.

###  2.3 HST/STIS Spectroscopy

Observations of SDSS J0823$+$2418 using the $HST$/Space Telescope Imaging
Spectrograph (STIS) were obtained on 13 January 2022 (Program GO-16210; PI: X.
Liu). A 0.2″ slit was used for a net exposure time of 897 s. at a position
angle of 48.6∘ E of N to cover both nuclei simultaneously. The 5030Å bandwidth
ranges from 5236Å to 10266Å, with a spectral resolution of $\sim$790 at 7751Å.
The spatial resolution of PSF FWHM$\sim$0.08″ allows for a reasonable
separation of the two nuclei.

Figure 2: STIS spectral extraction. Top: 2D spectrum with the extraction
regions marked in white. We highlight the position of the bright Mg ii
emission line. Bottom: Collapsed 1D profiles and corresponding extraction
regions.

We reduce the spectra using the stis_cti reduction package (Anderson & Bedin,
2010), accounting for flat-fielding, combining exposures, and removing trails
from cosmic rays and other artefacts. In Figure 2, we show the 2D spectrum
highlighting the separation of the traces for the two nuclei, as well as the
collapsed 1D brightness profiles and the extraction apertures. At an angular
separation of 0.64″, we are able to obtain reasonable extractions for the two
nuclei using apertures with widths of 5 pixels, corresponding to 0.25″, shown
in Figure 2. After wavelength and flux calibration, we obtain the final 1D
spectra for each nucleus, shown in Figure 3.

Figure 3: Observed-frame optical spectra for the two quasars. Top: We show the
HST/STIS spectra for SDSS J0823$+$2418 as well as the Gemini/GMOS spectra,
indicated by the legend. The STIS spectra have been smoothed with a Gaussian
kernel with sigma = 1 for readability. Green shaded regions are heavily
affected in the GMOS spectra by telluric absorption. Bottom: Flux ratios
between the N and S quasars for the STIS and GMOS spectra. Dashed lines show
the average of the corresponding flux ratio, with values that are roughly
consistent with those found with photometry. The N and S spectra from both
instruments are quite similar, suggestive that the system is a single lensed
quasar where the N source is more magnified than the S component.

###  2.4 Gemini/GMOS Spectroscopy

As an independent verification of the STIS spectra, we obtained optical/NIR
spectra for SDSS J0823$+$2418 using the Gemini Multi-Object Spectrograph
(GMOS) on the Gemini North Telescope on the night of 20 April, 2022 (Program
ID GN-2022A-Q-139; PI: X. Liu). GMOS has wavelength coverage that goes further
into the blue (5000Å-10000Å) than STIS, allowing us to fully cover the Ciii]
quasar emission line (rest-wavelength $\lambda$1909Å). We used the R150
grating and 0.75″slit oriented at a position angle of 3.2∘ E of N to cover
both nuclei. The spectral resolution of GMOS in this configuration
($R\sim$631) is slightly lower than that of the STIS observation above, as is
the spatial resolution with PSF 0.39″$<$FWHM$<$0.49″. However, the deeper
observations allow for higher signal to noise compared to the STIS data.

Figure 4: GMOS spectral extraction. Top: 2D spectrum with the extraction
regions marked in white. We highlight the position of the bright Mg ii
emission line. Bottom: Collapsed 1D profiles and corresponding extraction
regions.

We reduce the raw spectrum using the pypeit spectral processing package
(Prochaska et al., 2020a, b), set to gemini_gmos_north_ham mode. The pipeline
handles the typical reduction steps correcting for flat-fields, bias, and
cosmic rays. We manually set apertures with widths of 7.43 pixels
(corresponding to 0.6″) for boxcar extraction of the 1D spectrum for each
nucleus. In Figure 4, we show the 2D spectrum, as well as the collapsed 1D
brightness profiles and the extraction apertures. Extraction is less
straightforward with the GMOS spectrum than the STIS spectrum above due to the
larger FWHM of the PSF. Thus, there is some degree of source blending in the
extraction apertures. Wavelength calibration is done using the CuAr spectra
acquired during observations. Flux calibration is done using observations of
the standard star EG 131. We do not use pypeit telluric corrections as they
are not yet finalized in the pipeline; we thus mask heavily affected regions
in our subsequent treatment. We show the two calibrated spectra from GMOS in
Figure 3.

###  2.5 VLA Imaging

Table 1: VLA Observations

Frequency | UT | IT | rms | Beam | PA
---|---|---|---|---|---
Band | date | min |  $\mu$Jy/beam | (arcsec) | (deg)
(1) | (2) | (3) | (4) | (5) | (6)
C (6 GHz) | 2020-12-13 | 30 | 8.4189 | 0.92$\times$0.35 |  $-$69
Ku (15 GHz) | 2020-12-10 | 30 | 7.3683 | 0.11$\times$0.10 |  $+$18

Note. — (1) Frequency band and nominal central frequency; (2) UT date of
observation (Y-M-D); (3) Integration time for the observation, encompassing
the kinematic pair; (4) rms noise for the cleaned image; (5) Restoring beam
size (maj$\times$min axes); (6) Restoring beam position angle.

Imaging SDSS J0823$+$2418 at radio frequencies has several advantages that
make it complimentary to optical and IR observations. The spectral slope of
the radio emission can be used as a tracer of synchrotron emission from jet
activity, and differences between the spectral index values of the two sources
could indicate a system of two distinct quasars instead of a single lensed
source.

We obtained radio imaging of SDSS J0823$+$2418 using the NSF’s Karl G. Jansky
Very Large Array (VLA) in two frequency bands: C-band (central frequency
5.9985 GHz, width of $\sim$4 GHz) and Ku-band (central frequency 14.9984 GHz,
width of 6 GHz). Observations were conducted as part of Program 20B-242 (PI:
X. Liu), and both bands were observed in the VLA’s A-configuration. The source
3C147 was also observed and used for flux and bandpass calibration via the
standard VLA calibration and flagging pipeline. We use the Common Astronomical
Software Applications (CASA) package version 5.6.1 for data reduction and
analysis of the Stokes I images. We reduce the images using the synthesis
imaging and deconvolution routine tclean in w-project mode to account for sky
curvature. We use Briggs weighting with robustness of 0.5 to achieve a balance
between resolving small scale structure (i.e., separating the two nuclei)
while also preserving low surface brightness features. The pixel scales are
set to 0.1″ and 0.03″ for the C-band and Ku-band images, respectively, so that
3$-$5 pixels are sampled by the restoring beams. We set cleaning to run until
either 1000 iterations are completed, or the rms residuals drop below 12
$\mu$Jy beam-1. The properties of the restoring beams and noise levels for
each image are given in Table 1. In the leftmost panels of Figure 9 we show
the cleaned radio continuum images at 6 GHz and 15 GHz.

###  2.6 Chandra Imaging and Spectroscopy

Observations of SDSS J0823$+$2418 using the Chandra X-ray Observatory Advanced
CCD imaging spectrometer (ACIS; Garmire et al., 2003) -S array were split into
two observation blocks under Program 23700377 (PI: X. Liu; ObsIDs 25385 and
26279, hereafter referred to as observations 1 and 2, respectively). These
correspond to MJDs 21 and 24 of January 2022. Both observations were conducted
with the target close to the aim-point on the S3 chip, for durations of 3.83
and 14.86 ks, respectively. We reduce and analyze the observations using the
Chandra Interactive Analysis of Observations software (ciao) Version 4.14,
with calibrations from caldb 4.9.8.

The datasets are reprocessed using chandra_repro with the energy-dependent
sub-pixel event repositioning algorithm (EDSER; Li et al., 2004), subsampling
the pixels to 1/2 the native size. We then perform an astrometric correction
by detecting sources in the level-2 event file via wavdetect, which we then
cross-match to optical sources from the Sloan Digital Sky Survey Data Release
9 (SDSS DR9) catalog. For each observation, we require at least 3 matches to
within a tolerance of 2″. We find that background flaring is within 3$\sigma$
of the mean background level and is thus negligible for all time intervals. In
the top row of Figure 10, we show the reduced data for the two observations in
the broad 0.5$-$7 keV band.

##  3 Results

Differences in fluxes at a given wavelength between the N and S quasars would
be suggestive of two unique sources constituting a dual quasar system.
However, if the flux ratio is constant as a function of wavelength, then the
difference might be due to wavelength-independent lensing magnification
instead. In this section we test the imaging and spectroscopic data for such
differences. From Paper III, we have $HST$ photometric magnitudes for two
optical bands, which we list in Table 2.

Table 2: General Properties of SDSS J0823$+$2418

Target | Redshift | F475W | F814W | log${\it M}_{\rm BH,vir}$
---|---|---|---|---
|  | mag | mag | ${\it M}_{\odot}$
(1) | (2) | (3) | (4) | (5)
N | 1.81354$\pm 0.00036$ | 18.49$\pm$0.04 | 17.71$\pm$0.04 | 9.53$\pm$0.50
S | 1.81328$\pm 0.00036$ | 18.88$\pm$0.04 | 18.03$\pm$0.04 | 9.55$\pm$0.54

Note. — (1) Target SDSS J0823$+$2418; (2) redshift measured from best-fit STIS
Mg ii emission line center; (3)-(4) archival photometric magnitudes (AB) from
HST (Chen et al., 2022); (5) Black hole mass derived from virial mass
estimates using Mg ii.

###  3.1 IR Image Decomposition

To model the surface brightness profiles of each point source quasar, we first
generate a PSF model. The $HST$ PSF is known to vary during oribits due to
“breathing” ($e.g.,$ Lallo et al., 2005). Another complication is that our
F160W observation of SDSS J0823$+$2418 does not have any field stars in frame
to use as the PSF basis. We therefore construct an average PSF using 8 field
stars from 6 other VODKA observations that were taken within $\sim$1 month of
SDSS J0823$+$2418 and use this as our baseline PSF model.

Figure 5: $HST$ F160W imaging spatial decomposition with GALFIT. a: imaging
data. b: Model 3, where the legend gives the two quasar point sources and lens
galaxy Sérsic profile. c: the image with only the point sources subtracted,
leaving residuals due to the lensing galaxy but also the strong contribution
from the S host galaxy. d: the full model subtracted from the data, where the
residuals are scaled by $\sigma$ significance above the remaining background.
While the PSF models leave residuals corresponding to the host galaxies, the
central lensing galaxy is well constrained. e: radial profile starting from
the center of the N PSF. The lower sub-panel shows the scaled residuals as a
function of radius. Model 3 does a reasonably good job out to a radius of
$\sim$1″, beyond which the brightness from the extended host galaxies falls
off rapidly. The S host component closely follows the S PSF component, while
the N host component is heavily suppressed. Major tick marks in panels a-d are
0.5″.

We use the software galfit (Peng et al., 2010) to perform a 2D spatial
decomposition of the NIR image from $HST$. Our zeroth-order model, given as
Model 0 in Table 3, includes the two PSF components for the N and S sources
where the locations are allowed to vary, as well as a constant background. We
also include components for two sources that are $>$7″ to the NW of the N
quasar (just outside of the frame in Figure 5). The inclusion of these other
galaxies does not affect the best-fit parameters of the quasars, so we do not
discuss them further; however, their inclusion improves the overall fit
statistics, so they are included in all of our models going forward. These two
external sources are modeled with a simple Sérsic profile:

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

where $\Sigma(r)$ is the pixel surface brightness as a function of radial
distance $r$ from the center, $\Sigma(r_{\rm e})$ is the pixel surface
brightness at the “effective” half-light radius $r_{\rm e}$, $n$ is the Sérsic
index, and $\kappa$ is a parameter related to n via complete and incomplete
Gamma functions. The radius $r$ also encapsulates the ellipticity and position
angle of the profile.

The statistical uncertainties on the best-fit parameters are known to be
underestimated by galfit ($e.g.,$ Häussler et al., 2007), and the variations
in the $HST$ PSF likely contribute large systematic uncertainties. To more
accurately assess the PSF contribution, we perform the model fitting again
using each of the 8 individual field star PSF models. For each fitted
parameter, we then take the standard deviation across the 8 individual runs
and baseline model, and add this to the statistical uncertainties output by
galfit to quote the total uncertainties in Table 3.

Model 0 fails to capture the diffuse emission that extends $\sim$1.5″ out from
the center of the system. The residuals appear to be particularly concentrated
between the two point sources instead of showing a diffuse extended
morphology. The location of the residuals is suggestive of a faint central
galaxy between the two quasars. Under the supposition that the two quasar
targets might actually be a single quasar strongly lensed by a central
foreground galaxy located between the N and S source locations, we define a
more complex Model 1 which includes all of the components of Model 0 plus the
lens galaxy component. The brightness distribution of this lens galaxy
component is modeled using a Sérsic profile.

Table 3: HST (F160W) Spatial Modeling

Parameter | SDSS J0823$+$2418 N | SDSS J0823$+$2418 S
---|---|---
Model 0 (2 PSFs): $\chi^{2}_{\nu}$ = 7.777
mPSF (AB) | 17.929$\pm$0.002 | 18.209$\pm$0.002
CPSF | 08:23:41.09+24:18:05.5 | 08:23:41.09+24:18:04.9
Model 1 (2 PSFs + 1 Sérsic (lens)): $\chi^{2}_{\nu}$ = 4.357
mPSF (AB) | 17.960$\pm$0.001 | 18.297$\pm$0.001
CPSF | 08:23:41.09+24:18:05.5 | 08:23:41.09+24:18:04.9
mlens (AB) | 19.522$\pm$0.009
Clens | 08:23:41.09+24:18:05.1
$R_{\rm lens}$ (kpc) | 3.27$\pm$0.05
nlens | 1.93$\pm$0.05
Model 2 (2 PSFs + 2 Sérsics (hosts)): $\chi^{2}_{\nu}$ = 4.167
mPSF (AB) | 17.967$\pm$0.001 | 18.882$\pm$0.091
CPSF | 08:23:41.09+24:18:05.5 | 08:23:41.09+24:18:04.9
mhost (AB) | 20.400$\pm$0.021 | 18.853$\pm$0.094
$R_{\rm host}$ (kpc) | 5.06$\pm$0.11 | 0.51$\pm$0.07
nhost | 0.82$\pm$0.05 | 3.43$\pm$0.40
Model 3 (2 PSFs + 2 Sérsics (hosts) + 1 Sérsic (lens)): $\chi^{2}_{\nu}$ =
3.636
mPSF (AB) | 17.967$\pm$0.001 | 18.916$\pm$0.134
CPSF | 08:23:41.09+24:18:05.5 | 08:23:41.09+24:18:04.9
mhost (AB) | 21.632$\pm$0.100 | 19.034$\pm$0.151
$R_{\rm host}$ (kpc) | 5.06$\pm$0.26 | 0.6$\pm$0.1
nhost | 0.33$\pm$0.05 | 0.2$\pm$0.28
mlens (AB) | 19.95$\pm$0.021
Clens | 08:23:41.09+24:18:05.1
$R_{\rm lens}$ (kpc) | 3.94$\pm$0.10
nlens | 1.94$\pm$0.05

Note. — Best-fit components for 4 different models fit to the F160W image
using galfit. Central coordinates (C) of model components are in
hh:mm:ss.ss+dd:mm:ss.s; the host components are fixed to the same location as
the corresponding best-fit PSF. The lensing galaxy component, when present, is
denoted with lens and given as a column between the N and S components. The
best-fit magnitude, effective radius, and Sérsic index for a component are
given by m, $R$, and n, respectively.

The residuals of Model 1 suggest that the underlying host galaxies of the
quasars are contributing to the extended brightness distributions to some
degree, as expected at optical/NIR wavelengths. This is also reflected in the
radial profiles where the enclosed brightness is well modeled out to a radius
of $\sim$1″ (from the center of the N source), but is not well constrained at
larger radii. We therefore form Model 2 using the same parameters as in Model
0 with the addition of 2 more Sérsic profiles to model the underlying N and S
host galaxies. We tie the host Sérsic positions with the corresponding PSFs to
reduce model degeneracy (leaving the host locations as free parameters does
not improve the fit). We omit the lens galaxy component in Model 2. While this
model shows an improvement from the previous two, the residuals suggest that
even with the host galaxies included there is some amount of lingering light
between the two quasars in the form of a compact clump detected at
significance $>10\sigma$. The position of the clump seems consistent with the
best-fit location of the lens galaxy component in Model 1.

To address this clump, we form Model 3 as a combination of Models 1 and 2:
Model 3 includes a PSF and Sérsic profile for both sources plus a central lens
galaxy Sérsic, shown in Figure 5. This model achieves the lowest reduced
$\chi^{2}$ statistic out of the four. While the model still leaves some
residuals, they are not distributed in a coherent way warranting an additional
model component. The residuals also are not detected at more than 50% above
the noise level, so we do not consider them robust indications of merger-
induced tidal structures.

Figure 6: Same as Figure 5, but for Keck AO imaging. We again show Model 3.
The extended host components appear more prominent in the $K_{p}$-band, but do
not contribute much to the radial profile beyond a radius of 1″. The irregular
shape of the PSF is due to the adaptive optics of Keck. Major tick marks in
all panels are 0.25″. In panel $c$, the faint blob above the S PSF could be a
faint foreground galaxy. There also appears to be an arc to the residual
structure of the N galaxy that is mirrored in the S galaxy. This overall
ellipsoid arrangement suggests some degree of lensing shear.

The best-fit model (Model 3) characterizes the N quasar as markedly brighter
than its S counterpart, with a flux ratio of North to South of $\frac{f_{\rm
N}}{f_{\rm S}}\sim 2.4\pm 0.4$. The fit results paint a different picture for
the host galaxies, where the N host galaxy is substantially fainter than the S
host (flux ratio of 0.09$\pm 0.07$). The N host is much more spatially
extended than the S host, which is at the lower limit for the fit parameter
such that it functions like a PSF profile. In the radial profiles for both
Models 2 or 3, the S host galaxy component appears to mimic the profile of the
PSF which results in suppressing the contribution of the PSF, and thus a
higher N/S flux ratio. In Model 3 particularly, the S host radius and Sérsic
index approach the minimum allowable values for the fit, which are unphysical.
This is likely due to model degeneracy between the PSF and Sérsic components
in the fainter S source, especially when the nearby lens component is added.
As a check of whether the S host galaxy is indeed detected, we run the fitting
without its inclusion resulting in a substantially poorer fit
($\chi^{2}_{\nu}=5.496$) than most of the other models. We consider the S host
galaxy a necessary model component, but limitations of the fitting prevent a
detailed comparison of the host galaxy properties. We therefore also consider
the results of Model 1 when discussing the flux ratio of the N and S quasars;
we obtain a value of $\frac{f_{\rm N}}{f_{\rm S}}\sim 1.36\pm 0.06$ which is
consistent with the previously found flux ratios for photometry using WFC3
filters F475W (1.43$\pm$0.07) and F814W (1.34$\pm$0.07). At the spatial
resolution of $HST$, it is difficult to say definitively whether the centrally
located brightness distribution is due to a central lens galaxy or the
extended underlying galaxies. However, based on the penalized log-likelihood
ratio test of the Akaike information criterion (AIC; $e.g.$, Akaike, 1974;
Ruffio et al., 2019), Model 3 is strongly statistically preferred to both
Models 1 and 2, indicating that all of the components are warranted.

Table 4: Keck ($K_{p}$) Spatial Modeling

Parameter | SDSS J0823$+$2418 N | SDSS J0823$+$2418 S
---|---|---
Model 0 (2 PSFs): $\chi^{2}_{\nu}$ = 3.581
mPSF (AB) | 18.413$\pm$0.002 | 18.644$\pm$0.002
CPSF | 08:23:41.13+24:18:06.0 | 08:23:41.10+24:18:05.4
Model 1 (2 PSFs + 1 Sérsic (lens)): $\chi^{2}_{\nu}$ = 1.295
mPSF (AB) | 18.484$\pm$0.001 | 18.749$\pm$0.001
CPSF | 08:23:41.13+24:18:06.0 | 08:23:41.10+24:18:05.4
mlens (AB) | 17.985$\pm$0.004
Clens | 08:23:41.11+24:18:05.7
$R_{\rm lens}$ (kpc) | 4.27$\pm$0.03
nlens | 0.46$\pm$0.01
Model 2 (2 PSFs + 2 Sérsics (hosts)): $\chi^{2}_{\nu}$ = 1.211
mPSF (AB) | 18.525$\pm$0.001 | 18.808$\pm$0.002
CPSF | 08:23:41.13+24:18:06.0 | 08:23:41.1o+24:18:05.4
mhost (AB) | 19.324$\pm$0.008 | 17.842$\pm$0.024
$R_{\rm host}$ (kpc) | 1.91$\pm$0.02 | 5.87$\pm$0.16
nhost | 0.51$\pm$0.02 | 1.65$\pm$0.04
Model 3 (2 PSFs + 2 Sérsics (hosts) + 1 Sérsic (lens)): $\chi^{2}_{\nu}$ =
1.106
mPSF (AB) | 18.524$\pm$0.001 | 18.795$\pm$0.001
CPSF | 08:23:41.13+24:18:06.0 | 08:23:41.10+24:18:05.4
mhost (AB) | 19.628$\pm$0.011 | 17.925$\pm$0.024
$R_{\rm host}$ (kpc) | 2.25$\pm$0.02 | 6.82$\pm$0.21
nhost | 0.29$\pm$0.02 | 1.35$\pm$0.04
mlens (AB) | 19.877$\pm$0.014
Clens | 08:23:41.11+24:18:05.7
$R_{\rm lens}$ (kpc) | 3.19$\pm$0.04
nlens | 0.20$\pm$0.02

Note. — Best-fit components for 4 different models fit to the $K_{p}$ image
using galfit. Central coordinates (C) of model components are in
hh:mm:ss.ss+dd:mm:ss.s; the host components are fixed to the same location as
the corresponding best-fit PSF. The lensing galaxy component, when present, is
denoted with lens and given as a column between the N and S components. The
best-fit magnitude, effective radius, and Sérsic index for a component are
given by m, $R$, and n, respectively.

We also model the Keck $K_{p}$-band AO image using galfit. For consistency, we
adopt the same model components as we used for the F160W image above. The
results of the 4 model fits are given in Table 4. We show results of Model 3
in Figure 6. Here we see stronger indications of underlying disk structure
because of the higher S/N compared to the HST image. The disks appear to be
aligned in a parallel geometry which suggests a lensing scenario. As with the
F160W image, Model 3 (2 PSF point sources, 2 host galaxies, and a lens galaxy)
achieves the lowest reduced $\chi^{2}$ value ($\sim$1.1); however, we do not
see the same degeneracy between the S source’s PSF and Sérsic components here,
as was the case above. We do not detect any faint extended structure in the
residual image suggestive of tidal tails. For the quasars, we find a
qualitatively similar result to the F160W fit, where the N quasar is brighter
than the S quasar. The flux ratio of 1.284$\pm 0.002$ for the $K_{p}$-band
quasars is roughly consistent with the values above. Here, we are quoting the
flux ratio derived from Model 3, although the value derived from Model 1
(1.276$\pm$0.002) is also roughly consistent. The flux ratio of the host
galaxies is less extreme (0.208$\pm 0.025$) than found for the $HST$ image.
However, the best-fit sizes of host components appear to be reversed from the
F160W fitting above, such that the S host here is much more extended (beyond
the physical separation of the two nuclei) compared to the N host; this once
again suggests the limitations of our spatial analysis with such a blended
pair and model degeneracy.

Because of the lower sensitivity of Keck compared to $HST$, it is not as
obvious in panel $c$ of Figure 6 whether there is in fact a cluster of bright
residuals constituting a central lens galaxy; however, there appears to be a
faint blob detected at $\sim 4\sigma$ significance towards the S source that
is spatially consistent with the location of the lens galaxy in the $HST$
image modeling, and its inclusion is statistically warranted. As with the HST
modeling, the AIC values indicate that Model 3 is strongly statistically
preferred over Models 1 and 2, bolstering the need for including both the
lensing galaxy and the host galaxy components.

###  3.2 Optical/IR Spectral Fitting

At first glance, the N and S components appear to be similar, with a nearly
constant flux ratio. The flux calibration of the GMOS spectra yields
systematically lower fluxes compared to the STIS spectra. Since both nuclei
appear equally affected, it does not influence our bulk results of comparing
the spectra for differences. To check for differences, we normalize the GMOS
spectrum of the S source by the average flux ratio observed between the two
sources (shown with dashed line in Figure 3), and then take the difference to
reveal any obvious unique features. The GMOS spectra are remarkably similar,
apart from several isolated features, some of which are artefacts. On the
other hand, the STIS spectra show several differences between the N and S
sources. Subtle variations in the broad shapes of the STIS spectra might be a
result of the poor S/N. In particular, the STIS-N spectrum appears to exhibit
two strong narrow emission lines around $\lambda_{\rm obs}=6895{\rm\ and\
}6715$ ÅÅ. Given the redshift of the object, these two features do not
correspond to any typical strong quasar emission lines. Considering that we
see a similar feature in the STIS-S spectrum at $\sim 9200$Å, we believe these
features are due to uncleaned cosmic rays, and so we have masked them in
Figure 3, as well as with gray regions in Figure 7. The STIS-S spectrum
appears to have absorption features not present in the N spectrum. These
disparate features could be indications that the two spectra originate from
two non-lensed quasar cores. However, we do not see evidence of these same
features in the corresponding GMOS spectra, which is much less noisy than the
STIS observations. An important caveat here is that the GMOS extraction
regions are not as well separated as those used for STIS as a consequence of
the PSF and seeing, so it is difficult to conclude whether those features are
truly absent in the GMOS spectra or have been washed out due to some degree of
blending between both sources in each individual spectrum.

To constrain the properties of the quasars’ emission lines and the underlying
continua, we employ the spectral fitting code PyQSOFit (Guo et al., 2018). We
model each spectrum separately because of differences in spectral resolution.
We mask out regions with persistent artefacts of the spectral reduction. The
spectra are first corrected for extinction using the Milky Way extinction law
of Cardelli et al. (1989) and drawing upon the dust maps from Schlegel et al.
(1998), and then de-redshifted to rest-frame. We model each galaxy as a
pseudo-continuum that is a linear combination of a power law, Fe ii emission,
and a low-order polynomial for quasar continuum, combined with Gaussian
components for broad and narrow emission lines. Local fits are then done for
the emission line regions after subtracting the continuum elements. We show
the spectral fits for the STIS and GMOS data in Figures 7 and 8, respectively.

While we obtain reasonable fits for the GMOS spectra, there are several
windows where residual telluric absorption features need to masked out; one of
these windows truncates the broad Mg ii emission line, and so we rely
predominantly on the fit results of the STIS spectra for this emission line
fit. Uncertainties in fitted parameters are calculated using Markov Chain
Monte Carlo simulations. We determine the systematic redshifts of the nuclei
based on the velocity offsets observed for the emission lines. The best-fit
redshifts are 1.81354$\pm 0.00036$ and 1.81328$\pm 0.00036$ for the N and S
nuclei, respectively, which agree within the errors. We determine the slopes
of the underlying continua to be -1.030$\pm$0.014 and -1.031$\pm$0.005 for the
N and S sources. The slopes are consistent, suggesting that the quasars are
the same object doubly imaged.

Figure 7: Spectral fitting of the $HST$/STIS spectra for SDSS J0823$+$2418\.
The top and bottom sets of panels show the fits for the N and S targets,
respectively. The spectra are fit with a pseudo-continuum model composed of a
power law combined with a low-order polynomial (orange curve) and the Fe ii
template. Gray regions are masked out during the fit due to lingering cosmic
rays. The lower inset panels show zoom-ins on the local fits of broad emission
lines after continuum subtraction.

Figure 8: Same as Figure 7, but spectral fitting of the Gemini/GMOS spectra.
Gray regions are masked out primarily to avoid regions that are heavily
affected by telluric absorption bands. The extended range towards bluer
wavelengths allows for a better fit of the Ciii] emission line than is
possible with the STIS spectra. Table 5: Emission Line Properties of SDSS
J0823$+$2418

Emission Line | SDSS J0823$+$2418 N | SDSS J0823$+$2418 S
---|---|---
Measurement |  |
Ciii] Flux | 692.1 $\pm$ 8.2 | 688.1 $\pm$ 33.1
Mg ii Flux | 1433.8 $\pm$ 59.7 | 949.8 $\pm$ 58.4
Ciii] EW | 17.0 $\pm$ 0.2 | 21.0 $\pm$ 1.0
Mg ii EW | 34.6 $\pm$ 1.4 | 29.3 $\pm$ 1.8
Ciii] FWHM | 7142.1 $\pm$ 99.5 | 7424.8 $\pm$ 162.6
Mg ii FWHM | 5400.4 $\pm$ 268.4 | 5871.2 $\pm$ 560.9

Note. — Emission line properties measured from the spectral fitting of the
$HST$/STIS and Gemini/GMOS spectra in the observed-frame optical. Flux is
reported in units of 10-17 erg s-1 cm-2, rest-frame equivalent width (EW) is
in units of Å, full width at half max (FWHM) is in units of km s-1. All values
are for the broad component of the line. All errors are quoted at the
1-$\sigma$ level from Monte Carlo simulations.

The resulting best-fit parameters of the UV emission lines are given in Table
5. The fluxes of the Mg ii lines are different between the two sources by
$\sim$50%, which is again a roughly consistent flux ratio to other wavebands.
The equivalent width and FWHM values are roughly consistent. Because the Ciii]
emission line is not fully covered by the STIS spectra, we rely on the GMOS
spectral fitting for measurements in this case. Again, we find that the line
properties are roughly consistent between the N and S nuclei. The flux ratio
of the continuum at rest-frame 2500 Å is $\sim$1.35 which is consistent with
the results above from photometry.

We estimate the black hole mass for each nucleus using the empirical mass
estimator from Vestergaard & Osmer (2009) using Mg ii emission:

$M_{\rm BH}=10^{6.86}\left[\frac{{\rm FWHM(MgII)}}{{\rm 1000\ km\
s^{-1}}}\right]^{2}\left[\frac{{\rm\lambda}L_{\rm\lambda}(3000\
{\rm\AA})}{10^{44}\ {\rm erg\ s^{-1}}}\right]^{0.5}\ M_{\odot},$ (2)

where ${\rm\lambda}L_{\rm\lambda}$ is the monochromatic luminosity taken at
3000 Å obtained from our continuum fitting. We list the resulting black hole
masses in Table 2. Owing to the similarity of the spectra, the derived masses
are consistent with each other. The $\sim$50% differences in Mg ii flux, which
could be due to lensing magnification, do not impact the derived masses as
much as the FWHM values, which are consistent to within 9%.

Table 6: Radio Properties of SDSS J0823$+$2418

Radio Designation | $S^{\rm peak}_{\nu}$ | $S^{\rm int}_{\nu}$ | log $\nu L_{\nu}$ | Maj | Min | PA
---|---|---|---|---|---|---
J2000 | ($\mu$Jy/beam) | ($\mu$Jy) | (erg s-1) | (arcsec) | (arcsec) | (deg)
(1) | (2) | (3) | (4) | (5) | (6) | (7)
C-Band ($\nu$ = 6 GHz)
082341.08$+$0241805.7 | 335.8$\pm$6.4 | 365.0$\pm$14.0 | 42.131$\pm$0.001 | 0.333$\pm$0.105 | 0.069$\pm$0.034 | 106$\pm$11
082341.08$+$0241805.0 | 237.1$\pm$6.4 | 254.0$\pm$13.0 | 41.556$\pm$0.002 |  $<$0.969$\pm$0.043 |  $<$0.357$\pm$0.006 | 112$\pm$1
Ku-Band ($\nu$ = 15 GHz)
082341.08$+$0241805.7 | 61.5$\pm$6.9 | 60.0$\pm$12.0 | 41.745$\pm$0.009 |  $<$0.137$\pm$0.019 |  $<$0.080$\pm$0.007 | 160$\pm$7
082341.08$+$0241805.0 | 59.6$\pm$7.9 | 98.0$\pm$19.0 | 41.54$\pm$0.005 | 0.113$\pm$0.027 | 0.046$\pm$0.033 | 103$\pm$21

Note. — (1) J2000 coordinate of the identified radio source core, where we
list the N nucleus first; (2) Peak flux density at the central frequency of
the band within the source extent; (3) Integrated flux density at the central
frequency of the band within the source extent; (4) logarithm of the rest-
frame luminosity density at the central frequency of the band; (5$-$8) Best
fit beam-deconvolved source sizes (FWHMs in arcsec) along the major and minor
axes and the position angle of the major axis (degrees east of north). In
cases of unresolved point sources (denoted with $<$), the values for are for
the un-deconvolved regions. All errors quoted represent 1-$\sigma$
uncertainties which are derived via the correlated noise prescription of
Condon (1997). The error of the photometry does not include the 3% uncertainty
in the VLA flux density scale (Perley & Butler, 2013).

###  3.3 Radio Imaging

Figure 9: Top row: VLA Ku-band image decomposition of SDSS J0823$+$2418 at 15
GHz. We show from left to right: the continuum image, the model using imfit,
and the residual map. The restoring beam is shown in gray in the bottom right
corner of each panel. Contours begin at 3 times the rms background and
increase exponentially to the peak value. The legend in the right panel
indicates the N and S radio components, as well as the compact central blob
that might correspond to the foreground lensing galaxy. We overlay the $HST$
F160W best-fit coordinates, aligned by the N source, which show agreement to
within 0.1″. Major tick marks are spaced by 0.5″. Bottom row: same as the top
row, but for the VLA C-band continuum imaging at 6 GHz. The best-fit locations
between the two bands are consistent to within 0.1″. While the N source is
clearly brighter and more spatially extended at 6 GHz, the reverse situation
is seen in the 15 GHz image, indicating different spectral slopes between the
N and S quasars. We do not find evidence for a lensing galaxy component in the
6 GHz frame

We use the CASA task imfit to quantify the spatial extents and flux densities
of the sources, simultaneously fitting with a Gaussian component for each
source. The resulting best-fit quantities are given in Table 6. Figure 9 shows
the cleaned imaging, model, and residual maps for both the C and Ku-band data.
The spatial extents are slightly larger than the restoring beam sizes, and
show an interesting dichotomy: the N source appears more spatially extended at
6 GHz while the S source is more extended at 15 GHz. This could be suggestive
of two separate quasars. Specifically, we find that SDSS J0823$+$2418 N is
unresolved in the Ku-band image (beam FWHM: 0.11$\times$0.10″), while SDSS
J0823$+$2418 S is only marginally resolved in the C-band image (beam FWHM:
0.92$\times$0.35″). This is reflected in the markedly different integrated
flux values, and thus the spectral indices, between the two sources. SDSS
J0823$+$2418 N has $\alpha_{\rm 6-15\ GHz}$ = $-$1.971$\pm$0.222, while SDSS
J0823$+$2418 S has an index of $\alpha_{\rm 6-15\ GHz}$ = $-$1.039$\pm$0.219.
Both sources exhibit steep negative values consistent with non-thermal
optically-thin synchrotron emission. The observed-frame flux ratio between the
N and S source is 1.437$\pm$0.064 at 6 GHz, which is consistent with that
measured in the optical and NIR. However, due to the disparity of the spectral
indices, the flux ratio at 15 GHz is 0.612$\pm$0.280, which seems to argue
against a lensing scenario.

The Ku-band image also shows an intriguing central clump of residuals. At
first glance, its location appears to be roughly consistent with the putative
foreground lensing galaxy as seen in the Keck imaging. We add an elliptical
component during the model fitting to incorporate this blob, as shown in the
top middle panel of Figure 9. The additional component does not alter the
best-fit solutions for the N and S sources, and can be successfully modeled.
However, this unresolved central source is slightly smaller than the restoring
beam which suggests that it could be noise, and it has an integrated flux
density consistent with the background within its uncertainties yielding a
detection at only the 2-$\sigma$ level.

The faintness of the blob is not surprising giving the faintness of the lens
component in the NIR bands; the observed low level of radio emission is likely
not due to an AGN, although could be attributed to star formation processes. A
component between the N and S sources cannot be successfully modeled for the
C-band image; adding the additional component results spurious best-fit
positions and non-physical flux values. While it is possible that a faint
source does exist between the two sources, they are blended at the location
where the additional fainter source would be located.

###  3.4 X-ray Spectral Fitting

Given the spatial resolution of $Chandra$ ($\sim$0.5″native pixel size) and
the fact that the SDSS J0823$+$2418 system appears somewhat blended even in
the higher-resolution IR imaging, we begin our analysis by confirming that
there are indeed two X-ray sources detected. To disentangle the sources from
the detector effects, we utilize the Bayesian AnalYsis of Multiple AGN in
X-rays (baymax) algorithm, which is described in detail in Foord et al. (2019,
2020). Briefly, Baymax returns the likelihood of the two-source model over the
one-source model via the Bayes factor $\mathcal{B}$, and probabilistically
assigns each X-ray event to one of the sources in the two-source model.

In Figure 10, the counts are color-coded based on whether they are most likely
originating from either the N or S source. For both observations, we restrict
baymax to a 20″$\times$20″ region centered on SDSS J0823$+$2418\. While
visually obvious, we confirm that observation 1 is more consistent with a
single source model with $ln(\mathcal{B})=0.6\pm 1.4$, where all 33 counts are
associated with the N source. For observation 2, we obtain
$ln(\mathcal{B})=17.3\pm 1.8$, which is strongly in favor of two source model.
Here, 57 counts are associated with the N source, and 103 are associated with
the S source. The uncertainty on $\mathcal{B}$ is the 68% confidence interval,
taken from the statistical errors on $\mathcal{B}$, returned from the nestle
package.

Figure 10: Chandra X-ray observations and image decomposition. Top row:
Imaging of observations 1 and 2 (left and right panels) for SDSS J0823$+$2418
where the color bar indicates observed broad-band counts. In observation 1, no
source is detected for the S quasar. In observation 2, we mark the best-fit
locations of the two sources from baymax. The location of the brightest pixels
has shifted noticeably between the two observations. Bottom left: baymax
assessment of joint posterior distribution for separation and flux ratio of
two sources which is roughly consistent with our galfit modeling results.
Bottom right: Decomposition of the broad-band counts for observation 2,
assigning each count to the most likely source.

In Figure 10, we show the posterior distribution of the separation ($r$) and
count ratio ($f$) derived by baymax via MCMC sampling for observation 2. We
find a median separation of 0.81″${}^{+0.19}_{-0.22}$, which is consistent
within the quoted 3$\sigma$ uncertainties with the separation determined from
Gaia. However, the median count ratio between the weaker N source and the
dominant S source during observation 2 is 0.51$\pm 0.56$, which is quite
different from the flux ratios determined above.

Figure 11: Chandra X-ray spectral fitting for SDSS J0823$+$2418 using baymax.
The target is denoted as either North (N) or South (S) in the label, and
observation number is given in parentheses. The SW nucleus is not detected in
the first observation, and then dominates in the later observation. The
spectral shapes are fairly similar, although the N source appears to have a
less steep cutoff in the soft X-ray regime than the S source for both
observations.

The extracted spectrum for each source is composed of the counts which have
been spatially decomposed by baymax, and the spectral fits are evaluated using
Cash statistics. The background spectrum is similarly determined for the field
(gray points in Figure 10). The spectra are fit using the combined XSPEC model
phabs * zphabs * zpow, representing a redshifted powerlaw with photoelectric
absorption by local material. Milky Way obscuration for phabs is set to the
value of Milky Way neutral Hydrogen column density in the direction SDSS
J0823$+$2418 ($N_{\rm H}=3.77\times 10^{20}$ cm-2), found using the propcolden
tool (Dickey & Lockman, 1990). We fix the photon index to the canonical value
for AGN of $\Gamma=1.8$. We did attempt leaving $\Gamma$ as a free parameter,
but this consistently resulted in poorer fits. Sampling from the posterior
distributions, baymax creates 100 spectral realizations for both the N and S
source, and then fits them via XSPEC.

Table 7: X-ray properties of SDSS J0823$+$2418

Target | Counts | $F_{\rm 2-7\ keV}$ | HR | $L_{\rm 2-7\ keV}$
---|---|---|---|---
(1) | (2) | (3) | (4) | (5)
SDSS J0823$+$2418 N(1) | 33 | 12.0${}^{+0.2}_{-0.5}$ | 0.21${}^{+0.03}_{-0.12}$ | 13.0${}^{+0.6}_{-1.0}$
SDSS J0823$+$2418 N(2) | 57 | 4.5${}^{+2.4}_{-2.1}$ | 0.38${}^{+0.12}_{-0.24}$ | 4.6${}^{+4.1}_{-2.1}$
SDSS J0823$+$2418 S(2) | 102.5 | 8.6${}^{+2.0}_{-2.6}$ | 0.54${}^{+0.05}_{-0.1}$ | 12.0${}^{+3.0}_{-3.0}$

Note. — X-ray properties derived using baymax on the individual $Chandra$
ACIS-S observations. (1) Target and frame of observation (the S source is not
detected in observation 1); (2) broad-band net photon counts; (3) observed
flux in units of 10-14 erg s-1 cm-2; (4) hardness ratio: HR = (H-S)/(H+S),
where H and S are the net counts in the hard and soft X-ray bands,
respectively; (5) rest-frame unabsorbed luminosity from best-fit spectral
model (assuming fixed power law photon index $\Gamma=1.8$), in units of 1044
erg s-1. Values quoted are the medians resulting from 100 spectral
realizations, and all errors are quoted at the $\sim$3-$\sigma$ confidence
level.

We show all the spectral realizations for each source and observation in
Figure 11. We give the resulting derived quantities in Table 7. The average
intrinsic luminosities of $L_{\rm X}\sim 10^{45}$ erg s-1 are securely within
the quasar regime. The intrinsic unabsorbed flux ratio between the two sources
during observation 2 is $0.38^{+0.35}_{-0.20}$, which is starkly different
from the flux ratios found in the optical and UV regimes. The spectra for both
observations of the N source exhibit some degree of soft X-ray counts while
the S source shows a steeper decline below 2 keV. This is reflected in the
median values of the hardness ratio, HR $=\frac{H-S}{H+S}$, where and $S$ and
$H$ are the counts in the soft (0.5$-$2 keV) and hard (2$-$7 keV) bands,
respectively. The value of HR is consistent within the uncertainties between
the two observations for the N source (0.21${}^{+0.03}_{-0.12}$ vs.
0.38${}^{+0.12}_{-0.24}$), despite its flux diminishing by $\sim$50%. We find
a higher HR (0.54${}^{+0.05}_{-0.1}$) for the S source, albeit within the
scatter of the N source for observation 2. Were this a bona fide dual quasar,
the HR values could be interpreted as a higher degree of obscuration in the S
source. Differential extinction could also be due to different lines of sight
through the foreground galaxy. For a quadruply lensed quasar, Glikman et al.
(2023) found low degrees of differential extinction ($<$0.15 mag in
NIR/optical), which would also contribute to differential X-ray absorption.
Given that the two sources are unobscurred broad line quasars in optical/UV
and the previous evidence in support of the lensing scenario, this temporal
evolution of the two apparent sources instead hints that perhaps some
intrinsic variability of is at play, or a lensing-based time lag, or a
combination of both effects.

##  4 Discussion

We first combine the findings from the various lines of evidence in our multi-
wavelength analysis. Our primary goal is to assess whether SDSS J0823$+$2418
is a dual AGN or if it is a single AGN that has been gravitationally lensed by
a foreground galaxy. Our IR image decomposition has shown that inclusion of a
central galaxy component between the two sources, even when accounting for the
host galaxies, yields the lowest-residual model fits that are likewise
statistically preferred. In the Keck image, the location of this central
galaxy appears marginally discernible above the background; and while not
directly resolved in the HST image, the residual flux centrally concentrated
between the two quasars and hosts is significant at $>5\sigma$. There is also
tentative evidence of the same central galaxy in the Ku-band radio imaging,
although it is not resolvable in the other band. The strikingly similar
UV/optical spectra of the two quasars are naturally explained as lensed images
of the same source, such that the emission line and continuum properties are
in good agreement within their uncertainties. Finally, the optical and IR
photometry, as well as the continuum UV luminosities yield consistent flux
ratios between the N and S sources indicating that their underlying spectral
energy distributions (SEDs) are the same modulo a scaling factor that can be
attributed to magnification due to lensing.

We plot the SEDs of SDSS J0823+2418 based on our multi-wavelength photometry
in Figure 12. We overplot several template SEDs of typical quasars from the
literature for comparison (Richards et al., 2006b; Shang et al., 2011). The
SEDs of the N and S source are qualitatively similar across the rest-frame
optical/UV regime, and both show a steeper slope than the template SEDs
towards the Keck $K_{\rm p}$ filter. We show the differing values for the S
source between models 1 and 3 in the F160W filter, where the Model 1 value
produces a flux ratio that is consistent with the other UV/optical bands. Both
sources appear to be underluminous in X-ray and radio bands compared to the
optically bright quasar templates. The large uncertainty for the X-ray
luminosity of the N source during observation 2 renders the value roughly
consistent with the value from observation 1.

Figure 12: SEDs for SDSS J0823+2418. The top and bottom panels show the SEDs
for the N and S source, respectively. We also plot optically luminous (OL),
optically dim (OD), and the mean (All) SEDs of optically selected quasars from
SDSS as a comparison (Richards et al., 2006b). Typical SEDs of optically
bright non-blazar radio-loud (RL) and radio-quiet (RQ) quasars are also
plotted for reference (Shang et al., 2011). The literature SEDs are normalized
to the $HST$ F814W luminosities of the target. For the N source, we plot the
X-ray luminosity from observation 1 using a colorless star. The pink stars
refer to observation 2 luminosities in both panels. For the F160W fluxes, the
red squares refer to Model 3, while the colorless squares refer to Model 1 (N
value is nearly identical in both models).

Since both sources have radio luminosities near the radio-quiet quasar
template, we compute the radio loudness parameter, $R_{\rm 6\ cm/2500\
\mathring{A}}$, as the ratio of the flux densities at rest-frame 6 cm and 2500
Å. We scale the observed 6 GHz flux densities to 6 cm, assuming a power law
with the spectral index values $\alpha_{\rm 6-15\ GHz}$ observed for each
source. Radio-loud quasars are typically defined as having $R_{\rm 6\ cm/2500\
\mathring{A}}>10$ (Kellermann et al., 1989). We find that the N source is
above this threshold with $R=19.32\pm 5.27$, while the S source is radio-quiet
with $R=5.33\pm 1.57$. This difference is mostly due to the two-band spectral
indices. If we instead compute $R$ assuming the canonical value of
$\alpha\sim-0.7$ for AGN due to synchrotron emission, then we find that values
for the N and S sources ($R=4.12\pm 1.12$ and $R=3.53\pm 0.97$, respectively)
are consistent and definitively in the radio-quiet regime.

Despite the large amount of evidence pointing towards a lensed source, there
remain several pieces of evidence which appear to be more suggestive of two
distinctive quasars. While the C-band flux ratio is consistent with those at
other wavelengths, the higher-resolution Ku-band flux ratio is reversed such
that the S source is brighter at 15 GHz. The different inter-band spectral
index values between the two radio sources suggest unique AGN signatures.
There are also several differences in absorption features between the two STIS
UV/optical spectra, although some of these are artefacts due to uncleaned
cosmic rays and noise and do not appear in the GMOS spectra. Finally, in the
X-ray regime, the time-series of the two observations reveals the system to be
dynamic, where the N source heavily dominates at first and then the S source
flares, also leading to a flux ratio different from the other bands. SDSS
J0823$+$2418 falls into the newly-discovered category of so-called nearly
identical quasars ($e.g.$, Chan et al., 2022), where the properties of two
closely separated quasars suggest lensing, but a foreground lensing galaxy is
not apparent.

Some of the apparent differences between the two quasars might be due to
lensing-based flux anomalies caused by the foreground lensing galaxy. Xu et
al. (2012) find that dark matter halos and matter substructure ($e.g.$, from
microlensing by stars) in simulated lensing galaxies can introduce a 20-30%
flux ratio anomaly for a lensed quasar at $z=2$. Observationally, evidence for
the effects of foreground substructure has been been found using ALMA on the
lensed system SDP.81 (Hezaveh et al., 2016). Glikman et al. (2023) have found
flux anomalies between lensed images at different wavelengths. Based on
simulations, they were not able to reconcile the observed wavelength-dependent
flux anomalies with microlensing; however, their model of dark matter
substructure in the lensing galaxy close to one imaged quasar was able to
reproduce the high magnification needed to cause the flux anomaly in that
wavelength. Since the apparent foreground galaxy in SDSS J0823$+$2418 is only
marginally detected in NIR and 15 GHz radio and appears to be compact, it is
difficult to say whether foreground substructure plays a part in the observed
wavelength-dependent flux anomalies of the system.

Radio interferometers at high angular resolution suffer from spatial
brightness sensitivity filtering and thus introduce biases into component flux
density measurements as compared to single dish observations, depending on the
source substructure and scale of the interferometer. For poorly sampled $uv$
coverage, this can amount to a deficit in the observed flux densities. This
could be contributing to the flux anomaly observed specifically in the higher-
resolution Ku-band image.

Combining all of the evidence, we suggest that the most likely interpretation
of SDSS J0823$+$2418 is that it is a single quasar which is gravitationally
lensed into two images by a marginally detected foreground host galaxy,
although this conclusion comes with the caveats noted above. In the remainder
of the text, we probe the nature of the SDSS J0823$+$2418 system under this
framework.

###  4.1 Lensing Tests

Figure 13: Lens Mass modeling for SDSS J0823$+$2418\. The columns give the
decompositions for the HST F160W fitted results based on Model 1 (left) and
Model 3 (right). Top: image plane, where red triangles denote the observed
relative positions of the N and S sources, the central blue dot is the
location of the lensing galaxy as determined from our GALFIT analysis, and the
blue ellipse is the critical curve. Bottom: decomposed source plane using
glafic (Oguri, 2010), where the red triangle shows the position of the single
quasar relative to the cut and caustic. Despite the differing moth models can
be well-described by strong lensing using a simple SIE model.

We assess the feasibility of lensing to produce the observed source positions
and magnitudes by conducting lens mass modeling using the software glafic
(Oguri, 2010). We test the effects of the differing flux ratios produced by
the F160W Models 1 and 3. In Figure 13, we plot the de-projected source
positions and the corresponding observed image positions. In both cases, the
observed images are consistent with strong lensing. The positions of the image
components relative to the lensing galaxy naturally explain the enhanced
magnification observed for the N source in most of our observations. Given the
on-sky separation of the images and the spectroscopic redshift of the quasar,
we compute the minimum brightness needed for the foreground galaxy to produce
the strong lensing using the fundamental plane relation of Faber & Jackson
(1976), finding a minimum brightness of m$\sim$21.4 AB. Our best-fit values of
m${}_{\rm lens}\sim 19.5-20.0$ AB are brighter than this threshold and thus
are consistent with a foreground lensing galaxy at redshift of either $\sim
0.3$ or $\sim 1.6$.

Given the rather different positions of each quasar in the image plane
relative to the source quasar, we expect the different light travel paths to
result in a lensing-based time lag between the S and N images. This is
possibly corroborated by the X-ray observations which show variability over
the course of $\sim$3 days between exposures. We obtain an estimate for the
time lag using a relation for time delays in glafic based on a generalized
isothermal potential (Witt et al., 2000; Oguri, 2007; Lieu, 2008):

$\Delta t_{ij}=\frac{1+z_{\rm l}}{2c}\frac{D_{\rm ol}D_{\rm os}}{D_{\rm
ls}}(r^{2}_{j}-r^{2}_{i}),$ (3)

where the angular diameter distances $D$ are between the observer (o), lens
galaxy (l), and source (s), and $r_{i}$ is the distance between image
component $i$ and the center of the lensing galaxy. Assuming the lensing
galaxy redshift of $z\sim 1.6$, we find that the S source experiences a time
delay of $\Delta t\sim 112$ (36) days behind the N source based on Model 1
(Model 3). If the lensing galaxy is actually at $z\sim 0.3$, this time delay
is much shorter, $\Delta t\sim 2.8$ (0.9) days. Clumps of matter along the
line of sight and microlensing should not introduce any further time lags
(Lieu, 2008).

Figure 14: Archival SDSS spectrum of the SDSS J0823$+$2418 system, which has
more wavelength coverage towards the observed-frame blue than the STIS or GMOS
spectra. The presence of strong absorption lines around observed-frame 4500Å
and 4200Å consistent with Mg and Fe suggest that an intervening absorber is
located at $z\sim$0.62. These lines could be caused by relatively low column
density gas in the intergalactic medium, but could be consistent with the
foreground lensing galaxy.

The archival SDSS spectrum for the system, shown in Figure 14, appears to
exhibit absorption features that are not consistent with the redshift of the
quasar system. Instead, the complexes at observed-frame $\sim 4530$Å and $\sim
4200$Å seem to be consistent with cold ISM absorption lines of Mg ii 2796 and
2803$\lambda\lambda$ and Fe ii 2586 and 2600$\lambda\lambda$, respectively
(Bowen et al., 2000). These absorption features suggest that the redshift of
the absorber is $z\sim 0.62$, which in turn implies time delays of $\Delta
t\sim 7.5$ and 2.4 days for Models 1 and 3, respectively. Given the error on
the best-fit magnitude of the lensing galaxy as well as the scatter of the
Faber-Jackson relation, the lens galaxy redshift of 0.62 is not incompatible
with the predictions of 0.3 and 1.6 from our lens mass modeling. At $z=0.62$,
the lens galaxy would have an absolute magnitude of -22.88, similar to the
elliptical galaxy M87. The Faber-Jackson relation predicts a minimum absolute
AB magnitude of -21.73, which the observed lens galaxy magnitude satisfies.

The Mg ii and Fe ii absorption features might not originate in the foreground
lensing galaxy. There could be intervening intergalactic medium gas along the
line of sight, and these low-ionization metal absorption features require only
modest column densities to be produced ($10^{16}<$ N(Hi) $\leq 10^{22}$ cm-2,
Kacprzak & Churchill, 2011). However, these Mg ii absorption lines are known
as a common marker in the spectra of lensed quasar systems. In a study of
$\sim$266,000 SDSS quasars, Raghunathan et al. (2016) found $\sim$37,000
systems across a range of redshifts with evidence of foreground Mg ii
absorption. As a rough estimate of how common these absorption features are in
systems like SDSS J0823$+$2418, we visually inspected $\sim$16 SDSS spectra of
quasars at similar redshift and magnitude ($z\sim 1.8$, $r$ mag$\sim 17$), and
found 4 quasars with strong absorption lines similar to those in the SDSS
spectrum of SDSS J0823$+$2418\. We estimate the probability that a given
system exhibiting Mg ii absorption lines is a lensed quasar as:

$\displaystyle p({\rm lensing}\ |\ {\rm MgII})$ (4)
$\displaystyle=\frac{p({\rm MgII}\ |\ {\rm lensing})p({\rm lensing})}{p({\rm
MgII}\ |\ {\rm lensing})p({\rm lensing})+p({\rm MgII}\ |\ {\rm dual})p({\rm
dual})}$ (5) $\displaystyle\sim\frac{p({\rm lensing})}{p({\rm lensing})+p({\rm
MgII}\ |\ {\rm dual})p({\rm dual})}$ (6) $\displaystyle=\frac{1}{1+p({\rm
MgII}\ |\ {\rm dual})p({\rm dual})/p({\rm lensing})}$ (7)
$\displaystyle\sim\frac{1}{1+p({\rm MgII}\ |\ {\rm dual})}.$ (8)

The first equation (lines 5 and 6) follows Bayes’ theorem. The second equation
assumes that $p({\rm MgII}\ |\ {\rm lensing})=1$, $i.e.$ lensing systems
always show Mg ii absorption. The third equation divides both the numerator
and the denominator by $p({\rm lensing})$ so one can see that the ratio of
$p({\rm dual})/p({\rm lensing})$ also plays some role. The last equation
assumes that $p({\rm dual})/p({\rm lensing})=1$, based explicitly on our
results so far from two VODKA targets, but also informed by estimates from
Chen et al. (2022) and Shen et al. (2023) that the two quantities are likely
of the same order of magnitude, resulting in our final expression. We lastly
assume that $p({\rm MgII}\ |\ {\rm dual})=p({\rm MgII}\ |\ {\rm normal\
quasars})$ since foreground absorption occurs randomly and the population of
known dual quasars at similar redshift is too low for a statistical
assessment. Of course, there are some uncertainties in the assumptions of
$p({\rm MgII}\ |\ {\rm lensing})=1$ (which may depend on S/N ratio and
wavelength coverage) and $p({\rm dual})/p({\rm lensing})=1$, but it serves as
a useful expression for us to estimate such probability. Having found $p({\rm
MgII}\ |\ {\rm normal\ quasars})\sim 0.25$ above, we roughly estimate that Mg
ii absorption lines found in other VODKA targets have $\sim$80% chance of
indicating lensed quasars, making this a potentially powerful diagnostic test
for ruling out dual quasar candidates when optical spectra are available.

A time lag of roughly 2$-$3 days might explain the observed X-ray variability
of the system. Given the $\sim$3 days between the two $Chandra$ observations,
the N source is fairly luminous at first, perhaps exhibiting the onset of an
X-ray enhancement. By the time of the second observation several days later,
this signature has begun to fade by a factor of $\sim 2.8$ in the N image
location. Meanwhile, the S image has not yet seen the X-ray enhancement at the
time of observation 1 and is thus undetected, but subsequently exhibits a
brightening several days later to a luminosity comparable to the N image
during observation 1, but $\sim 2.6\times$ brighter than the N source during
observation 2. The fact that $L_{\rm X}$ is slightly lower for S(obs 2) than
N(obs 1) is likely a combination of the timing and the higher magnification
for the N image.

The X-ray variability observed for SDSS J0823$+$2418 is similar compared to
other several other AGN. Niu et al. (2023) have found X-ray flux variations in
M81* on the scale of 2 - 5 days, and these flares were on the order of a
factor of 2 brighter than the quiescent periods. M81* is smaller mass than the
black hole in SDSS J0823$+$2418 by several orders of magnitude, and the X-ray
variability of AGN is known to be anti-correlated with black hole mass (Ludlam
et al., 2015, and references therein), suggesting a longer expected timescale.
Huang et al. (2023) found a $\sim 5\times$ increase in the X-ray flux of
J081456.10+532533.5 between Chandra and XMM-Newton observations taken $\sim$30
days apart. During this time, they observed no noticeable variability in the
optical, concluding that X-ray variability in AGN typically occurs on shorter
timescales and at higher amplitudes than variability in the UV/optical range,
as evidenced by a handful of archival type 1 quasars with strong X-ray flux
variations and a lack of UV/optical variability (Huang et al., 2023, and
references therein). The STIS rest-frame UV spectra of SDSS J0823$+$2418 were
taken $\sim$8 days apart from the X-ray observations, so it is plausible that
the lack of UV deviations from the flux ratios observed in lower energy bands
is not tied to the observed level of X-ray variability.

Lastly, we investigate further whether the faint blob observed in the VLA Ku-
band can be explained by star formation radio emission from the lensing
galaxy. Assuming all of the observed emission at 15 GHz is due to star
formation, we extrapolate to the flux at 1.4 GHz assuming a power law with
spectral index $\alpha=-0.8$. We convert the flux to a $k$-corrected
luminosity density for each likely redshift (Zakamska et al., 2016), as
determined from the lens mass modeling and the SDSS spectrum. We then compute
the expected star formation rate (SFR) due to nonthermal radio emission using
the empirically-derived prescription of Bell (2003).

We obtain values of SFR = 15.3, 81.8, and 770 $M_{\odot}$ yr-1 at $z$ = 0.3,
0.62, and 1.6, respectively. Based on NIR samples from 3D-HST and the Cosmic
Near-IR Deep Extragalactic Legacy Survey, typical star-forming galaxies at
1$\lesssim z\lesssim$2.5 have peak SFR $\sim 10-15\ M_{\odot}$ yr-1, which
decreases to $\lesssim 2\ M_{\odot}$ yr-1 at $z=0$ (van Dokkum et al., 2013;
Patel et al., 2013). If the blob is indeed the lensing galaxy and its radio
emission is entirely due to star formation, then it is implausible that the
galaxy resides at $z=1.6$. While the SFR at $z=0.62$ is also high, the actual
rate could be lower if some of the observed radio emission is partially due to
a weak AGN in the lensing galaxy. The lensing galaxy redshift of $\lesssim$
0.62, and thus a time delay of $\sim$3 days, is thus not ruled out.

###  4.2 Comparison to Other Lensed Quasars at High Z

Until very recently, the search for gravitational lenses had resulted in a
only modest population of strongly lensed quasars. Some of the first well
constrained examples were characterized using similar approaches of HST image
decomposition, yielding 8 systems at $z>1.5$ which have prominent foreground
lensing galaxies (Lehár et al., 2000). Large systematic searches for lensed
quasars have been conducted in the optical, such as the SDSS Quasar Lens
Search (SQLS; Oguri et al., 2006; Inada et al., 2008). Their sample of 11
lensed quasars contains 6 system at $z>1.5$, however they are all at source
separations greater than 1″ by design; the 10 additional lensed quasars that
they studied are mostly at $z>2$, with several at smaller separations.
Interestingly, SQLS also identified 30 genuine quasar pairs (sep$\lesssim$1″)
of which 21 had been previously unrecognized. This unintended consequence
highlights the complimentary nature of dual quasar and lensed quasar searches.

More recently, the Gaia DR2 has been exploited to construct a catalog of 86
confirmed and 17 probable lensed quasars, some of which are at cosmic noon
(Lemon et al., 2018, 2022). However, the majority of these targets are also at
larger separations than SDSS J0823$+$2418 (0.78$-$6.23″). The larger
separations result from a component of the sample selection, where galaxy
catalogs are matched to Gaia detections within a broad matching radius and the
large PSF of the unWISE data. The authors find that, compared to mock
catalogs, 56% of the lensed quasar population at $z>1.5$ likely remains
undiscovered and constitutes systems with small separations (sep $<1.5$″)
where the foreground lensing galaxy is faint and difficult to resolve from the
quasar images (Lemon et al., 2022).

Paper III employed the varstrometry technique to identify promising candidates
of dual quasars at smaller separations than those noted above. A related
technique, the Gaia Multi Peak selection method (Mannucci et al., 2022, 2023),
has recently uncovered several high-$z$ dual quasar candidates. Of the four
candidates further investigated by Ciurlo et al. (2023), two systems had
strong evidence for dual quasars based on spatially resolved H$\alpha$
emission from Keck OSIRIS. The J1608+2716 system shows three resolved nuclei
at $z\sim 2.6$, with separations between components under 0.3″. While one of
the H$\alpha$ profiles of the components appears to be convincingly distinct,
the remaining two line profiles have consistent line centers. Ciurlo et al.
(2023) suggest that the differing equivalent widths of the matching components
could be due to microlensing effects, but despite the suggestive extended arc
morphology of the system, a putative lensing foreground galaxy is not detected
in their ground-based data. Similar spatially resolved spectroscopic follow-up
could be used to distinguish lensed/dual quasars in the remaining VODKA
sample.

Future campaigns using the Rubin Observatory Legacy Survey of Space and Time
(LSST) will provide complimentary searches for lensed quasars. A mock catalog
of simulated LSST images of lensed quasars by Yue et al. (2022) demonstrates
that the reasonable average seeing of the survey ($\sim 0.7\arcsec$) will more
than double the number of observable lensed systems with $z<7.5$ to roughly
2.4$\times 10^{3}$ (although Lemon et al. (2022) find that their mock catalog
underpredicts the number of currently observed lenses by 44%, so the final
number of LSST-observable lenses might be higher). While most of their catalog
contains quasars with image separations greater than that of SDSS
J0823$+$2418, smaller separation lenses can potentially be identified with
other wide-area survey programs, such as with the Euclid Telescope, which will
have a spatial resolution of 0.1 arcsec/pixel in optical wavelengths (0.3
arcsec/pixel for NIR imaging) across a 15,000 deg2 area of sky (Laureijs et
al., 2011). The Nancy Grace Roman Space Telescope will also have spatial
resolution comparable to HST (0.11 arcsec/pixel), but with $\sim 200\times$
larger field of view and greater sensitivity (down to $\sim$27 AB mag) in NIR
across a 2,000 deg2 survey designed partially to target gravitational lensing
(Spergel et al., 2015). Identified candidates should then be followed up with
multi-wavelength verification, similar to our approach above, to firmly
separate out dual quasars from lensed systems.

##  5 Summary and Conclusions

In this work, we have focused on one of the most promising dual quasar
candidates from the VODKA sample that is at a redshift near cosmic noon. To
uncover the true nature of SDSS J0823$+$2418, we rely on multi-wavelength
imaging and spectroscopy. Based on the similarities between the flux ratios of
the N and S quasars across several different bands, combined with the
similarity of their spectra and the tentative detection of a central
foreground galaxy, we conclude that the system is likely a single lensed
quasar. This result is bolstered by our subsequent lens modeling, which shows
that the best-fit positions and magnitudes from image decomposition are
consistent with strong lensing. This assessment comes with several caveats due
to conflicting results of our radio and X-ray imaging; however, these flux
anomalies can potentially be explained by contributions from systematic
effects, intrinsic variability, and lensing-based time delay of roughly 2$-$3
days. SDSS J0823$+$2418 is thus one of only $\sim$100 verified gravitationally
lensed quasars at high-$z$; however, owing to our varstrometry selection
technique, it is one of the smallest separation cases of these systems
currently known.

As a case study, our analysis here gives us a road map for the rest of the
VODKA candidates. When combined with the verification of SDSS J0749+2255 as a
bona fide dual quasar (Paper IV), we now have a series of tested metrics for
robustly discerning between the dual versus lensed quasar scenarios.
Specifically, we have found that several crucial pieces of multi-wavelength
evidence are necessary for a conservative follow-up program. Spectroscopy is
essential for probing the rest-frame optical or UV wavelengths where quasar
emission lines are easily diagnosed while concurrently allowing for robust
determinations of redshifts. This should be combined with imaging at optical
or NIR wavelengths that can capture both faint distinguishing features, such
as tidal tails, while also being able to spatially resolve potential
foreground lensing galaxies. Our complimentary pairing of Keck and HST imaging
in nearly identical bands has allowed us to capitalize on these two necessary,
but often mutually exclusive advantages. Radio and X-ray imaging (and the
corresponding general spectral features such as power law indices) can offer
additional signs that two AGN are distinct. Based on the two VODKA targets
studied in detail so far resulting in one dual quasar and one lensed, the
VODKA selection method therefore appears to be just as likely to identify
candidate lensed quasars at close separations as it is to find duals.
Therefore, the dual population overall might be overestimated. Further multi-
wavelength confirmation of candidates similar to this study is necessary to
shed light on the true fraction of lensed versus dual quasars.

We thank A. Kemball for helpful discussions on strong gravitational lensing,
as well as Junyao Li and Ming-Yang Zhuang for helpful discussions on PSF
modeling and image decomposition. We also 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. HCH appreciates the discussion with Ting-Wen Lan on the
foreground Mg ii absorption. This work is supported by NSF grant AST-2108162.
YS acknowledges partial support from NSF grant AST-2009947. Support for
Program 23700377 (PI: X. Liu) was provided by the National Aeronautics and
Space Administration through Chandra Award Number GO2-23099X issued by the
Chandra X-ray Center, which is operated by the Smithsonian Astrophysical
Observatory for and on behalf of the National Aeronautics Space Administration
under contract NAS8-03060. Support for Program number HST-GO-15900 (PI: H.
Hwang), HST-GO-16210, and HST-GO-16892 (PI: X. Liu) 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. This work was supported by JSPS KAKENHI Grant
Numbers JP22H01260, JP20H05856, JP20H00181. Based in part on observations
obtained at the international Gemini Observatory (Program ID 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. 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. This work is based in part on observations obtained with the 6.5 m
Magellan-Baade telescope located at Las Campanas Observatory, Chile, and with
the Apache Point Observatory 3.5 m telescope, which is owned and operated by
the Astrophysical Research Consortium. The National Radio Astronomy
Observatory is a facility of the National Science Foundation operated under
cooperative agreement by Associated Universities, Inc. The scientific results
reported in this article are based in part on observations made by the Chandra
X-ray Observatory. Based on observations made with the NASA/ESA Hubble Space
Telescope, obtained from the Data Archive at the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are
associated with program GO-16210;270 (PI: X. Liu), and program GO-16892;
PI:225 (X. Liu).

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, 061102
  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Anderson & Bedin (2010) Anderson, J., & Bedin, L. R. 2010, PASP, 122, 1035
  * Arenou et al. (2018) Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17
  * Barrows et al. (2021) Barrows, R. S., Comerford, J. M., Stern, D., & Assef, R. J. 2021, The Astrophysical Journal, 922, 179
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
  * Bell (2003) Bell, E. F. 2003, ApJ, 586, 794
  * Bowen et al. (2000) Bowen, D. V., Roth, K. C., Meyer, D. M., & Blades, J. C. 2000, ApJ, 536, 225
  * Capelo & Dotti (2016) Capelo, P. R., & Dotti, M. 2016, Monthly Notices of the Royal Astronomical Society, 465, 2643
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, The Astrophysical Journal, 345, 245
  * Chan et al. (2022) Chan, J. H. H., Lemon, C., Courbin, F., et al. 2022, A&A, 659, A140
  * Chen et al. (2022) Chen, Y.-C., Hwang, H.-C., Shen, Y., et al. 2022, The Astrophysical Journal, 925, 162
  * Chen et al. (2023) Chen, Y.-C., Liu, X., Foord, A., et al. 2023, Nature, 616, 45
  * Ciurlo et al. (2023) Ciurlo, A., Mannucci, F., Yeh, S., et al. 2023, arXiv e-prints, arXiv:2301.03091
  * Condon (1997) Condon, J. J. 1997, Publications of the Astronomical Society of the Pacific, 109, 166
  * De Rosa et al. (2019) De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New Astronomy Review, 86, 101525
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Dickey & Lockman (1990) Dickey, J. M., & Lockman, F. J. 1990, Annual Review of Astronomy and Astrophysics, 28, 215
  * Dotti et al. (2012) Dotti, M., Sesana, A., & Decarli, R. 2012, Advances in Astronomy, 2012, 940568
  * Faber & Jackson (1976) Faber, S. M., & Jackson, R. E. 1976, ApJ, 204, 668
  * Foord et al. (2020) Foord, A., Gültekin, K., Nevin, R., et al. 2020, The Astrophysical Journal, 892, 29
  * Foord et al. (2019) Foord, A., Gültekin, K., Reynolds, M. T., et al. 2019, The Astrophysical Journal, 877, 17
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Garmire et al. (2003) Garmire, G. P., Bautz, M. W., Ford, P. G., Nousek, J. A., & Ricker, George R., J. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4851, X-Ray and Gamma-Ray Telescopes and Instruments for Astronomy., ed. J. E. Truemper & H. D. Tananbaum, 28–44
  * Glikman et al. (2023) Glikman, E., Rusu, C. E., Chen, G. C. F., et al. 2023, ApJ, 943, 25
  * Goulding et al. (2019) Goulding, A. D., Pardo, K., Greene, J. E., et al. 2019, ApJ, 879, L21
  * Gross et al. (2023) Gross, A. C., Fu, H., Myers, A. D., et al. 2023, The Astrophysical Journal, 945, 73
  * Guo et al. (2018) Guo, H., Shen, Y., & Wang, S. 2018, PyQSOFit: Python code to fit the spectrum of quasars, Astrophysics Source Code Library, record ascl:1809.008
  * Häussler et al. (2007) Häussler, B., McIntosh, D. H., Barden, M., et al. 2007, ApJS, 172, 615
  * Hernquist & Mihos (1995) Hernquist, L., & Mihos, J. C. 1995, Astrophysical Journal, 448, 41
  * Hezaveh et al. (2016) Hezaveh, Y. D., Dalal, N., Marrone, D. P., et al. 2016, The Astrophysical Journal, 823, 37
  * Holgado & Ricker (2019) Holgado, A. M., & Ricker, P. M. 2019, The Astrophysical Journal, 882, 39
  * Huang et al. (2023) Huang, J., Luo, B., Brandt, W. N., et al. 2023, arXiv e-prints, arXiv:2304.07323
  * Hui et al. (2017) Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541
  * Hwang et al. (2020) Hwang, H.-C., Shen, Y., Zakamska, N., & Liu, X. 2020, The Astrophysical Journal, 888, 73
  * Inada et al. (2008) Inada, N., Oguri, M., Becker, R. H., et al. 2008, The Astronomical Journal, 135, 496
  * Kacprzak & Churchill (2011) Kacprzak, G. G., & Churchill, C. W. 2011, The Astrophysical Journal Letters, 743, L34
  * Kellermann et al. (1989) Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B., & Green, R. 1989, AJ, 98, 1195
  * Kochanek et al. (2006) Kochanek, C., Schneider, P., & Wambsganss, J. 2006, in Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro (Springer-Verlag: Berlin)
  * Kormendy & Richstone (1995) Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581
  * Lallo et al. (2005) Lallo, M., Makidon, R. B., Casertano, S., Gilliland, R., & Stys, J. 2005, HST Temporal Optical Behavior & Current Focus Status, Instrument Science Report TEL 2005-03, 20 pages
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Lehár et al. (2000) Lehár, J., Falco, E. E., Kochanek, C. S., et al. 2000, The Astrophysical Journal, 536, 584
  * Lemon et al. (2022) Lemon, C., Anguita, T., Auger-Williams, M. W., et al. 2022, MNRAS
  * Lemon et al. (2018) Lemon, C. A., Auger, M. W., & McMahon, R. G. 2018, Monthly Notices of the Royal Astronomical Society, 483, 4242
  * Li et al. (2004) Li, J., Kastner, J. H., Prigozhin, G. Y., et al. 2004, The Astrophysical Journal, 610, 1204
  * Lieu (2008) Lieu, R. 2008, The Astrophysical Journal, 674, 75
  * Ludlam et al. (2015) Ludlam, R. M., Cackett, E. M., Gültekin, K., et al. 2015, Monthly Notices of the Royal Astronomical Society, 447, 2112
  * Mannucci et al. (2022) Mannucci, F., Pancino, E., Belfiore, F., et al. 2022, Nature Astronomy, 6, 1185
  * Mannucci et al. (2023) Mannucci, F., Scialpi, M., Ciurlo, A., et al. 2023, arXiv e-prints, arXiv:2305.07396
  * McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184
  * Niu et al. (2023) Niu, S., Xie, F.-G., Wang, Q. D., et al. 2023, MNRAS
  * Oguri (2007) Oguri, M. 2007, ApJ, 660, 1
  * Oguri (2010) —. 2010, PASJ, 62, 1017
  * Oguri et al. (2006) Oguri, M., Inada, N., Pindor, B., et al. 2006, AJ, 132, 999
  * Patel et al. (2013) Patel, S. G., Fumagalli, M., Franx, M., et al. 2013, The Astrophysical Journal, 778, 115, doi: 10.1088/0004-637X/778/2/115
  * Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, Astronomical Journal, 139, 2097
  * Perley & Butler (2013) Perley, R. A., & Butler, B. J. 2013, Astrophysical Journal Letters, 204, 19
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Powell et al. (2023) Powell, D. M., Vegetti, S., McKean, J. P., et al. 2023, arXiv e-prints, arXiv:2302.10941
  * Prochaska et al. (2020a) Prochaska, J., Hennawi, J., Westfall, K., et al. 2020a, The Journal of Open Source Software, 5, 2308
  * Prochaska et al. (2020b) Prochaska, J. X., Hennawi, J., Cooke, R., et al. 2020b, pypeit/PypeIt: Release 1.0.0, Zenodo
  * Raghunathan et al. (2016) Raghunathan, S., Clowes, R. G., Campusano, L. E., et al. 2016, MNRAS, 463, 2640
  * Richards et al. (2006a) Richards, G. T., Strauss, M. A., Fan, X., et al. 2006a, AJ, 131, 2766
  * Richards et al. (2006b) Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006b, ApJS, 166, 470, doi: 10.1086/506525
  * Rosas-Guevara et al. (2019) Rosas-Guevara, Y. M., Bower, R. G., McAlpine, S., Bonoli, S., & Tissera, P. B. 2019, Monthly Notices of the Royal Astronomical Society, 483, 2712
  * Ruffio et al. (2019) Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200
  * Sahu (2021) Sahu, K. 2021, in WFC3 Data Handbook v. 5, Vol. 5 (Baltimore: STScI), 5
  * Satyapal et al. (2017) Satyapal, S., Secrest, N. J., Ricci, C., et al. 2017, Astrophysical Journal, 848, 126
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, The Astrophysical Journal, 500, 525
  * Shajib et al. (2022) Shajib, A. J., Vernardos, G., Collett, T. E., et al. 2022, arXiv e-prints, arXiv:2210.10790
  * Shang et al. (2011) Shang, Z., Brotherton, M. S., Wills, B. J., et al. 2011, The Astrophysical Journal Supplement Series, 196, 2
  * Shen et al. (2019) Shen, Y., Hwang, H.-C., Zakamska, N., & Liu, X. 2019, The Astrophysical Journal, 885, L4
  * Shen et al. (2021) Shen, Y., Chen, Y.-C., Hwang, H.-C., et al. 2021, Nature Astronomy, 5, 569
  * Shen et al. (2023) Shen, Y., Hwang, H.-C., Oguri, M., et al. 2023, ApJ, 943, 38
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757
  * van Dokkum et al. (2013) van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35
  * Vestergaard & Osmer (2009) Vestergaard, M., & Osmer, P. S. 2009, The Astrophysical Journal, 699, 800
  * Wassenhove et al. (2012) Wassenhove, S. V., Volonteri, M., Mayer, L., et al. 2012, The Astrophysical Journal, 748, L7
  * Witt et al. (2000) Witt, H. J., Mao, S., & Keeton, C. R. 2000, The Astrophysical Journal, 544, 98
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Xu et al. (2012) Xu, D. D., Mao, S., Cooper, A. P., et al. 2012, Monthly Notices of the Royal Astronomical Society, 421, 2553
  * Yu & Tremaine (2002) Yu, Q., & Tremaine, S. 2002, MNRAS, 335, 965
  * Yue et al. (2022) Yue, M., Fan, X., Yang, J., & Wang, F. 2022, The Astronomical Journal, 163, 139
  * Zakamska et al. (2016) Zakamska, N. L., Lampayan, K., Petric, A., et al. 2016, Monthly Notices of the Royal Astronomical Society, 455, 4191
  * Zhang et al. (2021) Zhang, T.-C., Guo, Q., Qu, Y., & Gao, L. 2021, Research in Astronomy and Astrophysics, 21, 212