Varstrometry for Off-nucleus and Dual Sub-kpc AGN (VODKA): Very Long Baseline Array Searches for Dual or Off-nucleus Quasars and Small-scale Jets

Dual and off-nucleus active supermassive black holes are expected to be common in the hierarchical structure formation paradigm, but their identification at parsec scales remains a challenge due to strict angular resolution requirements. We conduct a systematic study using the Very Long Baseline Array (VLBA) to examine 23 radio-bright candidate dual and off-nucleus quasars. The targets are selected by a novel astrometric technique ("varstrometry") from Gaia, aiming to identify dual or off-nucleus quasars at (sub)kilo-parsec scales. Among these quasars, 8 exhibit either multiple radio components or significant (>3$\sigma$) positional offsets between the VLBA and Gaia positions. The radio emission from the three candidates which exhibit multiple radio components is likely to originate from small-scale jets based on their morphology. Among the remaining five candidates with significant VLBA-Gaia offsets, three are identified as potential dual quasars at parsec scales, one is likely attributed to small-scale jets, and the origin of the last candidate remains unclear. We explore alternative explanations for the observed VLBA-Gaia offsets. We find no evidence for optical jets at kilo-parsec scales, nor any contamination to Gaia astrometric noise from the host galaxy; misaligned coordinate systems are unlikely to account for our offsets. Our study highlights the promise of the varstrometry technique in discovering candidate dual or off-nucleus quasars and emphasizes the need for further confirmation and investigation to validate and understand these intriguing candidates.


