VaDAR: Varstrometry for Dual AGN Using Radio Interferometry

Binary and dual active galactic nuclei (AGNs) are an important observational tool for studying the formation and dynamical evolution of galaxies and supermassive black holes. An entirely new method for identifying possible AGN pairs makes use of the exquisite positional accuracy of Gaia to detect astrometrically variable quasars, in tandem with the high spatial resolution of the Karl G. Jansky Very Large Array (VLA). We present a new pilot study of radio observations of 18 quasars (0.8 ≤ z ≤ 2.9), selected from the Sloan Digital Sky Survey DR16Q and matched with the Gaia DR3. All 18 targets are identified by their excess astrometric noise in Gaia. We targeted these 18 quasars with the VLA at 2–4 GHz (S band) and 8–12 GHz (X band), providing resolutions of 0.″65 and 0.″2, respectively, in order to constrain the origin of this variability. We combine these data with ancillary radio survey data and perform radio spectral modeling. The new observations are used to constrain the driver of the excess astrometric noise. We find that ∼44% of the target sample is likely to be either candidate dual AGN or gravitationally lensed quasars. Ultimately, we use this new strategy to help identify and understand this sample of astrometrically variable quasars, demonstrating the potential of this method for systematically identifying kiloparsec-scale dual quasars.


INTRODUCTION
In the canonical model of galaxy evolution, the formation of more massive galaxies proceeds via a series of major and minor mergers of their smaller counterparts, in addition to gas and dark matter accretion (e.g., Schweizer 1982Schweizer , 1996;;Toomre & Toomre 1972; Barnes & Hernquist 1992;Hibbard & van Gorkom 1996;Somerville & Davé 2015).Most massive galaxies contain a central supermassive black hole (SMBH; Kormendy & Richstone (1995)), generally of mass 10 6 − 10 10 M ⊙ .Thus, galaxy mergers should result in pairs of gravitationally-bound, synchronously feeding SMBHs that travel to the center of their host merger via dynamical friction, evolving from dual systems to binary systems that are expected to coalesce (Kormendy & Richstone 1995;Barnes & Hernquist 1992;Hopkins et al. 2008;Volonteri et al. 2016).Given the apparent scaling relations between the SMBHs and their host galaxy bulge properties, a potential co-evolution exists between an SMBH and its host galaxy (Ferrarese & Merritt 2000;Gebhardt et al. 2000;Heckman & Best 2014).*

NASA Postdoctoral Program Fellow
The timescale of a galaxy merger is on the order of hundreds of millions to a few billion years (Tremmel et al. 2018;Callegari et al. 2009).Dual SMBHs are precursor systems, defined as having two SMBHs with a separation beyond their mutual gravitational spheres of influence; i.e. separations of ≤ 100 kpc (Ellison et al. 2011;Liu et al. 2011;Burke-Spolaor et al. 2018).Binary systems are a more advanced stage of evolution, having separations within their mutual spheres of influence; i.e. separations ≤ 100 pc (Ellison et al. 2011;Liu et al. 2011;Rodriguez et al. 2006).
Characterization of these systems at all stages of evolution is critical for understanding SMBH evolution and co-evolution with host galaxies.This is made difficult by the lack of robust detections of dual and binary SMBH systems, particularly in the binary SMBH regime.During the course of the merger, gravitational torques drive gas towards the nuclei (e.g.Barnes & Hernquist 1992, 1996;Mihos & Hernquist 1996;Burke-Spolaor et al. 2018;Khan et al. 2020), which infalls and causes gas accretion onto the SMBHs.The SMBHs then ignite as active galactic nuclei (AGN; Hopkins et al. 2008;Hopkins & Quataert 2010;Blecha et al. 2018).
Of particular interest are candidates and confirmed systems with redshifts at which both the number density of luminous quasars and the global star formation rate density peak, often referred to as "cosmic noon," (1 ≤ z ≤ 3, Richards et al. 2006;Madau & Dickinson 2014).However, the dynamical evolution timescales of dual and binary AGN systems are comparatively short with respect to other evolutionary stages (Tremmel et al. 2018), and thus there are a substantially lower numbers of these systems predicted on the scale of tens of kpc, as they quickly (< 1 Gyr) advance to the next stage of evolution (Tremmel et al. 2018;Merritt 2013;Chen et al. 2020;Yu 2002).This issue is further compounded by the sample pre-selection criteria and the limited resolutions and sensitivities of current telescopes, leaving an observational gap between r p < 6 kpc at z > 1 (see Figure 1 in Chen et al. 2022a).At small separations (i.e.sub-arcsecond), the resolutions of radio interferometry, space-based optical imaging, or ground-based imaging with large apertures, adaptive optics, or optical/infrared interferometry are required.
The focus of this paper is to describe a method for using radio interferometric observations to select dual AGN candidates.The National Science Foundation's Karl G. Jansky Very Large Array (VLA) can reach the high angular resolution necessary to detect dual and binary AGN systems with separations on sub-arcsecond scales.Given their rarity, a pre-selection strategy to select promising targets for radio follow-up is crucial.
The paper is organized as follows.In Section 2, we introduce the varstrometry method.In Section 3, we define our target selection process.In Section 4, we describe the observations, data reduction, and analysis.We present our sample properties in Section 5 and radio spectral modeling in Section 6.Our findings are discussed in Section 7. Throughout this paper, all physical separations are the projected separation (r p ).A flat ΛCDM cosmology is adopted, with a Ω Λ = 0.69, Ω m = 0.31, and H 0 = 67.7 km s −1 Mpc −1 (Planck Collaboration et al. 2020).

