© The Authors 2026

1. Introduction

Quasars, as the most luminous class of active galactic nuclei (AGNs), lie at extreme distances and are compact in nature. They have traditionally been considered to exhibit negligible parallax and proper motions, making them ideal reference points for astrometric studies (Sovers et al. 1998; Ma et al. 1998; Gaia Collaboration 2018). However, recent high-precision astrometric observations, particularly those from the Gaia mission, have revealed that some quasars exhibit unexpected astrometric features. These anomalies include significant proper motions, parallaxes, or astrometric noise that are inconsistent with the traditional view of quasars as fixed points on the celestial sphere (Lindegren et al. 2018; Gaia Collaboration 2021). In a recent work, we identified a population of quasar candidates with anomalous astrometric properties by cross-referencing Gaia Early Data Release 3 astrometric data with spectroscopic data from the Sloan Digital Sky Survey (SDSS), which we refer to as Strange QUasar Candidates With ABnormal Astrometric Characteristics (SQUABs; Wu et al. 2022; Ji et al. 2023). Within the SQUAB population, a particularly intriguing subset comprises SDSS quasars that are associated with multiple Gaia counterparts within 1 arcsecond. We designate these as SQUAB sources with multiple Gaia detections (hereafter SMGDs). There are several compelling questions regarding the physical nature and possible astrophysical origins of SMGDs.

While chance alignments of unrelated sources at different redshifts could account for some of the SMGD cases (e.g., quasar plus a foreground star or AGN), alternative interpretations include physically associated quasar pairs or gravitationally lensed quasars (LQs). Each of these scenarios carries distinct astrophysical implications and observational challenges, highlighting the value of high-resolution and multiwavelength follow-up studies of SMGDs.

Among the possible scenarios of SMGDs, quasar pairs (or dual AGNs) are of particular importance. Dual AGN systems are kiloparsec-separated pairs of supermassive black holes (SMBHs) brought together through dynamical friction and interactions with surrounding stars and gas (e.g., Milosavljević & Merritt 2001; Mayer et al. 2007). They represent the precursors of supermassive black hole binaries (SMBHBs), where the separations shrink to parsec scales and the two black holes eventually become gravitationally bound. Identifying dual AGN systems is therefore a crucial step toward tracing the merger-driven evolution of SMBHs (Begelman et al. 1980; Volonteri et al. 2003; Khan et al. 2016).

According to hierarchical models of galaxy evolution, SMBHBs are a natural consequence of galaxy mergers. As galaxies collide and merge, their central SMBHs are expected to form a bound pair, eventually coalescing into a single, more massive black hole (Yu 2002; Kelley et al. 2017). This process is thought to be a key mechanism for the growth of SMBHs over cosmic time and is intimately linked to the coevolution of galaxies and their central black holes (Kormendy & Ho 2013).

Despite their importance, confirming the existence of dual AGNs and SMBHBs, especially from subparsec to a few-kiloparsec scales, has proven challenging (e.g., An et al. 2018). This difficulty arises from several factors. At cosmological distances, the angular separation of dual AGNs can be extremely small, often below the resolution limit of conventional optical telescopes (Yu 2002; Burke-Spolaor 2011). While double-peaked emission lines in optical spectra have been proposed to be a signature of SMBHBs, they can also arise from other phenomena, such as outflows or rotating disks (e.g., Popović 2012; Smith et al. 2010). The orbital periods of SMBHBs can be very long, making it challenging to detect periodic variations that might indicate binary motion (Graham et al. 2015). There may also exist biases in the current detection methods that make certain types of dual AGNs easier to find than others (Volonteri et al. 2009). These challenges have led to a situation where, despite dedicated efforts, the number of confirmed dual AGNs and SMBHBs remains small. The majority of known dual AGN systems are either at low redshifts or have large physical separations (> 20 kpc), with only a handful of dual quasars (DQs) with small separations found at high redshifts (Chen et al. 2020).

Very long baseline interferometry (VLBI) is increasingly recognized as a crucial tool for the study of dual AGNs and SMBHBs. Its unparalleled angular resolution on milliarcsecond scales makes it uniquely suited for directly resolving close pairs (see An et al. 2018 and the references therein). Over the years, VLBI observations have revealed a number of SMBH pairs with separations ranging from parsec to kiloparsec scales (see Liu et al. 2018 for a review). Until now, the only known parsec-scale SMBHB system had been detected via VLBI imaging (0402+379; see Rodriguez et al. 2006; Bansal et al. 2017). As for kiloparsec-scale dual AGNs, a few more cases were found serendipitously, including NGC 6240 (z = 0.025; Gallimore & Beswick 2004; Hagiwara et al. 2011; Sobolenko et al. 2022), NGC 5252 (z = 0.023; Mezcua et al. 2018; Yang et al. 2017), J1502+1115 (z = 0.39; Deane et al. 2014; Wrobel et al. 2014), J0749+2255 (z = 2.17; Chen et al. 2023a), and J1536+0441 (z = 0.389; Wrobel & Laor 2009; Bondi & Pérez-Torres 2010). A promising dual AGN candidate was also studied with combined Hubble Space Telescope and VLBI observations (J1220+1126 at z = 1.9; Glikman et al. 2023).