INTRODUCTION
Most massive galaxies are believed to harbor a central supermassive black hole (SMBH; Kormendy & Richstone 1995;Kormendy & Ho 2013), and as a result, binary SMBHs are expected to form from galaxy mergers (Begelman et al. 1980).Binary SMBHs are of ut-Corresponding author: Yu-Ching Chen ycchen@jhu.edumost importance to the studies of gravitational waves, black hole astrophysics, and galaxy evolution.Despite successful numerical simulations that have modeled the dynamical evolution of binary SMBHs down to parsec (pc) scales (Fiacconi et al. 2013;Mayer 2013;Souza Lima et al. 2017;Tamburello et al. 2017), direct observational evidence for binary SMBHs and sub-pc inspiraling SMBHs remains scarce.To date, the only known confirmed compact binary SMBH at pc scales is 0402+379, a flat-spectrum compact double radio source at z = 0.055 with a projected separation of 7 pc, which was discovered serendipitously (Rodriguez et al. 2006;Bansal et al. 2017).
A binary SMBH's orbital decay may slow down or halt at pc scales, a process commonly referred to as the "final pc problem", because close encounters with stars are insufficient to remove orbital energy (Begelman et al. 1980;Milosavljević & Merritt 2001;Yu 2002).While interactions with gas, triaxial galaxy structure, or a third SMBH may overcome this barrier (Gould & Rix 2000;Hoffman & Loeb 2007;Kelley et al. 2017), how efficiently binary SMBHs evolve from kpc to pc-scale separations remains uncertain.This uncertainty creates ambiguity in predicting the intensity of low-frequency gravitational waves by current facilities such as pulsar timing array (Arzoumanian et al. 2018).To address this issue, we need to observe the rate of dual SMBHs at subkpc scales and compact binaries at ∼pc scales.These observations could provide direct measurements and insights into final-pc physics and the interpretation of future low-frequency gravitational-wave detections (Bogdanović et al. 2022).
Quasars are luminous sources of electromagnetic radiation powered by the accretion of matter onto the SMBH in active galactic nuclei (AGNs).Following the eventual merger of a binary SMBH, the coalesced SMBH may receive a kick due to the anisotropic gravitational wave radiation because of momentum conservation (Baker et al. 2006;Campanelli et al. 2007).If the recoiled SMBH accretes from the ambient gas, it may emit radiation and be observed as an off-nucleus quasar (Loeb 2007).The observed statistics of off-nucleus quasars can be used to investigate the fueling processes of quasars in galaxy mergers, as well as the population of recoiling SMBHs that result from binary SMBH coalescence (Blecha et al. 2016;Tremmel et al. 2018).
The identification of sub-kpc dual/off-nucleus quasars has been challenging due to the strict spatial resolution required for their detection (Chen et al. 2022).Although an increasing number of sub-kpc dual AGN are being discovered serendipitously (Woo et al. 2014;Müller-Sánchez et al. 2015;Goulding et al. 2019;Koss et al. 2023), the systematic identification of sub-kpc dual/offnucleus quasars has remained inefficient.
We have developed a novel astrometric approach called "varstrometry" to identify sub-kpc unresolved dual/off-nucleus/lensed quasars throughout the entire sky (Hwang et al. 2020;Shen et al. 2019; see also Williams & Saha 1995;Liu 2015;Springer & Ofek 2021;Mannucci et al. 2022).With Gaia's coverage of the full sky and its depth reaching down to Gaia G∼21 mag (encompassing ∼10 6 AGN, Carnerero et al. 2022), varstrometry enables a systematic exploration of un-resolved dual/off-nucleus AGNs in the poorly studied range at projected physical separations of ∼10-1000 pc.The technique's effectiveness has been established using Gaia DR2 data on a sample of variable pre-main sequence stars, which have known unresolved close companions (Hwang et al. 2020).Hubble Space Telescope (HST) high-resolution dual-band optical imaging discovered dozens of quasars with kpc-scale double cores from varstrometry-selected targets (Chen et al. 2022), some of which were subsequently verified as kpc-scale dual or lensed quasars through extensive follow-up observations (Shen et al. 2021;Chen et al. 2023).Nonetheless, resolving sources at pc spatial scales remains challenging due to the angular resolution restriction of most existing facilities.Long-baseline radio or optical/infrared interferometers are a few instruments capable of directly imaging binary SMBHs with an angular resolution of a few milliarcseconds (e.g., Rodriguez et al. 2006;Burke-Spolaor 2011;Wrobel et al. 2014;GRAVITY Collaboration et al. 2017;Gravity Collaboration et al. 2018;Liu et al. 2018;Event Horizon Telescope Collaboration et al. 2019, 2022).
While radio interferometers are excellent tools to resolve sources that are not resolvable as binaries in optical surveys, finding two radio-loud sources in the same pair is statistically difficult if radio loudness is not affected by merging, because only ∼20% of quasars are radio-loud (Kellermann et al. 2016).Previous systematic searches for binary SMBHs with with Very Long Baseline Interferometer (VLBI) have found a few candidates (Burke-Spolaor 2011;Tremblay et al. 2016;Breiding et al. 2021).However, only one confirmed binary SMBH at pc scales, 0402+379, is known to date (Rodriguez et al. 2006).Gaia astrometry in comparison with VLBI astrometry is useful for studying jet systems of quasars (Kovalev et al. 2017;Petrov & Kovalev 2017;Makarov et al. 2017;Plavin et al. 2019;Petrov et al. 2019) or discovering binary/off-nucleus quasars (Skipper & Browne 2018;Breiding et al. 2021).A recent study using varstrometry with e-MERLIN revealed four targets with Gaia-radio offsets of ∼ 9-60 mas (Wang et al. 2023), highlighting the potential of varstrometry to uncover candidate dual/off-nucleus quasars.In addition to binary SMBHs, VLBI was also used to search for strong gravitational lenses on milliarcsec scales (Casadio et al. 2021), which might be detected by varstrometry.
In this study, we present observations of 23 radiobright quasars at redshifts 0.5<z<2.5,obtained with VLBA at X, C, and/or Ku-bands.These quasars were selected using Gaia data and the varstrometry technique (Shen et al. 2019;Hwang et al. 2020).While our earlier varstrometry papers (Shen et al. 2021;Chen et al. 2022)   focused on dual quasars at kpc scales, varstrometry can in principle probe ∼10-100 parsec scales.We discover six candidates with significant Gaia-VLBA offsets and three candidates with multiple radio cores at physical scales of ∼10-100 pc, which could be interpreted as binary quasars, off-nucleus quasars, or quasars with smallscale jets (see Figure 1).In this paper, we describe our methodology and sample selection in §2, and our VLBA observations, data reduction, and analysis in §3.Our main results are presented in §4, and we discuss their implications in §5.We summarize our findings and discuss future prospects in §6.Throughout this paper, all physical separations refer to projected separations.We adopt a flat ΛCDM cosmology with Ω Λ = 0.7, Ω m = 0.3, and H 0 = 70 km s −1 Mpc −1 .