VARSTROMETRY
A novel astrometric technique, previously used to detect unresolved stellar binaries (e.g., Makarov & Goldin 2016) through photometric variability-induced photocenter pseudo-motion, was applied to search for unresolved dual AGNs by Hwang et al. (2020).This technique, dubbed "varstrometry" in the context of dual AGN searches by Hwang et al. (2020), leverages the unprecedented astrometric precision of Gaia (Gaia Collaboration et al. 2016), which has mapped the positions, parallaxes, and proper motions of billions of stars in the Milky Way.It is sensitive enough to provide unprecedented sensitivity to the positions of hundreds of thousands of distant quasars (Brown et al. 2021).While this sensitivity has confirmed the presence of wavelength-dependent positional offsets in jetted AGNs (e.g.Makarov et al. 2017Makarov et al. , 2019)), it has also revealed the presence of astrometrically variable AGN (e.g.Shen 2021; Chen et al. 2022a).Because Gaia is progressively source-confused for separations less than ∼ 2 ′′ (2 ′′ corresponds to > 12.3 kiloparsecs for z > 0.5; Fabricius et al. 2021), it is not capable of discerning a close secondary AGN or other extended phenomena associated with the immediate AGN environment, such as jet production.However, the sub-milliarcsecond astrometric precision of Gaia allows for the motions of unresolved quasars to be detected.This makes astrometric variability a new discovery space for multi-AGNs and quasars.
The orbital periods of binary and dual AGNs are respectively hundreds to millions of years (e.g., Figure 2 in Dorland et al. 2020), and thus their positions are essentially fixed on the sky.This precludes a direct motion measurement by an astrometric mission.However, AGN exhibit intrinsic, stochastic variability measurable on timescales as short as hours (Sesar et al. 2011;MacLeod et al. 2012).In a dual AGN with a separation less than the effective angular resolution of Gaia, two components with varying brightnesses will appear to Gaia to have a shifting centroid (Hwang et al. 2020).
Gaia's resolution limits are such that, for an AGN pair system, individual light curves cannot be observed.Instead, a joint variability light curve encodes the stochastic variability of the system (see Figure 1).In cases where the apparent photocenter of the AGN is different from that of its host galaxy (e.g., a disturbed/merger system or a system in which the AGN obscuration level changes rapidly), the joint variability lightcurve of that system will be indistinguishable from that of a dual AGN system.Additionally, the varstrometry technique is not only sensitive to AGN pair systems but can also select samples contaminated with star+quasar superposition systems (Pfeifle et al. 2023), lensed quasar systems (e.g.Mannucci et al. 2022;Ciurlo et al. 2023;Inada et al. 2012Inada et al. , 2014;;O'Dowd et al. 2015), and any other morphol- ogy that might also drive the excess astrometric noise seen in AGN pair systems (e.g.single AGN in a host galaxy with bright stellar features).
In the case of AGN pair systems, as the mutual photocenter of the AGN pair wanders between the individual AGNs (again depending on their variability), Gaia measures a positional "jitter", which is representative of the astrometric variations.This measurement also provides a lower limit on the possible physical separation for the AGN pair (σ astro ; Hwang et al. 2020).For a pair of quasars with similar mean flux densities and rms fluctuations, the expected variability amplitude σ astro is where D is the separation of the component sources, and f and ⟨∆f 2 ⟩ are the total mean and rms fluxes of the unresolved system.Assuming a typical fractional rms of ∼10% (e.g.MacLeod et al. 2012), one can calculate the expected astrometric "jitter" to be ∼ 10 milliarcseconds for every 0.2 ′′ of angular separation.Note that this measurement is a lower limit on the angular separation.
High spatial resolution follow-up has been successful for varstrometry-selected samples.Chen et al. (2022a) followed-up a sample of 84 Gaia-identified AGN pair candidates with the Hubble Space Telescope (HST) and Gemini GMOS optical spectroscopy.Their search revealed two previously unknown dual AGN candidates (Chen et al. 2022a;Shen 2021;Chen et al. 2022b).They conclude that ∼ 40% of their HST resolved pairs are likely to be physical quasar pairs or gravitationallylensed quasars (Chen et al. 2022a).
In this work, we present the first results from a pilot VLA program of varstrometry-selected sources (Hwang et al. 2020).Each source has excessive astrometric variability as measured by Gaia.Our primary goal is to place constraints on the morphology of the targets and fully characterize the sample using high-resolution radio interferometric follow-up observations.In particular, we aim to identify any potential AGN pair candidate systems, to robustly constrain the morphologies, luminosities, and radio spectra of the sample in order to understand the sample that the radio+varstrometry method selects for, constrain the driver of the excess astrometric variability, and to gauge the efficiency of the radio+varstrometry method as a pre-selection for more systematic searches for AGN pairs.

TARGET SELECTION
The pilot sample of quasars for this study was selected from the SDSS quasar catalog for DR16 (DR16Q; Lyke et al. 2020), which contains the most complete sample of bona fide spectroscopically-confirmed quasars from SDSS/BOSS.DR16Q was cross-matched to Gaia's Early Data Release 3 (EDR3; Gaia Collaboration et al. 2016) to within 1.5 ′′ , ensuring spectroscopic fiber coverage of the Gaia counterpart (Fabricius et al. 2021).
There are several Gaia parameters that can act as indicators of positional noise in a target, as well as more that can be manipulated to search for AGN pair systems.For this sample, the important Gaia parameter (see Table 1) is astrometric excess noise (AEN; the previously mentioned astrometric "jitter"), and it is defined as the amount of statistical dispersion required such that Gaia's astrometric solution for the source leaves no unexplained variance (Lindegren et al.

2021).
Gaia also calculates the statistical significance of this value: astrometric excess noise significance (AENS).Limiting to sources with AENS > 5 ensured a highly statistically significant measurement of excess astrometric noise in the Gaia astrometric solution for each target.As Gaia's astrometry is inherently optimized for compact/unresolved sources (Makarov et al. 2012;Makarov & Secrest 2022), quasars with large, extended host galaxies will exhibit spurious astrometric excess noise, as Gaia's scanning window (∼18 pixels along-scan by ∼12 pixels across scan for the main astrometric field chips; Gaia Collaboration et al. ( 2016)) is filled by the host galaxy (Makarov & Secrest 2022;Souchay et al. 2022).This can be mitigated by making a cut on redshift, which we do here, of z > 0.5 (where 1 ′′ is ∼ 6.1 kiloparsecs).We also make a cut on the Gaia G magnitude of G < 20, as there is an over-abundance of sources near Gaia's sensitivity limit with high AEN, likely due to the current limitations of Gaia's astrometric pipeline (Hwang et al. 2020).After applying these limits on the sample, 148 systems remain.
It is generally accepted that between 1-10% of AGN are typically radio bright (Osterbrock 1993(Osterbrock , 1989)), so it was necessary to cross-match the 148 quasars remaining to various radio surveys before committing to pointed radio observations.The larger sample was matched to the VLA Sky Survey (VLASS; Lacy et al. 2019) in order to confirm the existence of a radio signature.The final sample was composed of 18 targets (see Table 2) exhibiting a VLASS detection.High sensitivity radio observations of all 18 targets were made with the VLA at S-band (2-4 GHz, central frequency 3 GHz), and X-band (8-12 GHz, central frequency 10 GHz), with the telescope in A configuration.This provides a resolution of 0.65 ′′ at S-band and 0.2 ′′ at X-band (Project Code 22A-384, PI Schwartzman).The observations reach a sensitivity of 18.1 µJy/beam at Sband, with a maximum recoverable scale of 18 arcseconds, and a sensitivity of 8.9 µJy/beam at X-band, with a maximum recoverable scale of 5 arcseconds.Table 3 lists the observational details for each target, including the primary flux density and complex gain calibrators, and synthesized beams.X-band targets were observed in repeated phase calibrator-target cycles, with a total on-source integration time of about 350 seconds per target, while S-band targets were observed in repeated phase calibrator-target-phase calibrator cycles, with a total on-source integration time of about 450 seconds per target.

Data Reduction
At S-band, the data were recorded with 16 spectral windows, each having 64 channels of 2 MHz width, and covering the 2-4 GHz band.For the X-band, the data were recorded with 32 spectral windows, each having 64 channels of 2 MHz width, and covering the 8-12 GHz band.The data were reduced and calibrated with the National Radio Astronomy Observatory's (NRAO) Common Astronomy Software Applications (CASA; CASA Team et al. 2022) VLA Pipeline 1.4.2, using CASA version 5.3.0.
The initial pipeline steps followed standard procedures; the data were Hanning smoothed, antenna position corrections were applied, in addition to ionospheric TEC corrections and requantizer gains.Calibration was performed similarly, using antenna delay, bandpass, and complex gain solutions.The flux densities for the primary calibrations were taken from the Perley & Butler (2017) extension to the Baars et al. (1977) scale.The gain solutions were then transferred to the target sources.
Once pipeline calibration was completed, radio frequency interference (RFI) was removed where necessary using a combination of the CASA tasks RFlag and TFCrop.To refine the gain calibration, several rounds of self-calibration were performed, beginning with phaseonly, and proceeding to amplitude and phase (see Table 2).In certain targets, self-calibration was limited (or not possible) due to low signal-to-noise.All imaging was completed with CASA, and was performed with a multi-term multi-frequency synthesis (MTMFS) deconvolver (Rau & Cornwell 2011) and Briggs weighting, with a robust parameter of 0.5.Two Taylor terms were used to model the frequency dependence of the sky, and to account for the VLA's wideband receivers.
When necessary, W-projection with 64 w-planes (Cornwell et al. 2005;Cornwell 2008) was employed to account for the wide-field errors.Clean masks were employed at all stages.