Many systematic VLBI-based searches for dual AGNs (or SMBHBs) failed to spot any multiple AGNs within a source. For example, Burke-Spolaor 2011 did not find any other flat-spectrum radio pairs from a sample of 3114 radio bright sources with multiband VLBI observations. Tingay & Wayth 2011 matched 11 radio AGNs from 87 promising dual AGN (or SMBH) candidates that displayed double-peaked optical emission lines. The 11 radio sources were imaged with VLBI at 1.4 GHz, revealing only 2 sources with a single, weak component (17 mJy and 1.8 mJy). Tremblay et al. 2016 made another attempt based on 1279 VLBA Imaging and Polarimetry Survey (VIPS) sources (Helmboldt et al. 2007), looking for double flat- or inverted-spectrum compact cores among a subsample of compact symmetric objects (CSOs) and CSO candidates, where no SMBHBs other than 0402+379 have been found. These results highlight that, in addition to VLBI imaging capabilities, the selection of an oriented and intrinsically bright sample is critical when searching for dual AGNs or SMBHBs.

The search for SMBHBs and kiloparsec-scale dual AGNs with different separations is essential for several reasons. They allow us to test galaxy merger models and our understanding of SMBH growth, to probe the different stages of the SMBHB merger process, and to constrain the timescales of SMBHB evolution, which have implications for the “final parsec problem” in SMBHB mergers; they also inform predictions for gravitational wave emission from these systems (Milosavljević & Merritt 2003). The Gaia-based multiple detection sample provides a bona fide and less biased parent sample that is well suited for VLBI follow-up studies. By combining the high-precision astrometric measurements from Gaia with the high-resolution imaging capabilities of VLBI, we have a unique opportunity to constrain the nature of these objects through a novel direct way.

In this paper we present the results of a pilot study that utilized this combined Gaia-VLBI approach. We focus on a subset of SMGDs with strong radio emission, hereafter referred to as radio-detected SMGDs (R-SMGDs). The target sources, namely J1433+4842, J1520+4211, and J0928+5707, were picked based on their radio flux densities determined via surveys such as the NRAO VLA Sky Survey (NVSS) and Faint Images of the Radio Sky at Twenty-centimeters (FIRST). Through VLBI observations of the three R-SMGDs, we obtained high-resolution radio images and accurate positions for comparison with Gaia detections. The combination of radio morphology and spectral analysis enabled us to search for dual AGNs and place tighter constraints on the intrinsic properties of these objects. Moreover, this approach provides a unique opportunity to clarify the origins of Gaia’s unusual astrometric measurements using the independent reference frame offered by VLBI.

We describe our sample selection based on the SMGD sample in Sect. 2. In Sect. 3 we present the data reduction and image analysis procedures. Section 4 provides our identification of the target sources and discusses the prospects of the Gaia–VLBI approach. Finally, Sect. 5 summarizes our main conclusions and outlines directions for future work. Throughout this paper, we apply the Λ cold dark matter cosmological model parameters with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Sample selection

The parent SMGD sample consists of 143 SDSS quasars. We extracted it from SQUAB sample A from our first SQUAB paper (Wu et al. 2022) by excluding any single stars. Our second SQUAB paper (Ji et al. 2023, hereafter J2023) concentrated on the examination of the SMGD sample; in that study, we conducted a series of optical identification methods to classify the sample. A straightforward version of the procedure from J2023 is outlined below:

1. Gaia detections of sources with obvious stellar features (e.g., proper motions or parallaxes greater than 5σ and/or spectral and color characteristics) were marked as quasar–star pairs.

2. After filtering out the quasar–star pairs, we searched for optical spectral features, and sources showing double-peaked emission line features were marked as DQ candidates.

3. A color-color analysis for the filtered source was conducted in parallel. Sources with colors similar to their Gaia counterparts were marked as LQ candidates.

4. Sources lacking adequate optical information (e.g., Gaia non-detections in several bands) for the above identifications were marked as "REST."

From the J2023 results, 65 quasar–star pairs, 2 DQ candidates, 56 LQs or LQ candidates, and 20 REST sources were compiled. However, using optical-only data for classification poses inherent limitations and biases. This concern arises from the following facts: (1) Confirmed sources remain scarce: only 13 LQs (out of the 56 LQ candidates) have been confirmed in the literature, and just 19 sources have high-resolution Hubble Space Telescope imaging available for more reliable verification. (2) Classification inconsistencies exist: for example, one reported DQ, J084129.77+482548.3 (Mannucci et al. 2022), is classified as a LQ candidate due to the lack of spectral confirmation. Two reported LQs, J081331.28+254503.0 (Jackson et al. 2015) and J234330.58+043557.9 (Krone-Martins et al. 2019), are classified as quasar–star pairs due to Gaia’s limitations. Two reported DQs and one reported LQ continue to carry the REST label due to the lack of necessary features.