SAMPLE SELECTION
The selection of dual/offset quasar candidates relies on a novel astrometric technique, varstrometry, which exploits the centroid jitters resulting from a quasar's brightness variability (Liu 2015;Shen et al. 2019;Hwang et al. 2020).In the case of a dual quasar, these jitters arise from the non-coherent variability of each nucleus, causing the centroid of the total flux to shift.For an offnucleus quasar, these jitters result from its brightness changing relative to a non-variable host galaxy, causing the centroid of the total flux to shift away from the center of the host galaxy as the flux of the quasar brightens.These astrometric jitter signals, induced by variability, can be detected by Gaia as "excess astrometric noise" beyond the best-fit astrometric solution, or even as an anomaly in parallaxes and proper motions (Shen et al. 2019;Hwang et al. 2020;Chen et al. 2022).
We focus on optically unobscured, broad-line quasars, for which the varstrometry technique is applicable.Our study began with a sample of 356,850 spectroscopically confirmed quasars from Sloan Digital Sky Survey (SDSS) Data Release 14 (DR14, Pâris et al. 2018) that have single matches within 2 ′′ radius in Gaia Data Release 2 (DR2).Among these, we identified 3,735 objects with peak flux densities >15 mJy in Faint Images of the Radio Sky at Twenty-Centimeters (FIRST, Becker et al. 1995).We then selected sources with significant Gaia DR2 astrometric excess noise (astrometric excess noise >3σ at any given Gaia magnitude, where σ is the standard deviation of AEN as a function of G) to obtain 27 targets (Figure 2).AEN selection also helps to avoid strongly variable sources such as blazars (Fig. 12 in Hwang et al. (2020)).We removed interlopers caused by extended host galaxies in PanSTARRS images (Chambers et al. 2016) and star superpositions by inspecting SDSS optical spectra that enclose all fluxes within a 2 ′′ -3 ′′ diameter circle.
Among the 27 candidates, we found that seven have archival VLBA C-band (5 GHz) observations (Taylor et al. 2005;Helmboldt et al. 2007).We removed 1 archival target that lacks an accurate VLBA position from Petrov & Taylor (2011).Of the remaining 20 candidates without archival VLBA images, three were excluded with FIRST-resolved jets, and then we conducted new VLBA observations for the rest 17 candidates.Our final sample consists of 23 targets that were observed in our new VLBA observations (18 targets) or have archival VLBA observations with accurate positions (6 targets).Only our target selection was based on Gaia DR2 AEN and magnitude (Figure 2); all the analyses in this paper are based on the new Gaia DR3 data.Gaia DR3 has notably improved its AEN measurements.The AENs decrease significantly from DR2 to DR3 for many of our targets (Table 1), especially those without multiple radio components and offsets, which would no longer pass the 3σ selection criteria if Gaia DR3 AEN is used.Properties of the 23 targets are listed in Table 1 including coordinate, redshift, magnitude, Gaia astrometric excess noise, Gaia flux variability, VLBA observed band, and FIRST flux density.

New VLBA Observations
Among the 23 targets, 18 (17 objects without VLBA archival images plus J1044+2959) were observed with VLBA in  between February 2020 and September 2020 (Program ID: BL276; PI: X. Liu).Additional observations were conducted in C-band (4.87 GHz) and/or Ku-band (15.2 GHz) for targets showing extended structures and/or offsets from the Gaia positions in the X-band images.The experimental setup for X-band observations included single polarization, 2 s integration times, and a total bandwidth of 512 MHz (16 spectral channels with a bandwidth of 32 MHz).The exposure time was chosen based on 1.5 Ghz FIRST flux density assuming a conservative spectral index of -1 for synchrotron radiation, a factor of 0.3 reduce in flux going from arcsecond to VLBA scales, a flux ratio of 1:5 between potentially resolved binary components, and at least a 6-sigma detection for each component.At least 8 stations were used during all observations.The antenna at FD (Fort Davis, TX) was used as the reference antenna.All targets were observed in the phasereferencing mode, with nearby bright quasars used as the phase-reference calibrators, except for J0840+4439 and J1044+2959, which were self-calibrated.A single bright quasar from the following list was used as the amplitude calibrator for each observing session: 3C454.3,0235+164, 4C39.25, and 3C345.The observation details are summarized in Table 2.

Archival VLBA Observations
Of the 23 targets, 6 (including J1044+2959) have archival VLBA C-band images from the VLBA Imaging and Polarization Survey (VIPS; Helmboldt et al. 2007) and have been re-calibrated with improved positions from Petrov & Taylor (2011).The VIPS observations were conducted in 2006 at a center frequency of 5 GHz, with a typical 1σ rms noise level of 0.2 mJy beam −1 .The improved absolute astrometry from Petrov & Taylor (2011) has a position uncertainty of 0.5 mas.Of the six targets, only J1044+2959 shows a significant VLBA-Gaia offset and has a compact core structure.Therefore, we propose additional X-band and Ku-band observations for J1044+2959 to its radio spectral indices.

Data Reduction and Analysis
The data reduction and analysis were conducted using the Astronomical Image Processing System (AIPS; Greisen 2003).We followed the standard reduction pipeline task vlbarun, which includes corrections to the Earth orientation parameters, ionospheric correction, amplitude correction, and phase correction.Subsequently, we used dbcon to combine the visibility data from different observing sessions for the same targets.Then, we created dirty images and used imagr to de-  (Helmboldt et al. 2007) b The redshift of 6.169 from the SDSS pipeline is not correct.convolve (clean) the images.To optimize the sensitivity, we used a natural weighting scheme.The typical minor axes of the synthesized beams were 1.5 mas for C-band, 1 mas for X-band, and 0.5 mas for Ku-band.The pixel scales were 0.3 mas for C-band, 0.2 mas for X-band, and 0.1 mas for Ku-band, which were chosen to have at least five pixels in the synthesized beam.Additionally, we created wide-field images to search for any radio cores up to 0. ′′ 6, which is the typical resolving power of Gaia.