Flux Density Measurements
Flux densities were measured with the Python Blob Detection and Source Finder package (PyBDSF; Mohan & Rafferty 2015).Every source in an image was fitted with a Gaussian, and the flux density properties were recorded.All unresolved, point-like sources were fitted with a single Gaussian.For extended sources, a multi-Gaussian source was fitted.In some cases, for a cluster of closely associated peaks within an island (sources with multiple components), a single-Gaussian source was fitted to each peak.
Flux density scale errors were also taken into account.Following the VLA Observing Guide, a flux density scale calibration accuracy of 5% was assumed for both bands.Flux density calibrators 3C48 and 3C138 are currently undergoing flares.While the 3C138 flare does not appear to impact measurements at the observed bands, the single target observed with 3C48 as a flux density calibrator was assumed to have a flux density scale calibration accuracy of 10%.

SAMPLE PROPERTIES
We present new VLA observations of 18 quasars with significantly high astrometric excess noise.Targets were initially categorized by radio morphology, though we stress that morphology is not necessarily equivalent to the driver of the astrometric variability.An unresolved source features a single source that remains unresolved at sub-arcsecond scales.A resolved source features some structure at sub-arcsecond scales, and those targets have been further categorized as multi-component or extended targets.A more detailed discussion of each category follows, as well as a more in-depth look at each individual target.

Unresolved Targets
Nine of the 18 targets display an unresolved pointlike morphology at scales probed with the new VLA observations.Four of these are further identified as star+quasar superposition via an analysis of the SDSS spectra (see Section 7.1).Table 4 presents the properties of these nine targets, including flux densities, luminosities, and source angular and projected physical size, which is equivalent to an upper limit on the separation of any possible components.The new VLA images are presented in Appendix A, Figure 8.
While unresolved at sub-arcsecond scales, these sources may each be hiding more detailed structure at smaller scales.Upper and lower limits on the separation of any possible components were calculated from the new VLA observations and the varstrometry method respectively.As these sources are unresolved at subarcsecond scales, we may calculate an upper limit on their separation as the full extent of the unresolved source.As described in Section 2, AEN can be used, in combination with other quasar properties, to calculate a lower limit on the separation of any possible components.
While we may place constraints on these targets, without smaller-scale very long baseline interferometry (VLBI) observations, the astrometric driver of each system will likely remain uncharacterized.Some unresolved targets will likely be single AGN with a slightly increased intrinsic variability driving the AEN.Others, with follow-up observations, may be identified as quasar pair systems with separations that are not resolvable at these current scales.It is also possible for these sources to exhibit a variety of sub-milliarcsecond scale structure (i.e.jets or other extended emission associated with a single AGN) that is similarly unresolved.All of these possibilities could give rise to astrometric variability.Higher angular resolution follow-up will be required in order to confirm the morphology of these sources.

Resolved Targets
Nine of the 18 original targets display a resolved morphology.This is defined as any source with two or more clear radio peaks at either 3 GHz, 10 GHz, or both, or otherwise extended or complex emission.In order to help sort out targets which display a morphology consistent with that of an AGN pair, this category can be further subdivided into multi-component and extended/complex targets, both of which are discussed below.

Multi-component Targets
Six of the nine resolved targets have been identified as multi-component targets, defined as any source with two or more clear radio peaks at either 3 GHz, 10 GHz, or both.Table 5 presents the properties of these targets, including flux densities, luminosities, and angular and projected physical separations of the components at 10 GHz (with the exception of J1128+2402 (see Appendix A, Figure 9e) for which the separation was measured at 3 GHz due to the lack of a secondary radio peak in the 10 GHz observations).The new VLA images are presented in Appendix A, Figure 9.
Fluxes of the multi-component targets were measured as integrated flux densities.Given the lower resolution of the 3 GHz observations, some targets with resolved components at 10 GHz remain unresolved at 3 GHz.In those cases, the total flux density of the 3 GHz signa-  ture has been reported, while the flux densities of each component at 10 GHz have been separately reported.
Our multi-component targets range in angular separation from 0.33 − 1.11 arcseconds, and in physical separation from 2.86 − 9.50 kiloparsecs, at redshifts of 1.247 to 2.166.While further analysis is still necessary, this range of separations at these redshifts is a parameter space that represents an observational gap in the known population of AGN pair systems (Chen et al. 2022a;Pfeifle 2023), and thus is of great interest.
Given the breadth of the definition used for multicomponent targets, it is expected that the majority of the sample will not be genuine dual quasar candidates.For example, J0951+2635 has been identified via spectroscopy as a gravitaionally-lensed quasar (Schechter et al. 1998).However, with the currently available data, we calculate that 44% of the entire sample are targets identified as either lensed quasar systems or bonafide AGN pair systems.While the new VLA observations have directly observed the systems identified here as likely candidate dual quasar systems, follow-up observations (i.e.high angular resolution optical imaging and spectroscopy with HST, as targets lack any nearby reference stars necessary for adaptive optics) will be necessary to confirm the nature of these systems.

Extended and Complex Targets
Three of the 18 targets display jets or other extended emission more commonly associated with a single AGN, though multi-AGN systems are capable of exhibiting jet activity (e.g.3C75; Owen et al. 1985).This morphology generally provides an excellent indicator of the driver of the astrometric variability (Hwang et al. 2020).Table 6 presents the sample properties for the extended and jetted targets, including flux densities, luminosities, and the measured size of any extended emission on the sky.Final flux densities were measured as discussed in Section 4.3.J1215+4529 features a pair of canonical radio jets (see Appendix A, Figure 10a).Table 6 presents the integrated flux densities of the core and both lobes, in addition to angular and projected physical sizes of the jet extensions and of the core.J1415+1129 and J1433+4842 both exhibit other forms of extended emission associated with the quasar.Given the complex nature of the targets, we have assigned four components, each coincident with a peak in the 10 GHz radio signature (see Appendix A, Figures 10b and 10c), and measured the integrated flux density for each.At 3 GHz, the integrated flux density of the entire source was measured.Finally, the full extent on the sky of each target was measured (see Table 6).