Given this context, in addition to the search for potential dual AGNs, we consider it essential to conduct high-resolution imaging of the SMGDs using the combined Gaia-VLBI approach. This can provide evidence that we can use to robustly constrain their physical nature. To identify radio-bright candidates suitable for VLBI observations within this sample, we performed a cross-match of the SMGDs with two major radio surveys: NVSS (Condon et al. 1998) and FIRST (Becker et al. 1995).

The search method was as follows: First, we separately cross-matched the SDSS positions of the sample with NVSS and FIRST. Considering the different beam sizes and position errors of the surveys, different searching radii were used for the cross-matching: The searching radius for SDSS-NVSS was 10 arcsec and for SDSS-FIRST 3 arcsec. Second, we combined the two matched groups and generated a final sample consisting of their positional information, radio flux densities, and the separations between their Gaia positions. Based on the cross-matching procedure, we built the R-SMGD sample, consisting of 18 radio sources with detectable radio emission in at least one of the surveys (see Table A.2; the classifications from J2023 are also listed in the last column).

As no prior VLBI studies have been reported for our R-SMGD sample (except for J1520+4211, which was included in sensitivity-limited astrometric surveys), we initiated new Very Long Baseline Array (VLBA) observations (Rodriguez et al. 2006) to investigate their physical nature. The high angular resolution of the Gaia−VLBI alignment allowed us to test their classifications with unprecedented precision and search for potential DQs. To optimize the observing strategy and minimize the risk of VLBA non-detections, the brightest three sources in the table were proposed and approved for a pilot session, with follow-up observations on weaker sources planned for a later stage.

3. VLBI observation and data analysis

3.1. VLBA observation and data reduction

The VLBA observations of the three selected sources were conducted in 2023 February and March (project code: BZ94; see Table 1 for more information). The observing frequencies were 1.6 GHz (L band) and 4.9 GHz (C band). The observations were performed using the digital downconverter system. Four 128 MHz baseband channels on both left and right circular polarization and 2-bit sampling were utilized as the basic observing configuration, resulting in a total data rate of 4 Gbit s⁻¹ for the high-sensitivity VLBI observation. Given that the target sources had not before been observed by VLBI, their high-resolution radio cores could be dimmer and more compact than their large-scale radio sources. Thus, the phase referencing technique was used for all three target sources. The phase calibrators were searched for using the astrogeo Calibrator Search Engine¹. The phase calibrators used for the three sources were all standard calibrators (with a “C” flag) and were within 2 degrees of the target source position. This guaranteed the accuracy and imaging quality of the phase that applied to the target source. We used a 5-minute phase-referencing cycle time (4 min on target and 1 min on calibrator) for each target source in both observing bands. Since the two Gaia positions of each source are not very distant (∼1 arcsec), we only used the first Gaia position (J2000) as the target position, i.e., $\mathrm{RA}={14}^{\mathrm{h}}{33}^{\mathrm{m}}33\stackrel{\text{s}}{.}0316$, $\mathrm{Dec}=+48°42\prime 27\stackrel{″}{.}7752$ for J1433+4842; $\mathrm{RA}={15}^{\mathrm{h}}{20}^{\mathrm{m}}39\stackrel{\text{s}}{.}7218$, $\mathrm{Dec}=+42°11\prime 11\stackrel{″}{.}4524$ for J1520+4211; and $\mathrm{RA}={09}^{\mathrm{h}}{28}^{\mathrm{m}}53\stackrel{\text{s}}{.}5310$, $\mathrm{Dec}=57°07\prime 35\stackrel{″}{.}9477$ for J0928+5707.