VLBA Continuum Images
We present the VLBA continuum images for the 23 targets.We divide these targets into two groups based on their characteristics.Figure 3 and Figure 4 display the VLBA images of the 8 targets that exhibit significant (> 3σ) VLBA-Gaia offsets or/and have multiple radio components.Figure 5 shows the VLBA images of the 15 targets that either show consistent positions between VLBA and Gaia, or are undetected.The images in Figure 3, 4, and 5 are 50 mas × 50 mas in size and are centered at the Gaia positions.However, for J1557+1236, we used a larger image size of 150 mas × 150 mas to accommodate the VLBA position.The majority of the images have a typical 1-σ noise level of approximately 0.1 mJy beam −1 .However, for the brighter sources (such as J0749+2255, J1044+2959, and J1137+4825), the 1-σ noise levels reach as high as approximately 1 mJy beam −1 .Further information about the individual images can be found in Table 3.

Source Detections and Properties
Of the 18 targets observed in our new VLBA observations, 16 are detected with a 5-σ detection threshold, and 2 are undetected.We search for sources within a 0. ′′ 6 radius and model them with 2-dimensional Gaussian components using jmfit in AIPS.Each image is fit using a single Gaussian component.For sources with extended emission or sub-structures, we fit only a single Gaussian to the brightest component.Most of tar-    gets exhibit single compact components, while a few targets (J1051+2119, J1110+4817 and J1259+5140) show second compact components or extended structures at the scales of dozens mas.We do not find any compact source at large scale (>100 mas), except for the known dual quasar J0749+2255, which has a separation of 0. ′′ 46 (Shen et al. 2021;Chen et al. 2023).The VLBA radio measurements and image properties for each target are summarized in Table 3.Previous analyses have suggested that differences in the u-v coverage of the observations may contribute up to 10% in systematic un-certainties.The coordinate positions obtained from the C-band and Ku-band images are consistent with those obtained from the X-band images, with differences of less than 0.1 mas, consistent with the astrometric accuracy.The 5-σ upper limits of the two undetected targets, J1318+3842 and J1528+1827, are also provided in Table 3.

Spectral Index
We have determined the spectral indices of the six sources for which we have detections at at least two frequencies.These six targets comprise the four tar-Figure 5. VLBA (C-band and/or X-band) continuum images of the 15 targets that either show consistent positions between VLBA and Gaia, or are undetected.Other information is the same as those in Figure 3 Table 4. Spectral indices, morphology, and classifications of the radio emission for the targets that have multi-band observations, exhibit significant VLBA-Gaia offsets, or show multiple radio components.
gets with significant VLBA-Gaia offsets and the two targets with extended structures.The spectral index, denoted as α, is defined as S ν ∝ ν α , where S ν is the radio flux density and ν is the frequency.We present the spectral indices (Table 4) and the radio spectra (Figure 6) for the six targets.Three targets, namely J0749+2255, J1044+2959, and J1110+3653, exhibit flat or rising spectra with α ≳ 0, suggesting synchrotron selfabsorbed cores.The remaining three targets have steep spectra with α ≲ −1, indicating optically-thin radio jets.
The spectral indices and radio emission classifications of these targets are also listed in Table 4.

VLBA-Gaia offset
We compare the optical positions from Gaia DR3 with the VLBA radio positions to identify any potential dual, offset, or lensed quasar candidates.The VLBA-Gaia offsets are calculated based on the X-band images.The VLBA position uncertainty is determined using the quadrature sum of the fitting error, the phase calibrator's position uncertainty, and the error associated with the phase-referencing technique (Pradel et al. 2006).We obtain the phase calibrator positions from the VLBA calibrator list (Charlot et al. 2020;Petrov 2021).The final VLBA position errors were within the range of 0.3 to 1.2 mas. Figure 7 displays the positional offsets' distribution for the 21 targets detected by the VLBA.Among the VLBA-detected targets, 6 have radio counterparts significantly offset from the Gaia positions using the 3-σ threshold.Further discussion of the origin of radio emission for these VLBA-Gaia offset targets is presented in Section 5.1.