RADIO SPECTRAL MODELING
We performed radio spectral modeling using data from a variety of public radio surveys (see Section 6.1).Radio spectral shapes, in combination with morphological information, reveal important information about the quasar environment, the physical conditions of the jets/lobes (e.g.absorbed vs. optically-thin), and the evolutionary stage of the radio source (Nyland et al. 2020;Bicknell et al. 1998;O'Dea 1998;O'Dea & Saikia 2021;Orienti & Dallacasa 2014).
The observed radio spectrum is a combination of absorption and emission processes.To quantify the distribution and the location of the spectral peak, least squares fitting was performed for the new VLA observations presented above, in addition to all available archival and commensal observations.Two basic synchrotron emission models were employed which will allow us to constrain important physical properties of each radio source.
• Standard Power Law: A standard power law model that is defined by S ν = S o ν α , where S ν is the flux density at frequency ν, α is the spectral index, and S o is the flux density at 1 GHz.The value and sign of the spectral index can be used to help distinguish between synchrotron and thermal emission (e.g., Patil et al. 2021Patil et al. , 2022)).
• Curved Power Law: A curved power law generated from the flux density at frequency ν, S ν , the spectral index α, and the flux density at 1 GHz S o , defined as S ν = S o ν α e q (lnν) 2 , where q represents the width of the peak, characterizes the degree of curvature, and is defined as ν peak = e −α/2q , where ν peak is the peak frequency.Significant spectral curvature is typically defined as |q| ≥ 0.2 (Duffy & Blundell 2012), and is indicative of absorption of low frequency emission.Absorption may be due to a high synchrotron optical depth within the source (i.e., synchrotron self-absorption (SSA); O'Dea & Saikia 2021) or free-free absorption from surrounding ionized gas (Bicknell et al. 1998).Understanding which process is in play can further constrain the radio source and environment properties.VLA 3 GHz (S-band) peak flux density ± error.Column 3: VLA 3 GHz (S-band) peak luminosity.Column 4: VLA 10 GHz (X-band) peak flux density ± error.Column 5: VLA 10 GHz (X-band) peak luminosity.Column 6: AGN pair separation where applicable, reported in arcseconds and kiloparsecs, measured at 10 GHz, with the exception of J1128+2402, for which the separation was measured at 3 GHz due to lack of a secondary peak in the 10 GHz observations.Column 7: Classification and component location.
Following the approach described in Patil et al. (2022), each source was fitted with a standard power law and a curved power law, making use of the Radio Spectral Fitting 1 tools (Patil et al. 2021).The new VLA data at 10 GHz were re-imaged and uv-tapered to produce similar resolutions in both the X-and S-bands.The uvtapered image properties are described in Table 7. Least squares fitting was performed with the values weighted by their measurement uncertainties.The results of the fitting are summarized in Table 8, and the radio spectra are presented in Appendix D, Figures 13, 14, and 15.Radio spectral classifications for the overall sample are discussed in Section 6.2.Important results from the 1 https://github.com/paloween/RadioSpectral Fitting radio spectral modeling are included in Section 7.3, as are spatially resolved spectral indices, where possible.

Archival Radio Surveys
Radio spectral modeling with broad frequency coverage required adding survey data to our new observations.Radio flux densities and their uncertainties were taken from the following archival surveys.The flux density values are taken directly from the published catalogs with corrections applied for calibration uncertainties.

LOTSS
The Low Frequency Array Two-meter Sky Survey (LOTSS; Shimwell et al. 2019Shimwell et al. , 2022) ) is an ongoing allsky survey at 144 MHz covering the northern sky above 34 degrees.The second data release covers about 27% of the sky with a point-source sensitivity of 83 µJy per beam.A calibration uncertainty of 10% was applied.

TGSS ADR
The Tata Institute of Fundamental Research (TIFR) Giant Metrewave Radio Telescope (GMRT) Sky Survey Alternative Data Release (TGSS ADR; Intema et al. 2017) is a 150 MHz survey covering about 90% of the sky, or about 36900 deg 2 .TGSS reaches a resolution of 25 ′′ , and has a 1σ sensitivity of between 3 and 5 mJy per beam.A calibration uncertainty of 10% was applied.

VLITE
The VLA Low-band Ionosphere and Transient Experiment (VLITE; Clarke et al. 2016;Polisensky et al. 2016) is a commensal system on the VLA operating at all configurations that records of 6000 hours of commensal data per year at 350 MHz.With the VLA in A configuration (as it was for all targets in this project), VLITE reaches a resolution of about 5.6 ′′ .A calibration uncertainty of 15% was applied.GHz survey covering 10,575 deg2 of the sky.The observations reach a resolution of about 5 ′′ and a sensitivity of ∼0.13 mJy per beam.A calibration uncertainty of 5% was applied.

VLASS
The VLA Sky Survey (VLASS; Lacy et al. 2019) is an ongoing 3 GHz continuum survey covering 33,885 deg 2 .Of the 3 planned epochs, the first two have been completed, and the third is ongoing.A source catalog 2 is now available (Gordon et al. 2021).The typical 1σ sensitivity of a single epoch image is ∼ 120 µJy per beam with an angular resolution of 2.5 ′′ .A calibration uncertainty of 3% was applied.

VIPS
The VLBA Imaging and Polarimetry Survey (VIPS; Helmboldt et al. 2007) observed 1128 flat-spectrum sources brighter than 85 mJy in the Northern Cap region of SDSS at 5 GHz, covering about 14% of the sky.The typical 1σ sensitivity of a VIPS image is ∼ 0.2 mJy per beam with an angular resolution of 1.4 milliarcseconds.A calibration uncertainty of 5% was applied.