After the observations, the raw data from each antenna were transferred to the VLBA correlation center in Socorro (New Mexico, USA) for correlation at the Distributed FX (DiFX) correlator (Deller et al. 2011) with a 2s integration time. For the second Gaia position of each source, we also requested the additional correlation of the same data using the second Gaia position as the target phase center, which helps reduce the smearing effect and increases the imaging quality for both Gaia positions. We then downloaded the correlated data and conducted additional data calibration and imaging using the VLBI pipeline at the China Square Kilometre Array (SKA) Regional Centre (An et al. 2019, 2022a). The pipeline takes advantage of the US National Radio Astronomy Observatory (NRAO) AIPS (Astronomical Image Processing System) software (Greisen 2003) and the ParselTongue tool to translate AIPS tasks into Python scripts (Kettenis et al. 2006) and automatically calibrates the data with minor human intervention. The calibration steps follow the standard calibration procedure of the AIPS Cookbook² and can be found in Zhang et al. (2022) and Baoqiang et al. (2023). After calibration, the calibrated visibility data were exported by the pipeline and loaded using the Caltech DIFMAP software package (Shepherd 1997) for imaging. We made use of the hybrid mapping procedure, which includes several loops of CLEAN and self-calibration steps in DIFMAP. After that, iterations of phase and amplitude self-calibration with decreasing integration time (divided by 2 in each iteration, down to 2 min) were conducted to generate the final images. Note that for source J1433+4842, only phase self-calibration was performed, to avoid spurious detections due to the weakness of the source (e.g., Martí-Vidal & Marcaide 2008). We performed this calibration procedure for both the first and second Gaia positions for each source. This was accomplished by separately processing the data whose correlation phase center corresponded to the respective Gaia position, to eliminate the bandwidth smearing effect, which would otherwise blur the detection of sources far from the imaging phase center. For source J1433+4842, we also did an extra shift-and-CLEAN procedure for a newly found extending component, which lies ∼280 mas from its first Gaia position (see Sect. 3.2 for more details). This extra step was done in AIPS with the task UVFIX to shift the imaging phase center to the component peak before any channel or time averaging of the data was conducted. Then the shifted visibility data were exported again for another imaging round in DIFMAP in order to get a better image of the extending component. The image pixel cell sizes were adopted following the recommendations of DIFMAP during the CLEAN process: 1.8 × 1.8 mas at 1.6 GHz and 0.5 × 0.5 mas at 5 GHz.

Following the imaging procedure, we performed model-fitting in DIFMAP on the self-calibrated visibility data to parameterize the VLBI components. During the fitting process, we initially adopted a circular Gaussian brightness distribution model centered on the component peaks of each detected source feature. If the residual peak remained beam-like in shape, an additional circular model was applied. In cases where the residual appeared too compact or exhibited an irregular structure, a point-source (delta-function) model was applied instead. The fitting process was iterated until no significant residuals above the 5σ level remained within the CLEANed area. Parameter descriptions from the fitting results are presented in Sect. 3.2. The uncertainties in the fitted component parameters were primarily estimated following the method described by Fomalont (1999). For flux density errors, we accounted for an additional 5% uncertainty introduced during the calibration process. In terms of VLBI positional uncertainties, we adopted the widely used $\frac{{\theta }_{\mathrm{beam}}}{2\mathrm{S}/\mathrm{N}}$ formula (where θ_beam is the full width at half maximum size of the beam and S/N is the signal-to-noise ratio of the component) as the measurement error for VLBI observations (e.g., Yang et al. 2010; An et al. 2020). Moreover, since we used phase-referencing, positional errors can also arise from the calibrator’s position. As a result, the uncertainties in the phase calibrator’s position are considered the dominant source of systematic error. The final positional uncertainties are thus the quadratic means of these two errors.

3.2. Image analysis

Our VLBI observations of the three target sources (J1433+4842, J1520+4211, and J0928+5707) yielded high-resolution images at both 1.6 GHz and 4.9 GHz. For each source, we present the CLEAN wide-field image at 1.6 GHz, covering both Gaia positions, along with the CLEAN zoomed-in images of individual components at 1.6 and 5 GHz, with astrometric positions clearly marked. Appendix B displays the large-scale images, where square boxes highlight the Gaia positions or components. The corresponding high-resolution zoom-ins are shown in Appendix C. In the various figures, we label the first and second Gaia detections as GAIA1 and GAIA2, respectively. The CLEAN image results are listed in Table 2.

For J1433+4842, from Figs. B.1 and C.1, we detected a compact radio component close to the GAIA1 position, while no significant radio emission is detected at its GAIA2 position. We also discovered an extended jet-like structure (hereafter EXT1) between the two Gaia positions, at both 1.6 and 5 GHz. The peak intensity of EXT1 is located at $\mathrm{RA}={14}^{\mathrm{h}}{33}^{\mathrm{m}}32\stackrel{\text{s}}{.}9876$ and $\mathrm{Dec}=48°42\prime 27\stackrel{″}{.}5872$, approximately 280 mas southwest of GAIA1. For sources J1520+4211 and J0928+5707, the imaging results are similar: only GAIA1 is detected, with no significant radio emission at the GAIA2 positions. At both 1.6 and 5 GHz, GAIA1 in J1520+4211 exhibits a core-jet morphology, while J0928+5707 appears as a compact source.

3.2.1. Flux density measurements

From the imaging and model-fitting procedures described above (see Sect. 3.1), we obtained high-resolution positions and flux density measurements of each source component. For each source, the VLBI core component associated with GAIA1 is well described by a single circular Gaussian model, which was used to determine the VLBI position (Table A.1) and to calculate the brightness temperature (Table 3).