Targets with multiple radio detections or significant VLBA-Gaia offsets
Out of the 16 detected targets, 8 targets exhibit either multiple radio detections or significant VLBA-Gaia offsets (Figure 3 and 4).Based on the spectral indices and morphology, we examine the potential origins of the radio emission for the 8 targets.If a target shows two flat-spectrum cores or a flat-spectrum core with significant Gaia-VLBA offset, we classify it as a dual quasar candidate.If a target shows steep spectral index or extended structure along the Gaia-VLBA offset, the radio emission is likely from a jet. a diffuse lobe structure.Therefore, the offset is most likely caused by a small-scale radio jet.

J0749+2255
J0749+2255 shows a 4.0±0.6 mas (corresponding to 33 pc) VLBA-Gaia offset.The radio component has a inverted spectral index of 0.24, indicating a selfabsorbed core.J0749+2255 was already identified as a 0. ′′ 46 dual quasar (Chen et al. 2022), with both nuclei detected in our VLBA images and reported in Shen et al. (2021).Gaia only reports one detection for J0749+2255, so the Gaia optical position may be affected by the blended second nucleus, leading to the offset.Nonetheless, the possibility of a third SMBH at a pc scale centered at the optical component cannot be ruled out.

J1044+2959
J1044+2959 shows a 3.2±0.5 mas (corresponding to 25 pc) VLBA-Gaia offset.The radio component is compact and has a flat spectral index of −0.13, suggesting a self-absorbed core.J1044+2959 is likely a candidate dual quasar, which could explain the offset if one nucleus is radio-bright but the other is optically bright (Panel b, Figure 1).Alternatively, both quasars could be optically bright, but only one is radio-bright (Panel a, Figure 1), resulting in an offset between the light centroid of the total optical flux from the radio centroid.

J1051+2119
J1051+2119 shows a consistent position between VLBA and Gaia, but it presents a second radio component located toward the east.The second component shows an elongated structure connecting to the primary component, which suggests that it is a small-scale jet.To confirm the jet scenario, a spectral index analysis derived from multi-band images is necessary, though we only have an archival C-band image.

J1110+3653
J1110+3653 exhibits a 3.7±1.0mas (corresponding to 25 pc) VLBA-Gaia offset.The radio component is compact and has a inverted spectrum (α C−X = 2.71 and α X−Ku = −3.02),consistent with a self-absorbed core.Similar to the discussion for J1044+2959, J1110+3653 is likely a candidate dual quasar.In addition, given its low redshift of 0.63, J1110+3653 could be a candidate off-nucleus quasar where the optical center is shifted because of emission from the host galaxy (Panel c in Figure 1).

J1110+4817
J1110+4817 exhibits an 11.9±1.0mas (corresponding to 87 pc) VLBA-Gaia offset.The C-band archival image reveals a possible bipolar jet that appears to originate from the Gaia position.Although the cause of the offset is likely due to the small-scale radio jet, confirmation of the jet scenario requires spectral index analysis derived from multi-band images.

J1259+5140
J1259+5140 shows a consistent position between VLBA and Gaia, but it presents several radio components positioned towards the west.The configuration of these components aligns with the direction of the small VLBA-Gaia offsets, indicating that they may constitute a small-scale jet.Nevertheless, confirming the jet hypothesis requires spectral index analysis derived from multi-band images, and we only have an archival C-band image.

J1557+1236
J1557+1236 exhibits a significant VLBA-Gaia offset of 62.0±1.2 mas (corresponding to 485 pc).The radio component is a weak detection at around 6σ and shows slight extension.Without deeper images or additional images in other bands, the origin of the offset is unclear.

Other possible origins of VLBA-Gaia offsets
In Section 5.1, we discuss individual cases that have significant VLBA-Gaia offsets or multiple radio components assuming the origins are dual/off-nucleus quasars and small-scale radio jets.Nevertheless, under a few other possible scenarios, single normal quasars could exhibit offsets between optical and radio positions.In the following sections, we explore these scenarios.

Optical Jet
The presence of optical jets at scales of hundreds of milli-arcseconds is considered to be one of the possible mechanisms for the VLBA-Gaia offsets (Petrov & Kovalev 2017;Petrov et al. 2019).These optical jets can cause a shift in the Gaia positions along the direction of the jet, away from the nucleus.To investigate the possibility of optical jets, we examine the optical spectra from SDSS and the optical images from the Dark Energy Camera Legacy Survey (DECaLS, Dey et al. 2019), which has an image quality of approximately 1 arcsecond.We look for any indications of optical jets that could contaminate the optical spectra with a featureless power-law continuum from synchrotron radiation.Among all eight candidates, we find that they show typical spectra of type 1 broad-line quasars.Given the typical radio flux density of our sample(∼10 mJy), we do not expect any optical synchrotron to contribute (Collinge et al. 2005).Additionally, based on the DECaLS images, we observe that seven out of the eight candidates, excluding J0749+2255, which has already been identified as a dual quasar (Shen et al. 2021;Chen et al. 2023), show a point-like morphology without any signs of extended optical jets.The point-like morphology suggests that the presence of strong optical jets at scales of arcseconds is unlikely.However, faint sub-arcsecond optical jets that are unresolved in DECaLS are still possible.