Radio Spectral Shape Classification
All radio spectra were inspected in order to best categorize their spectral classifications.The shapes of both the standard power law and the curved power law were considered, as were the reduced χ 2 values for each fit.We note that, in the case of SDSS J104406.33+295900.9 (see Appendix D, Figure 13), neither the standard power law nor the curved power law accurately model the target's spectrum, possibly due to flux variability on years-long timescales3 .Additionally, we note that for the systems identified as star+quasar superposition (see Section 7.1), there may be a stellar contribution to the radio spectrum (see Figure 1 in Güdel 2002).In particular, we note the inverted spectrum in SDSS J172308.14+524455.6,indicating potential stellar contributions at higher frequencies (Dulk 1985;Güdel 2002).Quantifying the contribution from the thermal regime is not possible without further high frequency observations.For each target, we adopt the approach described in Patil et al. (2022), and divide classifications into the following categories: • Standard Power Law: This is the standard power law defined in Section 6.The power law spans the full range of frequencies presented in our radio spectra, and the reduced χ 2 is lower for the power law model, in comparison to the curved power law model.In this case, a spectral index of < -0.7 is consistent with the optically thin synchrotron emission expected of most AGN.
• Curved Power Law: This is the curved power law defined in Section 6.The curved power law spans the full range of the frequencies presented in our radio spectra, but does not display a peak or turnover.The reduced χ 2 is lower for the curved power law model, in comparison to the standard power law model, indicating significant deviation from the standard power law.
• Peaked Spectrum: A visual inspection of the radio spectrum reveals a peak or turnover within the range of frequencies.
• Flat Spectrum: The spectral index as measured from the best-fit standard power law or curved power law is |α| < 0.5.
• Upturned Spectrum: A visual inspection of the radio spectrum reveals a upturned concave spectrum.
• Inverted Spectrum: The spectral index as measured from the best-fit standard power law or curved power law is α > 0.5.
In addition to the above classifications (see Table 8), the radio spectra were also characterized using a radio color-color plot (see Figure 3).Three flux densities measured at three sufficiently separated frequencies were used to calculate two spectral indices: the new 3 GHz and 10 GHz VLA observations, and either the FIRST (1.4 GHz) or the NVSS (1.4 GHz) survey values.These flux density measurements produce spectral indices α 3 1.4 and α 10 3 .The quadrants of the radio color-color plot represent classifications.
Figure 3 illustrates the distribution of our target spectral shapes and reveals that the sample is diverse in radio spectral shape.However, the majority of our targets, independent of the potential driver of excess astrometric noise, exhibit a spectral index that is consistent with, slightly flatter, or slightly steeper than that of optically thin synchrotron emission (lower left quadrant).We note that the impact of having sources with time variable behavior or extended emission that is not identically sampled at all frequencies (in this case, the 1.4 GHz survey fluxes) can result in apparent shifts in the source location in this color-color plot.In particular, sources exhibiting complex morphologies may have extended emission at 1.4 GHz that is resolved out in the higher frequency measurements, artificially shifting their location on the color-color plot to the left.We refer the .For example, the environment of the radio target could be a source of absorption (such as in the case of radio jet activity; Nyland et al. 2020;Patil et al. 2021Patil et al. , 2022)).Spectral indices steeper than what is expected from optically thin synchrotron emission may suggest ageing in the system, though further observations are required to fully characterize the underlying mechanism.
The spectral shape classifications were also compared to the radio morphologies investigated in Section 5.No preferential radio spectral classification is found for any of the radio morphologies investigated in this paper, nor for the sample as a whole.Overall, however, the sample is diverse in radio spectral shape.One possible alternative explanation for the unresolved targets is star+quasar superposition.To assess this possibility, we retrieved the SDSS spectra for the targets in this sample and searched the spectra for z = 0 absorption lines using Bayesian AGN Decomposition Analysis for SDSS Spectra (BADASS; Sexton et al. 2021) package. 4BADASS is an open-source spectral analysis program designed for deconvolution of AGN and host galaxy spectra via simultaneous fitting for a variety of spectral features, including but not limited to: the AGN power law continuum, individual broad, narrow, and absorption lines, and stellar host templates or (alternatively) the stellar line-of-sight velocity distribution.
For each target in the sample, we began by fitting host stellar population templates to the observed SDSS spectrum around common z = 0 stellar absorption features (Ca triplet, CaT, rest-frame 8500 Å, 8544 Å, and 8664 Å; Ca H+K lines, rest-frame 3934 Å and 3969 Å; NaD, restframe 5890 Å and 5896 Å, MgIb wavelength; Hα, restframe 6563 Å; Hβ, rest-frame 4861 Å).We did not incorporate any fitting of emission or absorption lines or the quasar continuum during these initial 'quick-look' fits.Each fit was then visually inspected for signs of absorption lines.In cases where absorption lines were identified, we swapped the host stellar template for an AGN power law continuum (to the model the background quasar continuum) and included narrow absorption lines to the fit of the observed spectrum; we centered these absorption lines on the expected z = 0 rest-frame wave- lengths of the lines (see above), left their amplitudes free to vary, and then refit the spectrum.Through this iterative fitting process, we found evidence for CaT in J155859.13+282429.3,J172308.14+524455.5, and J210947.09+065634.7,NaD lines in J024131.89-053139.6,J155859.13+282429.3,J172308.14+524455.5, and J210947.09+065634.7, and Hα in J172308.14+524455.5, but we do not find evidence for MgIb, Hβ, or Ca H+K lines in any targets.The lack of evidence for MgIb and Ca H+K lines is not particularly surprising; the rest-frame wavelengths for the Ca H+K lines reside in a spectral region of low throughput in SDSS while MgIb is not a particularly strong line relative to other lines.Figure 4 displays examples of our absorption line fits to the CaT lines in J155859.13+282429.3,J172308.14+524455.5, and J210947.09+065634.7 and the NaD lines in J024131.89-053139.6,along with inset panels showing the full SDSS spectra.These results demonstrate that 4 systems in this sample comprise quasar-star pairs caught in projection.Assuming no other systems within this sample comprise quasar-star pairs, our sample has a quasarstar contamination rate of 24.2 +10.4 −8.7 % (estimated using scipy.special.betaincinv).However, the absence of foreground stellar absorption lines does not necessarily preclude the presence of other quasar-star pairs within this sample; it has been shown even in the local Universe that the continuum of background AGNs can easily dominate the stellar continuum of foreground stars (even when the sources are separated by ∼ 6 ′′ , Pfeifle et al. 2023), and foreground, z = 0 stellar absorption features can be either diluted or completely elusive.Our quasar-star pair contamination rate is therefore a lower limit for this sample.

Dual Quasars versus Lensed Quasars
One of the primary goals of this study is to identify candidate dual AGN systems.For the portion of the sample identified as multi component targets, it is necessary to eliminate both star+quasar superpositions and lensed quasars contaminants.After eliminating foreground stellar contaminants, Shen (2021) argued that the abundance of high-redshift sub-arcsecond gravitationally lensed quasars is insufficient to account for most of the resolved pairs in a systematic search.This indicates that, while we can expect some additional lensed interlopers, the majority of our remaining sample will be bonafide AGN pairs.Four targets in the pilot sample have been identified as lensed quasars (two are shown in Figure 6: J1415+1129; left, and J0951+2635; right).In J1415+1129, the unresolved but extended source at 3 GHz resolves into at least four individual components at 10 GHz.Archival HST observations reveal that the target is a lensed quasar, an Einstein Cross type system known colloquially as the Cloverleaf System (Chartas et al. 2004).In-terestingly, the radio and optical peaks are not aligned, as was similarly observed at 8.4 GHz (Kayser et al. 1990), where an offset between the optical and modeled radio centroid was measured.However, CO(7-6) observations display good spatial alignment (Alloin et al. 1997) and a recent multi-wavelength study reports the discovery of a radio lobe (Zhang et al. 2022).
In J0951+2635, the double lens is resolved at both new radio frequencies.Archival HST observations reveal a similar optical structure.The target was identified as a lensed quasar via spectroscopy (Schechter et al. 1998), and followed-up by the GLENDAMA (Gil-Merino et al. 2018) and COSMOGRAIL (Sluse et al. 2012) projects.Finally, follow-up HST observations (Jakobsson et al. 2005) confirm the identification.
J0818+0601 was identified by Hutsemékers et al. (2020) as a gravitationally lensed quasar using spectropolarimetric observations.They detect an additional  foreground absorption system tentatively identified as the lensing galaxy.Finally, J1128+2402 has been identified by Inada et al. (2014) as a gravitationally lensed quasar based on the similar spectral energy distributions of the stellar components, in addition to the existence of an extended object in between the measured components.
For the other multi-component and unresolved targets, a lens cannot be ruled out.One possible avenue of further investigation would be follow-up high-resolution integral-field (IFU) spectroscopy of both peaks (Rusu et al. 2019;Lemon et al. 2018).While a nearly identical optical spectrum does not necessarily preclude a dual AGN morphology (given the general similarity of quasar spectra), similar changes in luminosity over time (requiring multi-epoch IFU) might be indicative of a lensed system (e.g.Ciurlo et al. 2023;Lemon et al. 2020Lemon et al. , 2023)).Another method might be to search for possible lensing galaxies in deeper, high angular resolution imaging of these systems (Lemon et al. 2018;Lemon et al. 2023).