In Table A.1, for each marked feature in Col. 2, Cols. 3–5 report the optical positions and errors from Gaia detection. Columns 6–11 report the imaging and model fitting results from our VLBI observation; Col. 7 represents the peak flux densities of each feature. For non-detections, the 5σ upper limits are presented. Columns 8 and 9 are the positional offsets between the Gaia detection and the VLBI core component position along RA and Dec, separately. Column 10 lists the local post-CLEAN rms noise derived from the pixels within the zoomed-in area (shown in Appendix C). Column 11 lists the total flux density of each radio feature. This was obtained by summing the fitted core component and the remaining model components associated with the feature derived during the model fitting process.

Table 3 presents the core parameters associated with GAIA1. Column 3 denotes the integrated flux densities of the fitted core components. The full width at half maximum size of the core is listed in Col. 4, where the size upper limits for unresolved components are derived based on Kovalev et al. 2005 (i.e., if the fitted component size is smaller than this value, this value is used as the upper limit for the component size). The brightness temperatures of the core component in Col. 5 were calculated as (e.g., Condon et al. 1982)

$$\begin{array}{c}\hfill T_{\mathrm{B}}=1.22\times {10}^{12}(1+z)\frac{S_{\nu }}{{\theta }_{\mathrm{comp}}^2{\nu }^2}\mathrm{K},\end{array}$$(1)

where S_ν is the component flux density (in Jy), ν the observing frequency of the image (in GHz), and θ_comp is the full width at half maximum of the circular Gaussian component diameter (in mas). If θ_comp is an upper limit, the calculated T_B was used as a lower limit. The fact that the T_B values of all the core components exceed 10⁸ K indicates that the GAIA1 features in our sample are predominantly powered by AGN activity.

3.2.2. Spectral index map

To gain further insights into the physical nature of these sources, we generated spectral index maps between 1.6 GHz (L band) and 4.9 GHz (C band). To avoid ambiguity, we defined the radio spectral index (α) according to the power-law relation S_ν ∝ ν^α, where S_ν is the flux density at frequency ν. In this convention, a negative α value indicates that the flux density decreases with increasing frequency.

The phase-referencing observations provide high-quality positions of the radio source at the milliarcsecond level, making the alignment of maps from different frequencies convenient and accurate. The positions in Table A.1 were used as image centers for the alignments. To account for differences in beam size and pixel scale between images at different frequencies, we adopted the mapping parameters (i.e., restoring beam size, pixel size, and map size) from the lower-frequency L-band images as a reference. The higher-frequency C-band images were then reconstructed using the same beam and pixel configurations as the referenced ones to ensure consistency in the mapping of spectral index values. Before computing the spectral index maps, we applied a masking procedure to exclude pixels with intensities below 3σ of the local image rms noise. This step helped eliminate spurious α values arising from noise-dominated regions.

Following the alignment and mapping procedures described above, spectral index maps for each radio feature were successfully generated, as shown in Fig. 1. To assess the reliability of these maps, we also present the corresponding spectral index error maps (see Fig. D.1), computed through numerical error propagation following standard methods (e.g., Rashid et al. 2024).

Fig. 1. Spectral index maps of the target sources. Contours denote the low-frequency (1.6 GHz) intensity map, with the corresponding synthesized beam shown in the bottom-left corners. The lowest contour corresponds to 4σ, and subsequent levels increase by factors of two. The color scale represents the spectral index distribution between 1.6 and 4.9-GHz (see Sect. 3.2.2).

Figure 1 (top panels) displays the spectral index maps for both the GAIA1 and EXT1 of J1433+4842. The radio core in GAIA1 shows a flat-spectrum feature (α ≈ 0.0). The spectral index map of EXT1 reveals a predominantly steep-spectrum jet structure with α ≤ −0.5 across most of its extent; it becomes steeper toward the western end. The spectral index map for J1520+4211 (Fig. 1, bottom-left panel) shows a typical core-jet structure. The core region exhibits a flat spectrum with α ≈ −0.2, typical of synchrotron self-absorbed emission from the jet base. Toward the southwest, the spectrum steepens progressively, reaching α ≈ −1.0 at the farthest extent of the detected emission. For source J0928+5707 (Fig. 1, bottom-right panel), a flat-spectrum core with α ≤ −0.5 is revealed.

4. Discussion

4.1. Source classification from VLBI identification

Based on our VLBI imaging results for the three target sources (J1433+4842, J1520+4211, and J0928+5707), the observed radio cores are in close alignment with their first Gaia detections (GAIA1). The compactness, high brightness temperature (e.g., ≥10⁸ K), and flat-spectrum characteristics of the radio cores strongly indicate a AGN nature for their GAIA1 features. Although no significant radio emission at their GAIA2 positions was detected, our VLBA observations provide a 5σ upper flux density limit of approximately 0.1 mJy/b. Based on the GaiaG-band detections, the optical flux ratios, F_GAIA1, opt/F_GAIA2, opt, for J1433+4842, J1520+4211, and J0928+5707 are 5.6, 1.0, and 0.3, respectively. On the basis of strong lensing laws, we would expect comparable flux ratios in the radio band. This corresponds to predicted radio peak intensities at 1.6 GHz for the GAIA2 counterparts of ∼0.2, 31.4, and 11.6 mJy beam⁻¹, respectively. These estimated values exceed our derived flux density thresholds for GAIA2, thereby allowing us to rule out the LQ scenario. Given that no additional compact radio components are observable aside from their GAIA1 features, the most plausible explanation for their optical GAIA2 detections is a foreground star.