Extended Host Galaxy
The astrometric solutions of Gaia detections can be affected by optical emission from host galaxies.According to a systematic study by Makarov et al. (2019), the fraction of sources with VLBA-Gaia offsets is notably higher at low redshifts (z ≲ 0.5), which can be attributed to the extended structures of host galaxies (Hwang et al. 2020).Our varstrometry-selected candidates generally have redshifts in the range of z = 1 − 3.Even for a few low-redshift targets, they still have redshifts of z ≳ 0.5.At the redshift range of most targets (z = 1 − 3), the host galaxy emission is considered negligible in the Gaia bandpass.Furthermore, the optical DECaLS images show point-like morphology for the candidates that exhibit significant VLBA-Gaia offsets or multiple radio components.Based on these observations, we conclude that the Gaia positions are primarily based on the emission from the quasars, and the host galaxies are not the cause of the offsets.

Reference coordinate system
Differences in reference frames between the VLBA and Gaia images can introduce systematic biases when comparing the two datasets.The celestial reference frame of our VLBA images is tied to the International Celestial Reference System (ICRS), and most of our phase reference calibrators have positions based on the 3rd realization of the International Celestial Reference Frame (ICRF3; Charlot et al. 2020).The Gaia celestial reference frame is officially aligned with ICRF3, with a largescale systematic uncertainty ranging from 20 to 30 µas (Gaia Collaboration et al. 2018;Lindegren et al. 2018).This systematic uncertainty of 20 -30 µas is much smaller than the VLBA-Gaia positional offsets observed in the eight candidates.Additionally, by analyzing the VLBA-Gaia positional offsets in the 21 VLBA-detected quasars, we did not observe any significant systematic error in positional offsets (see Figure 7).This lack of substantial systematic shifts supports a good alignment between the reference frames of Gaia and VLBA, indicating that the systematic uncertainties in the celestial reference frame do not account for the observed VLBA-Gaia offsets in the candidates.

Chance superposition of an unrelated radio source
The radio detection could be merely an unrelated radio source in the background or foreground.To assess the likelihood of chance superposition, we computed the probability using number statistics from the Very Large Array Sky Survey (VLASS, Lacy et al. 2020).Employing a 3σ detection threshold of ∼0.5 mJy at 8.37 GHz and assuming a typical spectral index α = 0.71 (Gordon et al. 2021), sources with a flux density greater than 1 mJy should be observable at 3 GHz.Based on VLASS epoch 1 images, there are approximately 1.7×10 6 radio sources with flux density exceeding 1 mJy across the entire 35,285 deg 2 footprint (Gordon et al. 2020(Gordon et al. , 2021)).Consequently, assuming a uniform distribution, the probability of encountering a random radio source with flux density > 1 mJy within a 0. ′′ 5 radius is ∼3×10 −6 .Therefore, it is unlikely that our radio detection is a result of a chance superposition with an unrelated radio source.

Fraction of VLBA-Gaia offset quasars
To assess whether varstrometry is successful in identifying candidate dual and offset quasars, we compare our sample with several quantitative studies focusing on normal single quasars that have both VLBA and Gaia matches.Several studies have Petrov et al. (2019) analyzed systematic offsets in the positions of 9,081 matched sources between Gaia DR2 and VLBI.They found that 9% of sources exhibit offsets statistically significant at the 99% confidence level.Makarov et al. (2019) investigated a sample of 3,413 extragalactic sources from ICRF3.After excluding potential contamination from double sources, confusion sources, extended sources, and sources with poor quality, they found that 20% of the sources had VLBA-Gaia offsets exceeding the 3σ significance level.
Our detection rate of 26±8% (6/23, assuming 1σ Poisson errors) is slightly higher but still broadly consistent (< 3σ) with the rates reported in the literature for normal single quasars (Petrov et al. 2019;Makarov et al. 2019).Figure 8 illustrates the distribution of VLBA-Gaia offsets for our sample compared to the offsets reported in Makarov et al. (2019).To ensure comparability with our sample, the quasars from the literature are selected using the criteria : G > 18.5 and 0.5 < z < 3. Our varstrometry-selected sample includes a higher number of targets with significant VLBA-Gaia offsets (the p-values of two-sample Kolmogorov-Smirnov test is 0.02), indicating the potential of the varstrometry technique to discover additional targets with such offsets (≳ 10 mas).If we only select the subset of the 21 targets that still pass varstrometry selection in Gaia DR3, the discrepancy between our sample and the control sample becomes even larger (p-values of 0.01).A larger sample with lower statistical uncertainty is needed to confirm this discrepancy definitively.Besides, the median flux density of the ICRF3 sources (∼0.13Jy at S band, Hunt et al. 2021) is one order of magnitude higher than those in our sample (∼0.02Jy at S band).Thus, the Makarov et al. (2019) sample, based on the ICRF3 sources, might contain a large portion of small-scale jets and be toward extremely radio-loud quasars.
The varstrometry technique, theoretically, is more sensitive to scales of dozens or hundreds of mas than just a few mas, given the current accuracy of Gaia.According to the varstrometry assumption (Equations 3 and 5 in Hwang et al. 2020), the typical astrometric signal for pairs with equal fluxes and 10% flux variability is approximately 5% of the pair separation.Therefore, for the typical Gaia astrometric excess noise of a few mas in our sample, the expected separations would fall within the range of dozens to hundreds of mas.Consequently, we believe that targets with offsets of around 10-100 mas could potentially be discovered using varstrometry, while those at a few mas scales are too small to be detected given the current astrometric accuracy of Gaia.(2019).We apply the criteria of G > 18.5 and 0.5 < redshift < 3 to the quasars from the literature to ensure comparability with our sample.We also show the subset of the 21 targets that still pass varstrometry selection using Gaia DR3 (in dark red).Our sample shows a notable excess in large VLBA-Gaia offsets.