Individual Targets
In this section, we present a more in-depth study of the existing ancillary multi-frequency data that exists for each target in the sample, in conjunction with the new VLA observations, radio spectral analysis, and SDSS spectra analysis.Final radio images are presented in Appendix A, Figures 8, 9, and 10.SDSS011114.41+171328.5:exhibits multiple bright components in both radio bands, with a separation of 0.94 arcseconds, or 7.96 kiloparsecs.This target could be tentatively identified as a candidate dual AGN with a relatively flat spectral index.However, we note that the flux ratio between the primary and secondary peaks is quite large, ∼ 20, possibly indicative of jet activity.SDSS024131.89-053139.6: is an unresolved point source at sub-arcsecond scales, with perhaps some slight extension to the NE.The SDSS spectrum for this target reveals evidence of a z = 0 stellar absorption Na doublet, and has thus been identified as a star+quasar superposition.
SDSS074922.97+225511.8:exhibits two radio peaks in both VLA bands, with a separation of 0.46 arcseconds, or 3.91 kiloparsecs.It is another candidate dual AGN, and has been previously identified as such by Chen et al. (2022a), with dual-band HST observations (see Figure 6; middle), in addition to Gemini GMOS spectra.Shen (2021) reported a separation of 3.8 kiloparsecs, essentially matching the separation measured from the new VLA observations.Finally, Shen (2021) reports initial follow-up VLBA observations that confirm the existence of and indicate morphological differences in the dual cores.SDSS081830.46+060138.0: is a faint, unresolved point source in both VLA bands.It has been identified by Hutsemékers et al. (2020) as a doubly imaged, gravitationally lensed quasar, and exhibits an upturned radio spectral shape with a steep spectrum component at low frequencies that may be caused by a previous phase of activity superimposed on more recent, ongoing activity from a compact region driving the flat-spectrum component that emerges at higher frequencies.Alternatively, the spectral shape may be tied to temporal variability and/or lensing effects.SDSS095031.63+432908.4:exhibits a single, extended radio source at 3 GHz that resolves into two component peaks at 10 GHz, with a separation of 0.34 arcseconds, or 2.86 kiloparsecs.We identify this target as another candidate dual AGN.The radio spectrum for this target is peaked, indicating SSA, with a peak frequency of 0.37 GHz, and a curved spectral index of 0.44.SDSS095122.57+263513.9:exhibits two radio peaks at both VLA bands.It has been identified as a lensed quasar via spectroscopy (Schechter et al. 1998) and followed-up by the Gravitational Lenses and Dark Matter project (GLENDAMA; Gil-Merino et al. 2018;Sluse et al. 2012).Follow-up HST (Jakobsson et al. 2005) observations confirm that this system comprises a lensed quasar with a separation of 1.11 arcseconds (see Figure 6; right).As the two components of this target are resolved in both bands, it is possible to calculate spectral indices for each component, in contrast to measuring the overall system.The northwest target exhibits a spectral index of 0.29, while the southeast target exhibits a spectral index of 0.16.SDSS104406.33+295900.9: is an unresolved point source at sub-arcsecond scales.It is bright in both bands, and has previous VLBA observations from the VLBA Imaging and Polarimetry Sky Survey at 5 GHz (VIPS; Helmboldt et al. 2007).Figure 5 shows the archival VLBA image, displaying a bright radio source still consistent with an unresolved point source at milliarcsecond scales.The radio spectrum for this target is likely curved, though there are archival flux densities significantly offset from the model that might indicate flux density variability over years-long timescales.Recently, Wang et al. (2023) found, from new e-MERLIN observations of this target, a significant and large Gaiaradio offset.SDSS112818.49+240217.4:exhibits two radio peaks at 3 GHz.The eastern peak is not visible at 10 GHz, leaving only the western peak remaining.This target VIPS 5 GHz VLBA observations of SDSS104406.33+295900.9 with a 0.0029 ′′ × 0.0018 ′′ beam.
has been identified by Inada et al. (2014) to be a doubly imaged, gravitationally lensed quasar.SDSS121544.36+452912.7:features a pair of canonical radio jets, between which exists a compact core.The northern jet is diffuse, with a separation from the core of 8.80 arcseconds, or about 74 kiloparsecs.The southern jet exhibits a knot at 7.55 arcseconds from the core, or about 64 kiloparsecs, and then continues on to a diffuse lobe at 11.53 arcseconds, or about 97 kiloparsecs.The spectral index of the target's core is ∼ -0.5, while for the northern lobe, it is -0.73, and for the southern lobe, it is -1.15.The target thus follows the canonical model of a radio jet with a flatter core with two steeper radio lobes.
SDSS141546.24+112943.4: is an unresolved point source at 3 GHz that resolves into at least four individual components at 10 GHz.The end-to-end extension of the target is 1.73 arcseconds, or about 14 kiloparsecs.Archival HST images reveal that the target is a lensed quasar (see Figure 6; left); an Einstein Cross known colloquially as the Cloverleaf System (Magain et al. 1988;Chartas et al. 2004).The target's radio spectral shape is best fit to a standard power law, with a spectral index of -1.2, slightly steeper than might be expected for a gravitationally-lensed quasar.SDSS143333.02+484227.7:exhibits a jet-like extension to the west of the radio core, and a smaller extension to the east.The extension is relatively small-scale: measured from end-to-end, the full length of the target is 1.39 arcseconds, or about 12 kiloparsecs.For this target, there is a secondary faint (magnitude of ∼21) Gaia detection close to the primary source (∼0.767 ′′ ).The secondary source is ∼18% of the brightness of the primary source.SDSS155859.13+282429.3: is an unresolved point source in both bands.The SDSS spectrum for this target reveals evidence for both a Ca triplet and a Na doublet absorption system, and has thus been identified as a star+quasar superposition.This target is the most obvious star-quasar pair; the SDSS spectrum also exhibits prominent TiO absorption bands, which are a common characteristic of M dwarfs (Morgan et al. 1943).
SDSS162501.98+430931.6:exhibits a single, extended radio source at 3 GHz, that resolves into two radio components at 10 GHz.This target is another candidate dual AGN, with a separation of 0.52 arcseconds, or 4.52 kiloparsecs.
SDSS172308.14+524455.5: is an unresolved point source at sub-arcsecond scales, with no clear extension in the fairly bright radio signature.The SDSS spectrum for this target reveals evidence for a z = 0 Ca triplet, Na doublet, and Hα absorption lines, and has thus been identified as a star+quasar superposition.
SDSS210947.09+065634.7: is an unresolved point source at sub-arcsecond scales.The SDSS spectrum for this target reveals evidence of a z = 0 Na doublet and a Ca triplet absorption line system, and has thus been identified as a star+quasar superposition.The radio spectrum for this target is curved, with a peak frequency of 1.67 GHz and a curved spectral index of 0.46.Recently, Wang et al. (2023) found, from new e-MERLIN observations of this target, a significant and large Gaiaradio offset.