4.2. Information on individual sources

4.2.1. J1433+4842

Wang et al. 2023 first studied J1433+4842 using Enhanced Multi Element Radio-Linked Interferometer Network (eMERLIN) observations, focusing on its astrometric jitter detected by Gaia. Their study, which used a selection method based on Gaia’s astrometric excess noise, detected a single compact radio source close to the Gaia position with a small Gaia-radio offset. They suggested that the observed astrometric jitter could be due to contamination from nearby sources rather than being indicative of a dual AGN system. Schwartzman et al. 2024 (hereafter S24) included J1433+4842 in their VaDAR (Varstrometry for Dual AGN using Radio interferometry) pilot study. Based on the jet morphologies and spectral index in S-band (2–4 GHz) and X-band (8–12 GHz) Very Large Array (VLA) observations, they identified J1433+4842 as a single jetted AGN rather than a dual AGN candidate.

Our high-resolution VLBI results for J1433+4842 are consistent with previous identifications and additionally reveal a distant jet feature, EXT1, that originated from the radio cores in GAIA1. By superimposing our radio features onto the VLA images from S24, we observe a strong concordance between the VLA images and our VLBI features of GAIA1 and EXT1 (see Fig. D.2).

4.2.2. J1520+4211

Our VLBI images of J1520+4211 reveal a core-jet morphology, which is characteristic of radio-loud AGNs. At the GAIA1 position, the radio core exhibits peak intensities of 32.50 ± 1.63 mJy/beam at 1.6 GHz and 25.96 ± 1.30 mJy/beam at 4.9 GHz. The total flux densities are measured to be 42.37 ± 2.12 mJy at 1.6 GHz and 28.90 ± 1.45 mJy at 4.9 GHz. At the GAIA2 position, no compact radio emission was detected. We derive 5σ upper limits of 0.11 mJy/beam at 1.6 GHz and 0.07 mJy/beam at 4.9 GHz.

4.2.3. J0928+5707

We identified a single, compact radio core with flux densities of 4.09 ± 0.21 mJy at 1.6 GHz and 1.03 ± 0.06 mJy at 4.9 GHz for J0928+5707 at the GAIA1 position. We established 5σ upper limits of 0.12 mJy/beam at 1.6 GHz and 0.09 mJy/beam at 4.9 GHz for its GAIA2 location.

4.3. Prospects of the Gaia-VLBI method

Xu et al. 2024 (hereafter Xu2024) performed a VLBI imaging study of four dual-AGN systems that had been confirmed by Fu et al. 2015a, both spatially and spectroscopically. The four dual AGNs show well-separated, spatially correlated double morphology in their VLA A-configuration 5 GHz observations (see Fu et al. 2015b), with each component showing steep spectra (α ≤ −0.5). In the Xu2024 results, none of them show detection (above 5σ) on both components down to a imaging sensitivity of ∼0.025 mJy/beam. Two of them are revealed to have a single VLBI component at one of the component centers, with flux densities of 0.2 and 0.4 mJy, respectively. The VLA-detected radio components suggest that radio emission from these dual AGNs are dominated by less energetic jets, AGN-driven outflows, corona and star formation.

Based on the “standard routine” of directly looking for dual AGNs or SMBHBs in the literature (see, e.g., Rubinur et al. 2018; Das et al. 2018, and our Sect. 1), discoveries of dual AGNs and SMBHBs heavily rely on optical imaging and spectroscopic data to highlight the spatial correlation between two objects. Searching for dual AGNs or binary systems purely through radio observations is still challenging, partly due to the current fact that most quasars in pair systems tend to be radio-quiet or lobe-dominated galaxies with weak radio cores.

Very long baseline interferometry plays a crucial role in identifying active radio cores in AGN systems and is indispensable for searching for dual AGNs at higher redshifts. For instance, at redshifts (z) > 0.5, a projected physical separation of 1 kpc corresponds to an angular separation of less than 0.2 arcsec, a scale at which only VLBI can directly resolve the two components in radio observations. In recent years, different approaches utilizing Gaia data have proven to be another prominent way of looking for dual AGNs, especially at high redshifts (i.e., z > 1); they are called the “varstrometry” method (e.g., Hwang et al. 2020; Chen et al. 2023b) and the “Gaia multi-peak” (GMP) method (see Mannucci et al. 2022, 2023). In the varstrometry method, sub-kiloparsec AGN pairs are searched for by looking at the variations in the photocenter (or “jittering”) from Gaia astrometric parameters (e.g., astrometric_excess_noise) of unresolved sources; this method has been successfully used to discover a few candidates and confirmed sources (e.g., Shen et al. 2021; Wang et al. 2023). The GMP method takes advantage of the ipd_frac_multi_peak parameter from the Gaia sources. This method can detect potentially multiple sources within a Gaia point spread function and hopefully help identify DQs or LQs at source separations in the range 0.1" to 0.7″ (see Mannucci et al. 2022).