CONCLUSIONS
We presented VLBA images and measurements of 23 radio-bright quasars selected by varstrometry, which have larger Gaia astrometric noise, aiming to search for candidate dual/off-nucleus quasars at pc scales.Out of the 23 quasars, 8 exhibit significant positional offsets between Gaia and VLBA or multiple radio components.
Out of the 8 candidates, 3 (J1051+2119, J1110+4817, and J1259+5140) exhibits multiple radio components.Based on the morphology and jet direction observed in single-band images, the radio emission is likely originating from small-scale jets, although confirmation through spectral index analysis using deep multiple-band observations is required.
The remaining 5 candidates with significant VLBA-Gaia offsets exhibit different characteristics.Three candidates (J0749+2255, J1044+2959, and J1110+3653) have compact radio cores with flat or inverted spectra, one candidate (J0108−0400) has extended radio emission with a steep spectrum, and one candidate (J1557+1236) has a large offset of 62.1 mas but lacks spectral index information due to only having a singleband image.Based on spectral indices and morphology, the VLBA-Gaia offset in J0108−0400 likely originates from small-scale radio jets.The offset in J0749+2255 could be attributed to blended optical emission in Gaia from a confirmed dual quasar on a larger scale of 0. ′′ 46 (Chen et al. 2023).J1044+2959 and J1110+3653 are possible candidate dual quasars at pc scales.The origin of the radio emission in J1557+1236 remains unknown due to weak detection and the lack of spectral index information.Deep multi-band follow-up imaging is necessary to confirm the properties of J1557+1236.
In addition to the radio jet and dual/off-nucleus quasar scenarios, we explored other possible causes of the VLBA-Gaia offset.Optical jets at scales of hundreds of milli-arcseconds are considered a potential cause of VLBA-Gaia offsets (Petrov & Kovalev 2017;Petrov et al. 2019).However, examination of the optical spectra and images for the eight candidates did not find any evidence of optical jets contaminating the spectra or extended jets at arcsecond scales.While strong optical jets at arcsecond scales are unlikely to be responsible for the offset, the presence of faint sub-arcsecond optical jets cannot be ruled out.Another potential factor is optical emission from host galaxies, which can affect Gaia's astrometric solutions.Most of our varstrometryselected candidates are at high redshifts (z = 1 − 3), so the host galaxy emission is negligible.Furthermore, differences in reference coordinate systems between VLBA and Gaia images can introduce systematic biases.The small systematic uncertainty of 20 − 30µas (Lindegren et al. 2018) and the absence of significant systematic shifts in the 21 VLBA-detected quasars, suggests that systematic uncertainties in the celestial reference frame is not a major contributor for the VLBA-Gaia offsets we observe.
We compare our varstrometry-selected sample and samples of normal single quasars from the literature (Petrov et al. 2019;Makarov et al. 2019) to assess whether varstrometry can select a higher fraction of quasars with VLBA-Gaia offsets.We find that the fraction of quasars with significant (> 3σ) offsets in our sample is slightly higher, but still broadly consistent with those in the literature, taking into account the uncertainty.Additionally, we examine the distribution of the VLBA-Gaia offsets and observe a greater number of large offsets in our varstrometry-selected sample, in comparison with normal single quasars, demonstrating the potential of the varstrometry technique to uncover more candidate dual or off-nucleus quasars.