Understanding the Varstrometry-selected Radio Parameters
In addition to understanding each individual source, it is equally important to understand the sample that the varstrometry method selects for.In order to characterize the overall target sample, a control sample was prepared.The control sample was generated alongside the target sample, and systems identified as star+quasar superpositions were identified by their SDSS spectra and were eliminated from the target sample.
The control sample was generated from a crossmatch of the SDSS DR16Q and Gaia EDR3 cata-logs, to within 1.5 ′′ .The sample Gaia G-magnitude (G < 20) and redshift (z > 0.5) cuts were applied.In this case, the important Gaia parameter, the astrometric excess noise significance was set to < 4, in order to generate a control sample of sources that do not display significant astrometric jitter.Finally, it was required that each control source have a VLASS signature so that overall control and target sample radio parameters could be compared.The control sample was then spatially matched to the target sample, so as to ensure that the controls covered the same portion of the sky as the targets to minimize systematics arising from the Gaia scanning law.The number of controls per target object was maximized and required to be the same for each target.Additionally, no controls were shared between targets, resulting in a control sample of about 120 sources.
The controlled parameters, Gaia G-magnitude and redshift, show probability values calculated using the Anderson-Darling, Kolmogorov-Smirnov, and Kuiper's tests that are >0.05,indicating that the null hypothesis (in this case, that both the target and control samples arise from the same parent population) cannot be ruled out.The same is true for the Gaia G-magnitude signalto-noise ratio, the peak radio luminosities, the total radio luminosities, the ratio of peak-to-total radio luminosity, and radio spectral indices.For the controlled parameter AENS, the p-values are <0.05,rejecting the null hypothesis, and indicating that the two samples vary by more than the parameter of interest.
Distributions are shown in Figure 7, illustrating the offset in the AENS.The controlled parameters exhibit the relationship driven by the sample selection parameters, but the important radio parameters (peak radio luminosity, total radio luminosity, ratio of peak-to-total radio luminosity, and radio spectral index) illustrate that the target sample is not radio-loud when compared to the matched control sample, likely eliminating the possibility of blazar jet activity driving a large percentage of the high AEN.The peak and total radio luminosities are those measured from VLASS.The spectral index is calculated between measurements from FIRST (at 1.4 GHz) and measurements from VLASS (at 3 GHz).The purpose of these plots is to compare the control and target samples for the important radio parameters.Gaia G-band magnitude, G-band magnitude SNR, and target redshift were controlled for in the creation of the control sample, and the p-values listed on the plots show the lack of disagreement between the two samples.Note: In Panel 4, the logarithmic values of AENS + 1 are used.In Panels 5 and 6, the logarithmic values of the peak radio luminosity are used.

SUMMARY AND CONCLUSIONS
We present deep, targeted VLA observations of 18 Gaia-unresolved quasars identified as astrometricallyvariable.
The high redshifts (0.8 ≤ z ≤ 2.9) of these quasars probe an interesting observational gap in the current population of confirmed and candidate kiloparsec-scale dual AGN near cosmic noon.This pilot study, which combines the varstrometry method (Hwang et al. 2020) with high-resolution radio observations, is the first in a series of studies designed to fully understand the radio+varstrometry sample.Our conclusions are as follows: • The new VLA observations reveal eight targets identified as multi-component (candidate dual AGN) or gravitational lenses, or ∼44% of the overall sample.Four of the eight are currently identified as gravitational lenses.The result is comparable to that of Chen et al. (2022a).Optical/infrared spectroscopy and imaging will be required to filter out other contaminants by identifying stellar features, modeling lensed quasar systems, and identifying lensing galaxies.
• Three of the 18 targets have been identified as systems exhibiting jet activity or otherwise extended emission.This population includes J1415+1129, which has been identified from existing HST observations as a lensed quasar.
• A total of four of the 18 targets have been identified as gravitationally lensed quasars.We note that, while the focus of this search was multiple AGN systems, gravitationally lensed systems are equally interesting, and are the focus of many ongoing studies, including in the high resolution radio regime (Spingola et al. 2019;Casadio et al. 2021).Additionally, gravitational lenses remain a focus of those mining Gaia's rich data set in order to create gravitational lens catalogs (Lemon et al. 2023;Lemon et al. 2018).Follow-up high resolution optical observations will be required to identify gravitational lenses with smaller separations.
• Nine of the 18 targets have unresolved radio signatures.After careful analysis of the SDSS spectra, four of these nine were identified as superpositions of foreground stars and background quasars, and thus were eliminated as candidate dual AGN.J0818+0601 was additionally identified as a gravitationally lensed quasar.There remain many astrophysical explanations for the excess astrometric noise exhibited by the remaining unresolved targets.Follow-up multiwavelength observations will be required to observe any existing structure at smaller scales and to confirm or reject all scenarios.For the AGN pair scenario, this includes highresolution optical/infrared imaging with HST or JWST to identify the two quasars and radio interferometric observations at smaller scales to identify two compact cores with the VLBA or the nextgeneration Very Large Array (ngVLA; Selina et al. 2018;Nyland et al. 2018).
• It is also possible for the unresolved targets to exhibit high astrometric variability caused by the properties of the host galaxy.If the host structure is significantly extended, as in the case of tidal features, this could contribute to the excess astrometric noise (Makarov & Secrest 2022, 2023).Despite the quality of the DECaLS imaging, due to the high redshifts of the targets, it is not possible to constrain the host galaxies with currently available optical observations from SDSS and DECaLS.Quasar host galaxy characterization will require follow-up optical/infrared imaging with higher angular resolution.
• We note that the excess astrometric noise could also be attributable to Gaia systematics, though this is unlikely given the careful sample selection.Future Gaia data releases will provide increasingly better overall astrometric measurements.Additionally, Gaia's pipeline continues to develop better treatment of extended sources, and continues to increase the reliability of the astrometric excess noise parameter, increasing the number of AGN pairs identifiable in future data releases.
• An in-depth study of the overall sample illustrated that the targets sample is not particularly radio loud in comparison to a matched control sample.This likely eliminates the possibility of blazar jet activity as the main driver of the high excess astrometric noise.
• Finally, a thorough review of the radio spectral shapes of each target reveals that the majority of the targets, no matter their spectral shape classification, exhibit a spectral index that is consistent with that of optically thin synchrotron emission.Targets with curved radio spectral shapes would benefit from follow-up simultaneous observations with broad frequency coverage above and below the peak frequency.
Overall, we demonstrate the potential of using ra-dio+varstrometry to systematically discover genuine kiloparsec-scale dual quasar candidates, including at redshifts greater than z ∼ 0.8.Additionally, these results will assist in refining the targeting strategy to facilitate a more efficient identification of dual and binary AGN specifics by providing radio constraints that can be applied to a larger, more systematic search of existing radio surveys.] SDSS J210947.09+065634.7 Figure 11.Observed SDSS spectra for the entire sample of 18 targets.These panels display the observed SDSS spectrum (∼4000 Å to ∼9000 Å) available for each system in the sample; these spectra are not redshift corrected.The designation for each system is listed in the top left corner of each panel.As noted in Section 7.1, four of these systems likely comprise a foreground star and background quasar caught in projection (J0241-0531, J1558+2824, J1723+5244, J2109+0656).Notice in several cases that the sources appear elongated rather than point-like, which may potentially point to underlying complex host morphologies and/or a multiplicity of sources, neither of which are not resolvable in these images.
Figure 13.Broadband radio spectra of unresolved targets showing the new, quasi-simultaneous VLA observations in purple squares, and the following archival measurements: VLASS (cyan, right-facing triangle), LOTSS (orange circle), FIRST (blue cross), VLITE (brown, left-facing triangle), NVSS (pink pentagon), VLBA (grey diamond), TGSS (hot pink star), and RACS (green, downward-facing triangle).For each source, the spectrum has been modeled with both a standard power law (red line) and a curved power law (dashed black line).The resultant spectral indices are listed.