Our approach of combining VLBI and high-resolution Gaia observations complements existing methods for searching for dual AGNs, particularly in excluding contaminants among dual AGN candidates with ambiguous optical features. Thanks to the milliarcsecond-scale angular resolution achievable through this method, dual AGNs can be directly resolved with projected separations down to several parsecs at moderate redshifts (assuming a resolution of ∼1 mas for redshifts above 0.5). These results motivate our ongoing VLBI follow-up studies of the remaining fainter R-SMGDs (see Table A.2).

5. Conclusion

In this study we conducted a dual-frequency VLBI imaging study of the three radio-brightest sources in our compiled R-SMGD sample: J1433+4842, J1520+4211, and J0928+5707. We have provided the first high-resolution images of these sources. The present work serves as a pilot investigation to evaluate the effectiveness of the SQUAB classification from J2023 and to search for potential dual AGN systems through a combined Gaia-VLBI approach.

For each target, we obtained VLBI images at 1.6 GHz and 4.9 GHz for both Gaia components. We parameterized the source radio components to estimate key physical properties (e.g., brightness temperature). We produced spectral index maps to help constrain the nature of the radio emission and shed light on the origin of their multiple Gaia detections. A summary of our observational results is provided in Table 4.

All three sources show clear evidence of AGN activity coinciding with their primary Gaia detections (GAIA1), while no significant radio emission was detected at the positions of their secondary Gaia components (GAIA2). In the case of J1433+4842, we also detected extended radio emission spatially between the two Gaia positions, which is interpreted as jet structure originating from the AGN at GAIA1. Based on these results, we confirm that all three sources are best explained as quasar–star pairs, in agreement with previous studies.

Our pilot study on the R-SMGDs proved the validity of the SQUAB criteria and demonstrates that the incorporation of Gaia and VLBI can be a promising method for searching for dual AGNs and for ruling out contaminants among dual AGN candidates. With the sensitivity achieved in our pilot observations, we are able to detect radio emission down to ∼0.1 mJy beam⁻¹, recovering flux densities as low as 10% of the VLA total flux densities for fainter sources. Encouraged by these results, our follow-up studies (Zhang et al. in prep.) on the remaining R-SMGDs – including several DQ and LQ candidates identified by J2023 – are expected to expand the census of dual AGN candidates. Future studies that use more sensitive optical facilities (e.g., JWST) and upcoming Gaia data releases will likely yield additional valuable candidates for expanding Gaia-VLBI investigations into DQ searches.

Acknowledgments