Figure 1 .
Figure 1.Illustration of different physical origins that can cause VLBA-Gaia offsets.The scales of offsets shown here are usually < 0. ′′ 2, which is the Gaia's point spread function (blue ellipse).The optical position from Gaia is an unresolved weighted center (blue diamond), while the radio emission (yellow hexagon) is resolved with VLBA's mas angular resolution.(a): Dual quasar: The optical emission of two quasars is blended in Gaia resulting in a weighted optical center, while the radio emission is centered at individual quasar because of VLBA's high angular resolution.(b): Dual quasar: If only one quasar is radio-bright, the binary quasar show an optical-radio offset.(c): Off-nucleus quasar: The optical center is a weighted position of an off-nucleus quasar and a host galaxy, while the radio position is centered at the quasar's position.(d): Radio jet from a single quasar: The optical center is centered at the quasar's position, while the radio position is shifted toward the jet direction.

Figure 2 .
Figure2.Sample selection and the target properties.We selected the targets based on their peak flux densities (>15 mJy) in FIRST and astrometric excess noise >3σ at a given Gaia magnitude.Left: FIRST flux density versus redshift for the final analysis sample of 23 targets (marked in red) compared to all FIRST-detected sources in the parent quasar sample (marked in black), with a blue horizontal line indicating the 15 mJy threshold.Right: Astrometric excess noise versus Gaia magnitude for the final analysis sample of 23 targets (marked in red) compared to all radio-bright (>15 mJy) SDSS quasars (marked in balck), with a blue curve indicating the >3σ threshold of astrometric excess noise at a given Gaia magnitude, based on all SDSS quasars.

Note-
Column (1): Target name.Column (2): Receiver band.Column (3) & (4): Right Ascension and Declination of the VLBA detected source.The position corresponds to the brightest VLBA component in the presence of multiple components.Column (5) & (6): Peak and total flux density of the detected source.The 5-σ upper limit is reported if a source is undetected.Column (7): 1-σ noise level.Column (8) & (9): Major and Minor axis of the synthesized beam.Column (10): Position angle of the synthesized beam.

Figure 3 .
Figure 3. VLBA (C-band, X-band, and/or Ku-band) images of the eight (four in this figure and four in Figure 4) targets that exhibit significant offsets between VLBA and Gaia positions or/and have multiple radio components.The images are oriented with north up and east to the left.The Gaia DR2 positions are indicated by red crosses, with their lengths representing 1σ rms errors, and the red ellipses show the 3σ rms position error.The VLBA positions are marked by blue crosses, and the blue ellipses indicate the 3σ rms position error.The synthesized beam sizes are shown in the bottom-right corners of the images.The solid contours correspond to 5, 10, 20, 40, 80, 160, and 320  times the 1σ rms level.The typical 1σ rms noise level is approximately 0.1 mJy beam −1 , but for bright sources, it can reach up to 1 mJy beam −1 .

Figure 4 .
Figure 4. VLBA (C-band, X-band, and/or Ku-band) images of the eight (four in Figure 3 and four in this figure) targets that exhibit significant offsets between VLBA and Gaia positions or/and have multiple radio components.Notations are the same as those in Figure 3.

Figure 6 .
Figure 6.Radio spectra of the six targets for which we have detections at at least two frequencies.

Figure 7 .
Figure 7. Distribution of positional offsets between the Gaia DR3 optical position and the VLBA radio source positions for the 21 VLBA-detected targets.Targets with offsets greater than 3σ are highlighted in blue.The error bars represent 1σ errors considering the fitting error, the phase calibrator's position uncertainty, and the error associated with the phase-referencing technique, as described in Section 4.4.The axes of the insets are in units of mas, identical to those in the the main plot.

Figure 8 .
Figure8.Distribution of VLBA-Gaia offsets for the 21 targets that have a detection in our sample (in red) compared to those for single normal quasars (in grey) inMakarov et al. (2019).We apply the criteria of G > 18.5 and 0.5 < redshift < 3 to the quasars from the literature to ensure comparability with our sample.We also show the subset of the 21 targets that still pass varstrometry selection using Gaia DR3 (in dark red).Our sample shows a notable excess in large VLBA-Gaia offsets.

Table 1 .
Basic properties of the final sample of 23 targets.

Table 2 .
Logs of new VLBA observations for the 18 targets (17 objects without VLBA archival images plus J1044+2959).The experimental setup for X-band observations includes single polarization, 2 s integration times, and a total bandwidth of 512 MHz (16 spectral channels with a bandwidth of 32 MHz).All targets were observed in the phase-referencing mode, with nearby bright quasars used as the phase-reference calibrators, except for J0840+4439 and J1044+2959, which were self-calibrated.

Table 3 .
Source and image properties of the 23 targets having the new and archival VLBA observations