Figure 1 .
Figure 1.The Varstrometry Technique: Top: A cartoon illustrating the stochastic variability of a dual AGN including the shifting of the light center between the two component AGN.Bottom: As an example, an AGN photometric lightcurve simulation, as seen by Gaia, described by a first-order continuous autoregressive process (e.g.damped random walk; MacLeod et al. 2010) with τ = 10 day.A typical variability amplitude for AGN of about 0.1 magnitudes was assumed (e.g.Sesar et al. 2011), as well as a pair angular separation of 30 milliarcseconds.The middle plot shows a joint variability lightcurve, which in itself indistinguishable from systems where the apparent photocenter of the AGN is different from that of its host galaxy.The bottom plot shows the astrometric "jitter:" the photocenter of the pair appears to wander back and forth.Top figures adapted from Shen (2021).

Figure 2 .
Figure 2. Target Selection: A flowchart summarizing the process of sample selection.
The new VLA images are presented in Appendix A, Figure 10.
6.1.4.RACSThe Rapid Australian Square Kilometer Array Pathfinder (ASKAP) Continuum Survey (RACS;Mc- Connell et al. 2020)  is the first all-sky survey conducted with ASKAP.The first epoch observations at 887.5 MHz use 36 12-m antennas to cover the entire sky south of +41 degrees, reaching a sensitivity of between 25 and 40 µJy per beam and a beam of about 15 ′′ .A calibration uncertainty of 5% was applied.6.1.5.NVSSThe NRAO VLA Sky Survey (NVSS;Condon et al. 1998) observed the entire northern sky about δ > -40 degrees at 1.4 GHz.The survey was conduction with the VLA in D configuration, resulting in an angular resolution of about 45 ′′ .A calibration uncertainty of 3% was applied.6.1.6.FIRSTThe Faint Images of the Radio Sky at Twenty-one centimeters (FIRST;Becker et al. 1995) survey is a 1.4

Figure 3 .
Figure 3. Radio color-color plot showing two-band spectral indices α 3 1.4 and α 10 3 , calculated using flux densities from the new 3 GHz and 10 GHz VLA observations, in addition to FIRST (1.4 GHz; VizieR 2014 catalog Helfand et al. 2015) values.Where no FIRST catalog source was recorded, NVSS (1.4 GHz; VizieR catalog Condon et al. 1998) flux densities were used.If neither survey was available, an upper limit was used of NVSS 1σ 0.45 mJy (Condon et al. 1998).The different symbols highlight classes of interest: star+quasar superposition systems are marked with stars, gravitational lenses are marked with diamonds, and multi-component systems (the most likely to be candidate AGN pairs) are marked as triangles.All remaining targets are marked with circles.The arrows used represent lower limits on the spectral indices due to the use of NVSS upper limit flux densities.readerto Section 12 for detailed spectral modeling on each individual source.There are many processes that might explain genuine absorption exhibited by some targets (i.e.synchrotron self-absorption (SSA; de Kool & Begelman 1989; Tingay & de Kool 2003; O'Dea & Saikia 2021), or free-free absorption(FFA;Kellermann 1966;O'Dea 1998;Bicknell et al. 1998)).For example, the environment of the radio target could be a source of absorption (such as in the case of radio jet activity;Nyland et al. 2020;Patil et al. 2021 Patil  et al.  , 2022)).Spectral indices steeper than what is expected from optically thin synchrotron emission may suggest ageing in the system, though further observations are required to fully characterize the underlying mechanism.The spectral shape classifications were also compared to the radio morphologies investigated in Section 5.No preferential radio spectral classification is found for any of the radio morphologies investigated in this paper, nor for the sample as a whole.Overall, however, the sample is diverse in radio spectral shape.
. Chance Superposition with Stars

Figure 4 .
Figure 4. SDSS optical spectroscopic absorption line fitting.Here we show the four cases of confirmed foreground starbackground quasar pairings (top: J0241-0531, top-middle: J1558+2824, bottom-middle: J1723+5244, bottom: J2109+0656) where we have identified foreground stellar absorption features consistent with z = 0.Each panel shows (1) the fits to the spectrum across a narrow spectral window straddling the observed absorption lines, (2) the residuals of the fit, and (3) an inset panel showing the full SDSS spectrum.Black lines represent the observed SDSS spectrum (uncorrected for the redshift of the quasar), red lines represent the model (quasar power-law continuum + gaussian absorption lines).We have marked the observed absorption lines with dashed vertical lines and denoted the specific line by name.The full SDSS spectra for these four likely star-quasar pairs are shown in Appendix 9 along with the SDSS spectra of the rest of the sample.

Figure 7 .
Figure 7.These histograms compare the control sample of 120 targets to a target sample of 12 targets (all with AENS > 5).The peak and total radio luminosities are those measured from VLASS.The spectral index is calculated between measurements from FIRST (at 1.4 GHz) and measurements from VLASS (at 3 GHz).The purpose of these plots is to compare the control and target samples for the important radio parameters.Gaia G-band magnitude, G-band magnitude SNR, and target redshift were controlled for in the creation of the control sample, and the p-values listed on the plots show the lack of disagreement between the two samples.Note: In Panel 4, the logarithmic values of AENS + 1 are used.In Panels 5 and 6, the logarithmic values of the peak radio luminosity are used.

Figure 8 .
Figure 8. New VLA observations for all targets identified as unresolved.

Figure 8 .
Figure 8. Continued from previous page.New VLA observations for all targets identified as unresolved.

Figure 8 .
Figure 8. Continued from previous page.New VLA observations for all targets identified as unresolved.

Figure 9 .
Figure 9. New VLA observations for all targets identified as having multiple components.

Figure 9 .
Figure 9. Continued from previous page.New VLA observations for all targets identified as having multiple components.

Figure 10 .
Figure 10.New VLA observations for all targets identified as extended.

Figure 12 .
Figure12.Optical tricolor ugz images drawn from the DeCaLs Legacy Viewer for the full sample.The SDSS designation for each object is listed in the top left corner of each panel, while the scale bar in the bottom right corner indicates 5".Notice in several cases that the sources appear elongated rather than point-like, which may potentially point to underlying complex host morphologies and/or a multiplicity of sources, neither of which are not resolvable in these images.

Table 1 .
Parameter Details

Table 2 .
Target Details

Table 3 .
Observation Details

Table 4 .
Sample Properties for Unresolved Sources

Table 5 .
Sample Properties for Multi-component Sources

Table 5 .
Note -Column 1: Coordinate names in the form of "hhmmss.ss±ddmmss.s"based on SDSS DR16Q are bolded and starred, and individual component location as measured from 10 GHz VLA observations are listed beneath.Column 2:

Table 6 .
Sample Properties for Extended Sources

Table 7 .
Tapered Image Properties

Table 7 .
Note -Column 1: Coordinate names in the form of "hhmmss.ss±ddmmss.s"based on SDSS DR16Q.Column 2: Restoring beam of uv-tapered imaging: major axis, minor axis, position angle.All imaging done with Briggs weighting and a robust factor of 1. Column 3: 10 GHz uv-tapered image source flux density ± error.Column 4: 10 GHz uv-tapered image 1-σ noise.

Table 8 .
Radio Spectral Modeling Properties

Table 8 .
Note -Column 1: Coordinate names in the form of "hhmmss.ss±ddmmss.s"based on SDSS DR16Q.Column 2: VLA 3 GHz (S-band) total flux density ± error.Column 3: VLA 10 GHz (X-band) uv-tapered total flux density ± error.Column 4: Spectral index from least squares fitting of radio spectral fitting with standard power law ± error.Column 5: Spectral index from least squares fitting of radio spectral fitting with curved power law ± error, where possible.Column 6: Classification.Column 7: Spectral shape classification.