We are deeply grateful to the anonymous referee for the profound feedback, which has greatly enhanced both the scientific merit and the clarity of this work. This work is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDA0350205. This work used resources of China SKA Regional Center prototype funded by the Ministry of Science and Technology of the People’s Republic of China and the Chinese Academy of Sciences (An et al. 2019, 2022b). YZ is supported by the National SKA Program of China (grant no. 2022SKA0120102), the China Scholarship Council (No. 202104910165), the Shanghai Sailing Program under grant number 22YF1456100. YZ thanks Ivy Wong at CSIRO Space & Astronomy, Australia, for her warm hospitality and helpful comments. YZ thanks the valuable suggestions and comments from Sándor Frey and Krisztina Éva Gabányi in Konkoly Observatory, Hungary. TA acknowledge the support from the Xinjiang Tianchi Talent Program. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). This work made use of public data from the data archive of the National Radio Astronomy Observatory (NRAO; https://data.nrao.edu/). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under a cooperative agreement by Associated Universities, Inc.

References

Appendix A: Additional tables

Appendix B: The VLBI big images of the target sources at 1.6 GHz

Fig. B.1. Naturally weighted wide view VLBI image of J1433+4842 at 1.63 GHz. It includes the first Gaia position (GAIA1) the second Gaia position (GAIA2) and the extended emission that lies between the two positions (EXT1). To note that the shape of EXT1 in the picture is not exact the same with the zooming-in images showing in Fig. C.1 due to the bandwidth smearing effect caused by the large distance between the source and the image phase center (I used the GAIA1 position as the big image’s phase center). The lowest contours represent 4σ (σ the post-CLEAN background noise) and contour levels increase by a factor of 2.

Fig. B.2. Naturally weighted wide view VLBI image of J1520+4211 at 1.63 GHz showing both GAIA1 and GAIA2 positions. The lowest contours represent 4σ and contour levels increase by a factor of 2. The single-source zooming-in image of each position are shown in Fig. C.2.

Fig. B.3. Naturally weighted wide view VLBI image of J0928+5707 at 1.63 GHz showing both GAIA1 and GAIA2 positions. The lowest contours represent 4σ and contour levels increase by a factor of 2. The single-source zooming-in image of each position are shown in Fig. C.3.

Appendix C: VLBI zoomed-in images of the target sources at 1.6 GHz and 4.9 GHz

Fig. C.1. Naturally weighted VLBI single-source images of each target position for J1433+4842 at 1.7 and 4.9 GHz. The white crosses in the GAIA1 and GAIA2 images mark the corresponding Gaia positions in each source position. For each CLEANed image, the synthesized beam shape are shown at the bottom left corner. The lowest contours represent 4σ and contour levels increase by a factor of 2. For the GAIA2 images, no significant signal can be detected above 5σ. For this source in the image, 1 mas represents 8.413 pc.

Fig. C.2. Naturally weighted single-source VLBI images of J1520+4211 at 1.7 and 4.9 GHz. The imaging notes are the same as in Fig. C.1. The scale factor of this source is 6.038 pc/mas

Fig. C.3. Naturally weighted single-source VLBI images of J0928+5707 at 1.7 and 4.9 GHz. The imaging notes are the same as in Fig. C.1. 1 mas represents 8.467 in the above images.

Appendix D: Additional figures

Fig. D.1. Spectral index error distributions corresponding to the spectral index maps for Fig. 1.

Fig. D.2. Naturally weighted large scale radio structure for J1433+4842. The upper two panels are derived from the NRAO FIRST (1.4 GHz) and VLASS (3 GHz) survey. The bottom two are from the VLA observations from S24 at 3 and 10 GHz. The stars marked the Gaia positions. The lowest contours are 4 times of the background noise and contours increase by a factor of two. The beams are shown at the bottom left corner of each image.

All Tables

All Figures

	Fig. 1. Spectral index maps of the target sources. Contours denote the low-frequency (1.6 GHz) intensity map, with the corresponding synthesized beam shown in the bottom-left corners. The lowest contour corresponds to 4σ, and subsequent levels increase by factors of two. The color scale represents the spectral index distribution between 1.6 and 4.9-GHz (see Sect. 3.2.2).
In the text

Fig. B.1. Naturally weighted wide view VLBI image of J1433+4842 at 1.63 GHz. It includes the first Gaia position (GAIA1) the second Gaia position (GAIA2) and the extended emission that lies between the two positions (EXT1). To note that the shape of EXT1 in the picture is not exact the same with the zooming-in images showing in Fig. C.1 due to the bandwidth smearing effect caused by the large distance between the source and the image phase center (I used the GAIA1 position as the big image’s phase center). The lowest contours represent 4σ (σ the post-CLEAN background noise) and contour levels increase by a factor of 2.

In the text

	Fig. B.2. Naturally weighted wide view VLBI image of J1520+4211 at 1.63 GHz showing both GAIA1 and GAIA2 positions. The lowest contours represent 4σ and contour levels increase by a factor of 2. The single-source zooming-in image of each position are shown in Fig. C.2.
In the text

	Fig. B.3. Naturally weighted wide view VLBI image of J0928+5707 at 1.63 GHz showing both GAIA1 and GAIA2 positions. The lowest contours represent 4σ and contour levels increase by a factor of 2. The single-source zooming-in image of each position are shown in Fig. C.3.
In the text

Fig. C.1. Naturally weighted VLBI single-source images of each target position for J1433+4842 at 1.7 and 4.9 GHz. The white crosses in the GAIA1 and GAIA2 images mark the corresponding Gaia positions in each source position. For each CLEANed image, the synthesized beam shape are shown at the bottom left corner. The lowest contours represent 4σ and contour levels increase by a factor of 2. For the GAIA2 images, no significant signal can be detected above 5σ. For this source in the image, 1 mas represents 8.413 pc.

In the text

	Fig. C.2. Naturally weighted single-source VLBI images of J1520+4211 at 1.7 and 4.9 GHz. The imaging notes are the same as in Fig. C.1. The scale factor of this source is 6.038 pc/mas
In the text

	Fig. C.3. Naturally weighted single-source VLBI images of J0928+5707 at 1.7 and 4.9 GHz. The imaging notes are the same as in Fig. C.1. 1 mas represents 8.467 in the above images.
In the text

	Fig. D.1. Spectral index error distributions corresponding to the spectral index maps for Fig. 1.
In the text

Fig. D.2. Naturally weighted large scale radio structure for J1433+4842. The upper two panels are derived from the NRAO FIRST (1.4 GHz) and VLASS (3 GHz) survey. The bottom two are from the VLA observations from S24 at 3 and 10 GHz. The stars marked the Gaia positions. The lowest contours are 4 times of the background noise and contours increase by a factor of two. The beams are shown at the bottom left corner of each image.

In the text