© The Authors 2025

1. Introduction

The nature of dark matter (DM) remains elusive. The current and most widely accepted cosmological paradigm, the ΛCDM model (see, e.g. Frenk & White 2012, for a review), has been remarkably successful at explaining a great number of observational properties of the Universe, especially observations related to the large-scale structure and the evolution of the Universe (Vogelsberger et al. 2014). However, several discrepancies have appeared on smaller scales. These include the problem of missing satellites (Klypin et al. 1999), the diversity problem of galaxy rotation curves (see, e.g. Oman et al. 2015), and the cusp-core problem (see, e.g. Oh et al. 2011). (See also Bullock & Boylan-Kolchin 2017 for a review of these topics.)

A variety of alternative DM models have been proposed to alleviate these problems, including warm DM (WDM; Viel et al. 2005), self-interacting DM (SIDM; Spergel & Steinhardt 2000), fuzzy DM (Kulkarni & Ostriker 2022), and primordial black holes (PBHs; Clesse & García-Bellido 2018). All the proposed DM models make diverse predictions on the abundance and density profiles of DM halos on different mass scales, with CDM predicting a larger number of sub-galactic halos than WDM, for example. The mass range between 10⁶ and 10⁹ M_⊙ turns out to be an especially crucial discriminating factor between the different models (see, e.g. Lovell 2020). To include both DM sub-halos and a possible population of (free-floating) supermassive black holes (SMBHs), which might be primordial, we refer to these objects as supermassive compact objects (SMCOs) from hereon. In particular, PBHs have recently gained increased attention with the latest results from the James Webb Space Telescope (JWST; see, e.g. Maiolino et al. 2024), but attention was also garnered earlier with the detection of high mass black hole (BH) mergers by the LIGO/Virgo collaboration, such as GW190521 (Abbott et al. 2020). Because DM halos on the mass scales described above are not expected to form galaxies (Benitez-Llambay & Frenk 2020), they remain largely unprobed.

A unique way to probe SMCOs on these mass scales is through gravitational lensing. Multiple studies have been investigating the effect that different DM sub-halos have on the gravitational lensing potential (see, e.g. Gilman et al. 2023) and thus also on the properties of the multiple-lensed images for the case of strong lensing by galaxy clusters (see, e.g. Vegetti et al. 2023, for a recent review). In addition, the lensing effect has been used in fast radio bursts (FRBs) to study the properties of DM on sub-solar mass scales (see, e.g. Leung et al. 2022) and in searches of lensing signatures in light curves of γ-ray bursts (GRBs; e.g. Kalantari et al. 2021).

When considering any DM halo within the mass scales given above along the line of sight, part of the halo should be able to act as a gravitational lens for a background source if that part of the halo is sufficiently dense (Press & Gunn 1973). That means it should have a surface mass density above the critical value of Σ_crit  ≃  1.6  ×  10³ M_⊙ pc⁻², assuming lens and background sources are at cosmological distances. SMBHs will naturally fulfil this criterion. For DM halos, it depends crucially on the concentration-mass relation (Wang et al. 2020), which dictates what fraction of the actual halo mass can effectively act as a gravitational lens (see also Loudas et al. 2022).

For lensing in the mass range between 10⁶ and 10⁹ M_⊙, the expected angular separation between lensed images should be of the order of ∼2  ×  10⁻⁶(M_SMCO/M_⊙)^0.5 arcseconds if the source and lens are at cosmological distances (Wilkinson et al. 2001). For our mass range of interest, that translates to ∼1 − 100 mas. These spatial scales can be directly probed by very long baseline interferometry (VLBI) if the lensed background source is a radio-loud, compact active galactic nucleus (AGN). We refer to such systems as milli-lenses from hereon. A previous search in a sample of 300 AGN found no such objects (Wilkinson et al. 2001), and thus the cosmological number density of uniformly distributed SMCOs in the mass range 10⁶ to 10⁸ M_⊙ could be constrained to make up no more than ∼1% of the closure density Ω_total = 1.

More recently, Spingola et al. (2019a) conducted a search for milli-lenses with image separations > 100 mas at 1.4 GHz, corresponding to lensing mass scales > 3  ×  10⁹ M_⊙. Their search resulted in finding one previously unidentified lens candidate. With this, they constrained the lensing rate to 1 : (318  ±  225). Other searches for milli-lenses have also been conducted, but mostly in singular cases. Peirson et al. (2022) found evidence for milli-lensing in the long-term light curves of a blazar, most likely caused by an intervening molecular cloud or an undetected dwarf galaxy with a massive BH.

In a previous paper (Casadio et al. 2021, from now on Paper I), we described the framework in which we aim to constrain the number density of milli-lenses in a large sample of VLBI images of radio sources. With this, Loudas et al. (2022) investigated how different DM models can be constrained based on the number of observed milli-lenses. We found that with a large enough sample with known redshifts, one can have significant constraining power over some viable DM models. This is one of the goals of the Search for MIlli-LEnses (SMILE) project.

The Cosmic Lens All-Sky Survey (CLASS; Myers et al. 2003), conducted with the Very Large Array (VLA), provides the flux-limited sample as the basis for the SMILE project. The survey includes all sources within [0; +75°] declination and galactic latitude |b|  ≥  10°, with S_5 GHz  ≥  30 mJy and spectral index α ≤ 0.5 between the flux density of the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) at 1.4 GHz and GB6 (Gregory et al. 1996) at 5 GHz, and it contains 11,685 sources. The CLASS catalogue is an ideal starting point for the SMILE project because of the flat spectrum selection criterion, which maximises chances of detection of radio nuclei, and the available VLA 8 GHz flux density measurements, which serve as a hint for the expected Very Long Baseline Array (VLBA) 4–8 GHz flux density. From the CLASS sample, we chose a subsample of 4968 sources with S_8 GHz  ≥  50 mJy to increase the probability of detection with VLBI. For the project, we will mostly use archival data from the NRAO archive and combine them with new observations conducted with the VLBA (e.g. Project ID: BP267; PI: Felix Pötzl).

In a pilot project described in Paper I, we searched for milli-lens candidates in a different sample of 13 828 sources taken from the Astrogeo VLBI image database¹ (Petrov & Kovalev 2025). Serving as a pathfinder for the SMILE project, in this paper we present a thorough investigation of the nature of the 40 best candidates presented in Paper I. We combine the available archival data with new observations obtained with the European VLBI network (EVN), which we describe in Sect. 2. In the following, we assume a standard cosmology according to Planck Collaboration VI (2020), with Ω_m = 0.315, Ω_Λ = 0.685, and H₀ = 67.4 km/s/Mpc.

2. Data

The 40 lens candidates presented here have been chosen from images of 13 828 unique sources in the Astrogeo database at 2.3 (S-band), 5 (C-band) and 8 GHz (X-band). The details of the selection process are described in Paper I, which also lists the coordinates of each source. The basic criteria for the selection were i) a minimum of two compact components were found in the Astrogeo images and ii) the surface brightness ratio (SBR) of the two components, based on Gaussian model-fitting, should be < 7. We chose this more conservative threshold in Paper I to account for the very heterogeneous Astrogeo sample with partly poor (u, v) coverage and short integration times. Based on our improved error analysis (see Sect. 4.1), we revise this threshold in Sect. 4.2. Of the 40 candidates, 38 had available X-band data, and 35 of those were in the declination limit for the EVN. For those sources, we then obtained follow-up observations with the EVN at 5 and 22 GHz (K-band) in phase-referencing (PR) mode, to help detect the possibly very weak sources and to retain the relative positions of the components between the two frequencies with respect to a compact calibrator. The sources were observed in five sessions, as summarised in Table 1. At C-band, calibrator and target sources were observed in cycles with scans of 90 s and 210 s length, respectively. At K-band, scans of 120 s for calibrators and 90 s for targets were chosen to capture the faster atmospheric variability at higher frequencies. Typical on-source integration times were 20-30 minutes at C-band and 60-70 minutes at K-band for the science targets. Two representative (u, v) coverages are shown in Fig. 1. Together with our new data, for 16 sources, we have data for more than two epochs; for 19 sources, we have data for exactly two epochs at C- or X-band.

Fig. 1. Two examples of (u, v) coverages for the EVN observations at C-band. Left: J2209+6442 (high declination target). Right: J2347-1856 (low declination target).

3. Calibration and imaging

The new EVN data were calibrated with the AIPS package (Greisen et al. 2003) in a standard manner. First, the a priori amplitude calibration and the parallactic angle correction has been applied by copying the respective calibration tables from the EVN pipeline, as is recommended by the EVN. Only a few stations provide opacity-corrected data or weather information (important at higher frequencies) for the EVN, so no further corrections have been applied. Uncertainties in the amplitude calibration are included in the gain factor in Eq. (4). In the next step, corrections for the ionospheric dispersive delay were applied to the data at C-band, while it is negligible at K-band. Following that, we removed the instrumental delays by doing a fringe-fit with the task FRING in AIPS on a short scan of a bright calibrator. Then, we applied a global fringe-fit with FRING on the calibrators, as well as a complex bandpass correction. The calibrators were then subsequently imaged in difmap (Shepherd 1997), with careful self-calibration first in phase only, and then also in amplitudes. We chose a natural weighting scheme for the mostly compact sources for increased sensitivity (uvweight 0,-1 in difmap). The obtained images were then re-read into AIPS to obtain improved self-calibration solutions in the task CALIB. Those solutions were then applied to the target sources, which were subsequently imaged as well. This procedure removes residual gain errors as well as the effects of residual source structure in the calibrators. Since the application of self-calibration on a noisy data set can introduce false flux and structure and even create significant flux from pure noise (see, e.g. Martí-Vidal & Marcaide 2008), we did not attempt self-calibration in images with S/N  ≲  5 in the dirty image. An example of final images at both 4.9 GHz and 22.2 GHz is shown in Fig. 2 for the source J2209+6442. In total, we detected 13 out of 35 sources that were observed with the EVN at 22.2 GHz (see Table 5). Of the 13 detections, seven have both of the components detected that were originally identified as possible lensed images at lower frequencies.

Fig. 2. J2209+6442, a source that has both components detected at 22 GHz and has been rejected as a milli-lens but is retained as a CSO candidate. Left: EVN 4.9 GHz data. The contours start at four times the image rms noise of 0.061 mJy/beam and increase by factors of two. The restoring beam size is 1.65 × 4.24 mas at −31.7° PA. Right: EVN 22.2 GHz data. The contours start at four times the image rms noise of 0.118 mJy/beam and increase by factors of two. The restoring beam size is 0.23 × 0.81 mas at 55.2° PA.

4. Analysis

4.1. Model fitting

To assess the total flux density S_tot, position (radial coordinates r and ϕ), and FWHM (full width at half maximum) size θ of components, we fit the visibilities directly with Gaussian brightness templates with the built-in modelfit function of difmap. To find the optimal number of components that describe our data, we tested if the reduced χ² of the fit significantly improves adding new components, as described in Schinzel et al. (2012). Furthermore, we checked if adding a new component significantly increased the fraction of the total flux recovered (as estimated from the clean component image), that is, adding more than ∼10% of the total flux. If either of these criteria was fulfilled, we deemed it a significant improvement of our model; otherwise, additional components were discarded.

We identify the two distinct areas in the images with region A and region B, which may consist of multiple Gaussian components. We denote these components with the region name followed by numbers, starting at 0, which designates the brightest component. Additional, unassociated components are named in alphabetical succession. We note that the components are not necessarily cross-identified between all frequencies and epochs; we simply named them in order for each dataset. Should the parameters for regions A and B turn out to be consistent with a lensing system, we also investigate the other regions in comparison, as they could potentially be associated with a higher number of lensed images. To assess the error budget of our fits, we utilised the following formulas:

$$\begin{array}{cc}\hfill \mathrm{S}/\mathrm{N}& =S_{\mathrm{p}}/{\sigma }_{\mathrm{rms}},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{p}}& ={\sigma }_{\mathrm{rms}}·\sqrt{1+\mathrm{S}/\mathrm{N}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {(\mathrm{S}/\mathrm{N})}_{\mathrm{p}}& =S_{\mathrm{p}}/{\sigma }_{\mathrm{p}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{t}}& ={\sigma }_{\mathrm{p}}·\sqrt{1+(S_{\mathrm{tot}}^2/S_{\mathrm{p}}^2)},\hfill \\ \hfill {\sigma }_{\mathrm{tot}}& =\sqrt{{\sigma }_{\mathrm{t}}^2+{(\mathrm{\Delta }g{S}_{\mathrm{tot}})}^2},\hfill \end{array}$$(4)

and

$$\begin{array}{c}\hfill {\theta }_{\lim }=\frac{4}{\pi }\sqrt{\pi log(2)b_{\mathrm{maj}}{b}_{\min }log\left(\frac{\mathrm{S}/\mathrm{N}+1}{\mathrm{S}/\mathrm{N}}\right)}.\end{array}$$(5)

Here, S_p is the component’s peak flux density, as read off in the clean map, σ_rms is the root-mean-square noise measured around the component position in the residual map, σ_p is the error of the peak flux density, Δg are the errors of the absolute flux calibration due to individual telescopes’ gain uncertainties and amplitude calibration, σ_tot is the error on the total flux density, and θ_lim is the minimum size of a Gaussian brightness template that the interferometer with given beam (b_maj and b_min for beam major and minor axis, respectively) and signal-to-noise ratio S/N can resolve (Lobanov 2015). Then,

$$\begin{array}{cc}& \theta =max({\theta }_{\mathrm{fit}},{\theta }_{\lim }),\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\sigma }_{\theta }=\theta /{(\mathrm{S}/\mathrm{N})}_{\mathrm{p}},\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\sigma }_r=\sqrt{b_{\mathrm{maj}}{b}_{\min }+{\theta }^2}/{(\mathrm{S}/\mathrm{N})}_{\mathrm{p}},\hfill \end{array}$$(8)

$$\begin{array}{cc}& \mathrm{and}\hfill \\ \hfill & {\sigma }_{\varphi }=arctan({\sigma }_r/r)·180/\pi .\hfill \end{array}$$(9)

In case θ of a Gaussian component is smaller than θ_lim, the size θ_lim will become an upper limit on the actual size. In the following, we describe how we use all these measured quantities.

4.2. Criteria for milli-lenses

A variety of different sources could in principle show two or more compact components, which we discuss in more detail below. However, gravitational lensing conserves certain properties of the source. So, a potential milli-lens system should fulfil a set of criteria that we define in the following:

1. The morphology of the sources should show two (or more) compact components at all frequencies and epochs studied. Exactly two compact images are expected for the core of a background AGN being lensed by an SMCO when the SMCO acts as a point-mass lens. If the components show extended emission or an elongation, it should be tangential to the critical curve along the Einstein radius θ_E (Henstock 1995). If any elongation is present along the line of separation between images, it has to be small compared to θ_E. Such an elongation can represent the source-intrinsic structure and is consistent with the lensing scenario (as is the case with the curious lensed blazar JVAS B0218+357 Patnaik et al. 1993; Biggs et al. 2003; Spingola et al. 2016).

2. The surface brightness ratio (SBR) between components in region A and B fulfils 0.25 < SBR < 4 in any epoch studied. These are conservative limits that take into account un-modelled systematic uncertainties in imaging and model-fitting (see also Browne et al. 2003). In principle, the SBR should be close to unity for a milli-lens. We tighten our constraints to 0.25 < SBR < 4 compared to Paper I, where we used a threshold of 1/7 < SBR < 7, because of our improved error analysis and better understanding of the data properties in this work.

3. The flux ratio (FR) between regions at a specific observing frequency should not deviate more than the 2-σ combined error between different observing epochs. This reflects the fact that since lensing time-delays are small for small image separations Δθ (Δt  ∼  5(Δθ/mas)² seconds, Press & Gunn 1973) compared to VLBI observing times, any time variability should be seen in both regions in the same way.

4. The spectral index α (defined as S_ν ∝ ν^α) between the components should be comparable (Δα < 0.46). It should be noted that this criterion corresponds to the one employed by Wilkinson et al. (2001), where they allowed a maximum of 30 % deviation between the FR at two different frequencies, considering frequencies at C- and X-band. Comparing FRs between S- and X-band, and C- and K-band, this corresponds to Δα < 0.22 and Δα < 0.17, respectively. While gravitational lensing conserves the shape of the spectrum between the images, it is in principle possible that light rays, travelling along different paths to produce multiple images, cross different plasma regions with different absorption coefficients, leading to differences in their spectra. While this is an important consideration at larger scales, we believe that at scales of < 100 mas considered here, this should be negligible, as it is highly unlikely to have an absorbing plasma screen with size ≪ 100 mas traverse the line of sight, affecting only one image.

5. There should be no significant radial proper motion between components detected (< 2-σ deviation considering the combined error between epochs). Apparent proper motion between putative lensed images can be caused either by motion of the lens or the background source (see, e.g. Garrett et al. 1997; Spingola et al. 2019b, for the latter case). For the first case, Wilkinson et al. (2001) calculated that for milli-lensing the expected change in image separation by a putative lens moving with 1000 km/s transverse to the line of sight is ∼0.01 μas/yr and thus negligible. For the second case, we have to consider two possibilities. First, the background source has a core+knot structure, leading to elongated images (as shown in Henstock 1995). In this case, one might expect source-intrinsic motion of the knot with respect to the core, which then leads to respective shifts within both elongated images. Since we calculate the distance only between the brightest components in each epoch, which we assume to be the stationary core, there should be no radial proper motion detected for a lens system. We note that Spingola et al. (2019b) use a detailed lens model to investigate the proper motion seen in multiple images, which is beyond the scope of this work. Second, the core itself could not be stationary. This is the case if a jet component moving close to or through the core is unresolved, which then leads to an observed shift of the core position (Plavin et al. 2019b). With the magnitude of this shift of the order of 0.2 mas/yr (Plavin et al. 2019b), we estimate that any radial proper motion between the lensed images would be negligible. This is because the separation between two lensed images in the point-mass lens approximation does not strongly depend on the lensing impact parameter, and is always of the order of 2θ_E. This observed motion of the core is also not expected to last longer than typical flares in AGN jets (Plavin et al. 2019b), and should rather manifest as a stochastic jitter of the core position. That means, if significant radial proper motion is detected between components over timescales ≳ 1 yr, it is against the lensing scenario. Here, we focus on criteria with which we can safely reject systems as milli-lenses. We note that proper motion can also potentially be used to confirm a system as such, if the observed proper motion reflects a specific symmetry that is related to the lensing configuration. A detailed analysis of this is dedicated to future work on the remaining lens candidates.

Any such system that does not fulfil any one of the above criteria can be discarded as a milli-lens system. We may introduce more criteria in the future to discard any source listed here that is not already discarded by criteria 1) – 5).

For calculating the surface brightness of a region, we sum up the flux densities and areas of all respective components in that region. The area of a Gaussian component is calculated as A = θ². The parameters are then calculated as follows:

$$\begin{array}{cc}\hfill \mathrm{FR}& =\frac{{\sum }_i{S}_{\mathrm{tot},\mathrm{A},\mathrm{i}}}{{\sum }_i{S}_{\mathrm{tot},\mathrm{B},\mathrm{i}}},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill \mathrm{SBR}& =\frac{{\sum }_i{S}_{\mathrm{tot},\mathrm{A},\mathrm{i}}/{\sum }_i{\theta }_{\mathrm{A},\mathrm{i}}^2}{{\sum }_i{S}_{\mathrm{tot},\mathrm{B},\mathrm{i}}/{\sum }_i{\theta }_{\mathrm{B},\mathrm{i}}^2},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill d& =\sqrt{r_A^2+r_B^2-2r_A{r}_{B}cos({\varphi }_1-{\varphi }_2)},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill \alpha & =\frac{log(S_{{\nu }_2}/S_{{\nu }_1})}{log({\nu }_2/{\nu }_1)},\hfill \end{array}$$(13)

where FR denotes the flux ratio, SBR the surface brightness ratio, d the distance between components, and S_ν1 and S_ν2 are the sum of the total flux densities of all components in a region at frequencies ν₁ and ν₂, respectively. Here, we define the spectral index α as S_ν  ∝  ν^α. To calculate the distance d between regions in case of multiple components making up one region, we calculated the distance between the brightest components. All errors on these quantities are calculated with standard error propagation. We did not propagate the gain error in Eq. (4) for the absolute flux calibration in the uncertainties for FR and SBR, since they should scale the same for both components and the ratios should be unaffected. The observed brightness temperature of the components and regions has then been calculated as

$$\begin{array}{c}\hfill T_{\mathrm{b},\mathrm{obs}}=\frac{c^2}{2k_{\mathrm{B}}{\nu }_{\mathrm{obs}}^2}\sum _i{S}_{\mathrm{tot},\mathrm{i}}/\sum _i{\theta }_{\mathrm{i}}^2,\end{array}$$(14)

where c is the vacuum speed of light, k_B the Boltzmann constant, and ν_obs is the observing frequency.

The most common contaminant object mimicking a potential milli-lens system is expected to be an AGN core-jet structure, which would reveal itself by not following criteria 1), 2) and 3). In addition, criterion 4) will not be fulfilled, as the spectral indices of the two components would be different because the expanding jet is expected to be more optically thin (α ≲ −1) and the core optically thick (α ≳ −0.5).

Furthermore, the sources studied here could fall under the class of compact symmetric objects (CSOs; see Readhead et al. 1993; Wilkinson et al. 1994), which are, as with our milli-lenses, also morphologically defined as having two (compact) components detected on VLBI-scales. These components might also be straddling a central core. The observed symmetry in these objects, combined with their low abundance among radio-loud AGN as well as their lack of evidence for relativistic beaming made previous authors conclude that these are either very slowly expanding old objects (see, e.g. the reviews by Orienti 2016; O’Dea & Saikia 2021) or very short-lived objects, with new evidence pointing to the latter scenario. This has been found in the most complete study of CSOs to date, conducted in a series of papers (Kiehlmann et al. 2024a,b; Readhead et al. 2024; de la Parra et al. 2024). Kiehlmann et al. (2024a) have set out the following criteria to define CSOs, complementing and streamlining earlier definitions: i) no projected radio structure larger than 1 kpc; ii) evidence of emission on both sides of the centre of activity; iii) no fractional variability greater than 20 %/yr; iv) no evidence for apparent superluminal motion in any jet component in excess of v_app = 2.5c. These criteria apply broadly to how we defined milli-lens systems above. There is no apparent reason that CSOs should fulfil criterion 2) for a given epoch and observing frequency; however, if they do, they will most likely fulfil this criterion for all epochs due to their lack of strong variability. In this case, the spectra of the components can become a discriminating factor, as CSOs are likely to have a gigahertz peaked spectrum (GPS; 82 % of sources in the sample of Tremblay et al. 2016), and expected to steepen at frequencies above a few gigahertz that we are investigating here. In contrast, milli-lenses are expected to be two images of the same flat-spectrum core. However, in a new study, de la Parra et al. (2024) find the full range of spectral types in a complete sample of CSOs, with about half of the sample of 17 sources showing a flat-spectrum between 5 and 8 GHz (α > −0.5), and the other half a steep spectrum (α < −0.5). Following this, they argue that CSOs and milli-lenses can only be distinguished for a putative milli-lens having α > −0.3. In addition, CSOs can show significant proper motion, as long as their apparent speed stays below v_app = 2.5c, while for milli-lenses, no proper motion is expected (Wilkinson et al. 2001, see again 5). A prime example of a good CSO candidate found in our study can be seen in Fig. 3, with the derived parameters shown in Table 3.

Fig. 3. J0616-1957 EVN 4.9 GHz data. The contours start at four times the image rms noise of 0.142 mJy/beam and increase by factors of two. The restoring beam size is 1.57 × 3.74 mas at 66.9° PA. The source is not detected at 22.2 GHz. Despite being rejected as a milli-lens candidate, it shows clear morphological evidence of being a CSO with its extended radio lobes, which were only revealed through our new EVN observations.

Initially, we restricted ourselves to flat-spectrum components as expected for milli-lenses, but as shown in de la Parra et al. (2024), CSOs can also show flat (α > −0.5) spectra at frequencies well above the spectral peak. Although it seems highly unlikely to observe a steep-spectrum compact component that is isolated, and then gravitationally lensed, it is not impossible if the core of the lensed background AGN is strongly self-absorbed. Currently, studies are underway to carry out a large survey of steep-spectrum CSOs, covering the unlikely scenario of the hotspot of a CSO being lensed as well.

In addition, steep-spectrum sources like CSOs are predominantly found at redshifts below one (Kiehlmann et al. 2024a), while lensed sources on arcsecond-scales typically have z > 1 (Browne et al. 2003). This remains indicative, and we do not rule out lens candidates based on their redshifts at this point.

A third alternative is supermassive black hole binaries (SMBHBs). These objects are expected as the result of galaxy mergers, where the central SMBHs of both galaxies will spiral towards one another (Merritt & Milosavljević 2005). If both are accreting and possibly also producing extended jets, we can expect to see the radio emission of both of them. The only known example of an SMBHB discovered this way with a separation of 7 pc is the case of 0402+379, described in Rodriguez et al. (2006). Another example of a close dual AGN system detected with VLBI is described in Deane et al. (2014) with a separation of 26 mas (138 pc). For SMBHBs, criterion 4) will likely be fulfilled, as both components are expected to show an optically thick radio core. However, there is no reason that our other criteria should all be fulfilled. Maness et al. (2004) demonstrated that monitoring over several years can constrain the variability of two flat-spectrum components to discriminate between the SMBHB and the lensing scenario. In addition, to confirm (or reject) a system as SMBHB, broad-band characterisation is necessary. For example, in an SMBHB, double-peaked emission lines are expected (Deane et al. 2014).

In conclusion, we believe that our criteria create a fine enough sieve to find and distinguish milli-lens systems. We categorise all our candidates based on the applicability of these criteria into: being still a viable milli-lens candidate; core-jet; CSO; SMBHB candidate; or other.

5. Results

Here we present the results of our study. Figure 4 shows the 4.9 GHz image for the source J2347-1856, as an example of a source that we cannot discard as a milli-lens candidate at this stage. Table 4 lists all the derived properties with their respective errors. The images and derived quantities for the other sources are available on Zenodo.

Fig. 4. J2347-1856 EVN 4.9 GHz data. The source is not yet rejected as a milli-lens but has strong indications to be a CSO. The contours start at four times the image rms noise of 0.296 mJy/beam and increase by factors of two. The restoring beam size is 1.55 × 4.44 mas at 68.91° PA. The source is not detected at 22.2 GHz.

We summarise our categorisation of lens candidates in Table 5. In total, we rule out 31 sources as milli-lens systems from our initial 40 candidates. We discuss the sources rejected as milli-lens candidates in more detail in Sect. 6. As can be seen from Fig. 5, the most discriminating factors in rejecting milli-lens systems are the spectrum test (23 sources rejected), the SBR test (23 sources rejected) and the morphology test (nine sources rejected). Less discriminating are the FR test (three sources rejected) and the separation test (one source rejected). Investigating both variability (FR test) and proper motion (separation test) is expected to be more discriminating with the addition of more epochs of observations, as we have only two epochs for 19 of our sources. In total, 19 out of 31 rejected lens candidates fail more than one of our criteria.

Fig. 5. Summary charts for milli-lens rejection criteria. Left: Venn diagram showing which combination of criteria lead to source rejection. Right: Bar diagram showing the total number of sources rejected for each criterion alone.

6. Discussion

We searched the literature and various databases for additional data on our 40 candidates. For 25 of the 40 sources, we found optical counterparts with available redshifts within 1″. Optical images are presented in Figs. C.1 and C.2. Six of those sources have spectroscopic redshifts available. For another six sources, we found optical counterparts in Gaia DR3 (Gaia Collaboration 2016, 2023).

6.1. Gaia-VLBI offsets

To get the offsets of the Gaia photo-centre to the VLBI images, we compared the J2000 source positions in the Radio Fundamental Catalog² (RFC), the International Celestial Reference Frame 3 (ICRF3; Charlot et al. 2020) and the VLBI Celestial Reference Frame Solution OPA2024a³ close to our observing epochs with the Gaia coordinates. The Gaia coordinates are given in J2016, which differs from J2000 only if significant proper motion is observed. While this correction is mainly necessary for stars, some extragalactic sources also show significant proper motion, which may be due to the motion of bright jet components compared to the jet base. Since we want to compare the positions of the VLBI core with the jet base seen in the optical, we do not apply any corrections to the J2016 coordinates taken from Gaia to convert them to J2000.

Gaia-VLBI offsets have been systematically studied in various works. Plavin et al. (2019a) studied the offset in over 1000 AGN, comparing VLBI positions from the RFC and Gaia DR2. They find a significant trend of the offset being either upstream or downstream of the jet direction, supporting the fact that the optical emission is either dominated by the accretion disc (upstream offset) or the jet (downstream offset). Lambert et al. (2024) find a similar trend, comparing model-fit positions from the MOJAVE survey to the photo centres from Gaia DR3. Both studies find a trend for downstream offsets to occur more often in flat spectrum radio quasars (FSRQs) compared to BL Lac objects, supporting the fact that the former have more powerful accretion discs. Variability in the optical jets could cause a jitter in the VLBI-Gaia offset (Petrov & Kovalev 2017) and may be used to constrain the location of optical flares. Proper motion of jet components in the optical may be traced with future data releases from Gaia. Offsets between VLBI and Gaia positions may also be used to investigate possible SMBHB systems. Chen et al. (2023) identified a few SMBHB candidates with VLBI follow-up from an initial selection of candidates based on variability-induced astrometric jitter observed with Gaia (Hwang et al. 2020).

Out of the six sources with measurable VLBI-Gaia offsets in this work, three (J0132+5211, J0237+1116 and JJ2114+4036) show offsets clearly along the direction of the radio jet. These systems are also more likely to classify as core-jet systems instead of being CSOs, and are rejected as milli-lens systems (see Sect. 6.2). The other three sources show different offsets, and their interpretation is beyond the scope of this paper. Aside from being caused by potential dual or offset AGN systems, the VLBI-Gaia offsets could be caused by a variety of effects, including contamination by the emission of the host galaxy, dust obscuration affecting the observed photo-centres, or even Gaia observing foreground stars.

6.2. Notes on individual sources

6.2.1. J0010-0740 (0008-079)

The source (Fig. A.1) is unclassified in the literature. We discarded this source as a milli-lens because of the different spectra of the components and the high SBR of 6.82 at C-band and 22.17 at X-band, respectively. The rejection as a lens candidate and classification as a CSO is also affirmed by the detection of a central component C between regions A and B. The spectra of the two regions A and B are notably different, with A showing a rather flat spectrum, while B has a steep spectrum. Since component C is not identified at the 2015 C-band epoch, we cannot calculate its spectrum to see if it could be a flat-spectrum core straddled by components A and B. The source is a good CSO candidate, the VLBI fractional flux variability being < 1 %/yr. J0010-0740 has an optical counterpart classified as a galaxy with 22.17 optical r-band magnitude and photometric redshift of z = 0.31  ±  0.09 in the Sloan Digital Sky Survey (SDSS DR17; Abdurro’uf et al. 2022). The linear distance between components A and B at C-band is then ∼61 pc, consistent with the CSO definition.

6.2.2. J0011+3443 (0008+344)

The source (Fig. A.2) has an optical counterpart with 23.2 mag at r-band with a spectroscopic redshift of z = 0.89 (Vardoulaki et al. 2010) and is associated with a compact, steep spectrum (CSS) radio source. This source was discarded as a lens candidate because of the low SBR at C-band of 0.09 in the 2020 epoch. The fractional variability in our data is about 9 %/yr, consistent with a CSO. The distance between components B and A0 also does not change over a period of more than 4 years. The spectra are remarkably similar, and steep. The linear distance between the components is ∼315 pc. All of this favours this source to be a CSO.

6.2.3. J0024-4202 (0022-423)

The source (Table A.3) has previously been classified as a GPS radio galaxy (e.g. Labiano et al. 2007), with the radio spectrum peaking at 1.6 GHz. GPS sources are likely to fall into the CSO category as well (Tremblay et al. 2016). The optical counterpart (20.3 r-band magnitude; de Vries et al. 1995) has a spectroscopic redshift z = 0.937 and is classified as quasi-stellar object (QSO; Labiano et al. 2007). The integrated radio spectrum of the source appears rather flat around frequencies up to a few gigahertz, but it seems to become considerably steeper at higher frequencies (Healey et al. 2007). D’Abrusco et al. (2014) classify the source as FSRQ based on matching known γ-ray blazars with their IR colours measured with the Wide-field Infrared Survey Explorer (WISE).

We did not observe this source with the EVN due to its low declination. The epochs we have from Astrogeo show two similar and steep spectra for both components. The variability between the first two X-band epochs is much more than 20%/yr, while for the other epochs, it is within about 10%. We found that, for the early X-band epochs in Astrogeo, the weaker component is found very close to the edge of the field of view in the map provided in Astrogeo, and thus we are cautious in interpreting the observed high variability in these datasets, and still judge this source as a CSO candidate. The linear distance between the components is ∼230 pc, consistent with being a CSO. The previous classification as blazar has to be further scrutinised with future observations in terms of radio variability and possible superluminal motion. Ultimately, we do not reject this source as a lens candidate at this point.

6.2.4. J0044+2858 (0041+287)

The source (Fig. A.3) is unclassified in the literature. We found an optical counterpart in SDSS classified as a galaxy with 21.57 r-band magnitude and photometric redshift of z = 0.75  ±  0.07. The source has two steep-spectrum components in the VLBI image, with almost no variability (< 4%/yr) observed over 5 years between the two available epochs. We discarded this source as a milli-lens because of the different component spectra. The linear distance between the identified components A and B is ∼38 pc. We classify it as a good CSO candidate.

6.2.5. J0052+1633 (0049+162)

For this source (Fig. A.4) we found an optical counterpart with 22.3 r-band magnitude in SDSS classified as a galaxy with photometric redshift of z = 0.43  ±  0.04. That means the linear distance between components A and B is ∼49 pc. The spectra are notably different, with component B having a much flatter (or even inverted) spectrum compared to component A. This trend continues for the spectrum between 5 and 22 GHz, although the difference is not as prominent. The source was discarded as a gravitational lens system because of the spectra, but also because of the too-large FR variations of 0.87 compared to the combined error of 0.81 at C-band, and the high SBR of 6.44 observed at K-band. The morphology also does not support the lensing scenario, since component A shows some elongated emission towards component B. For all frequencies and epochs, the variability is maximally about 20%/yr. The source might still classify as a CSO since for a core-jet system, the brighter component would be expected to be the flat-spectrum core and more variability is expected. The component separation of ∼49 pc also fits into the CSO scenario.

6.2.6. J0118+3810 (0115+379)

The source (Fig. A.5) is unclassified in the literature. There is an optical counterpart with r = 27.3 mag from the Blazar Radio and Optical Survey (BROS) catalogue (BROS J0118.1+3810; Itoh et al. 2020), but no available redshift. The two steep-spectrum components, as well as the lack of proper motion and variability (∼1(1 + z) %/yr) over 4 years, strongly suggest the source to be a CSO. It was ultimately rejected as milli-lens candidate because of the SBR being above the threshold of 4 within the uncertainties for all epochs.

6.2.7. J0132+5211 (0129+519)

The source (Fig. A.6) has been classified as an FSRQ in D’Abrusco et al. (2019) based on its IR colours in WISE, which were selected to match those of known γ-ray emitting blazars. An optical counterpart with photometric z = 0.98  ±  0.16 is found in Richards et al. (2015), and the r-band magnitude listed in SDSS is 19.97. The spectra of both components are steep but significantly different. Also, the SBR fluctuates between being below the threshold of 0.25 and above the threshold of 4, failing the source as a milli-lens candidate. Component A varies maximally by 4%/yr, component B by 21%/yr, and the total flux by 6%/yr. Surprisingly, the distance between the components has decreased more than the 2-σ combined error of 1.54 mas over the course of 4 years. That is also inconsistent with a milli-lens. In addition, the rejection as a milli-lens is supported by the position of the Gaia photo-centre (see Fig. A.6), found in Gaia DR3 with 19.85 G-band magnitude, whose position is very close to component A (2.5 mas away). This hints at the core-jet nature of the source, with component A being the core. The linear distance between the components is ∼1.2 kpc at the given redshift, which puts the source outside the CSO definition. The source is more likely a core-jet system, although in that case, the large size would still be puzzling, and more variability should be expected. The source remains a curious case and more data will shed light on the nature of this interesting source.

6.2.8. J0139+0824 (0137+081)

The source (Fig. A.7) is unclassified in the literature. There is an optical counterpart in SDSS, classified as a galaxy with 21.10 mag at r-band, and there is a photometric redshift z = 1.05  ±  0.11 from DES DR8 (Dark Energy Survey; Dey et al. 2019). In the VLBI image there are two steep-spectrum components, separated by ∼130 pc, with some indication of emission towards the centre between them. This and the lack of significant variability (∼4%/yr) of this source over four years suggests it to be a CSO. The difference between the spectral indices between components of Δα = 0.47 is larger than our defined threshold and thus does not fit a milli-lens system.

6.2.9. J0203+3041 (0200+304)

The source is associated with an optical counterpart classified as a galaxy in SDSS, and there is also a Fermiγ-ray source associated (Massaro et al. 2015). The source was discarded as a milli-lens because of the high SBR (≫ 4) observed at K-band. Also, a new component C between region A and B appears in the new EVN C-band image (Fig. A.8), and regions A and B slightly extend towards each other. This is highly indicative of the source being a CSO, although not enough to reject the source based on criterion 1. The redshift listed in Massaro et al. (2015) is z = 0.76, while Ackermann et al. (2015) list it at z = 0.96 (spectroscopic, originally in Shaw et al. 2012) and categorise it as FSRQ based on optical classification. The magnitude listed in SDSS is 21.84 at r-band. The NVSS image shows two sources separated by ∼1′. In the MOJAVE survey (Lister et al. 2018), there is a source about 0.7′ away from J0203+3041 with the same B1950 source identifier, which is probably why these MOJAVE data are falsely listed on the Astrogeo webpage for J0203+3041. This other source is also listed at z = 0.76 and is categorised as BL Lac object (Ackermann et al. 2015). We decided to adopt a value of z = 0.96 for J0203+3041. Then, the linear distance between components is ∼322 pc. The spectra of both components are steep. The variability is ≫20%/yr between the two close S-band epochs but is otherwise observed to be between 10 and 30%/yr. While the morphology and the spectra are suggestive of a CSO, the variability is too high to classify as such. The source is more likely a core-jet system, which is also more likely to be γ-bright. However, we note that in Kiehlmann et al. (2024a), three out of the 79 bona-fide CSOs are also detected in γ-rays, most likely originating from the radio lobes.

6.2.10. J0210-2213 (0207-224)

The source (Fig. A.9) was previously classified as a GPS source (Snellen et al. 2002), peaking at 1.5 GHz. The optical counterpart is located at redshift z = 1.491  ±  0.003 (spectroscopic; de Vries et al. 2007) with a magnitude of 23.52 at r-band (Snellen et al. 2002). The linear distance between components is ∼247 pc, falling within the CSO category. The two steep component spectra fit well as a CSO, but we also cannot reject the source as a milli-lens at this stage. The fractional variability is just within ∼20%/yr, barely passing as a CSO.

6.2.11. J0213+8717 (0159+870)

This source (Fig. A.10) has an unclassified optical counterpart with z = 1.5 and 19.6 r-band magnitude found in Milliquas v7.2. Since it is not listed in Milliquas v8 (Flesch 2023) any more, we decided to consider it as an unknown redshift. We rejected the source as a milli-lens candidate due to the high SBR of 19.8 seen in the EVN C-band data, and the very different spectra between components. Richards et al. (2014) studied the radio variability of a sample of mostly γ-ray detected blazars and found a modulation index of 0.077 (defined as the standard deviation of the flux density measurements in units of the mean measured flux density) in the 15 GHz OVRO light-curve for this source. The maximum fractional variability we found is ∼10(1 + z) %/yr. The spectra of the two components are different, with component B having a much flatter, or even inverted spectrum. While not significant within the errors, there is an indication of an increase in separation between the components over time. The observed component brightness temperatures are also partly in excess of 10¹¹ K, a typical value for the equipartition brightness temperature T_eq in AGN (Readhead 1994; Liodakis et al. 2018). Values above T_eq indicate Doppler boosting due to the source moving towards the observer close to the line of sight. All this indicates this source to be a core-jet system.

6.2.12. J0222+0952 (0219+096)

This source (Fig. A.11) is unclassified in the literature, with no optical counterpart. We observed the source at C-band in two sessions. It shows a peculiar structure with three main regions A, B and C, and some emission between them, that encompasses a significant fraction of the separation between the regions. While the source does not fulfil the milli-lens criteria for morphology, SBR, and spectra, it might still classify as both a CSO or core-jet. However, the variability we see between the two EVN epochs at C-band is much in excess of 20%/yr, excluding the CSO possibility. The source appears point-like on VLA scales with a rather flat spectrum of α  ∼   − 0.4 (Vollmer et al. 2010), also more suggestive of a core-jet system.

6.2.13. J0232-3422 (0230-345)

This source (Table A.13) is unclassified in the literature, with an optical counterpart found in DES DR2 (J023230.02-342203.8) with 24.51 mag at r-band, ∼0.3″ away from the VLBI position. We did not yet obtain follow-up data for the source due to its low declination which is inaccessible for the EVN. The spectral indices of the two components differ by 0.76, failing the source as a milli-lens system. Since there is only one epoch available, we cannot make further statements about the nature of the source at this point, and further analysis of archival data and new observations will shed light on this.

6.2.14. J0237+1116 (0234+110)

This source (Fig. A.12) was already briefly described and ruled out as a lens in Paper I. It is associated with an optical counterpart in Gaia with 17.91 G-band magnitude, and we found a spectroscopic redshift of z = 0.985  ±  0.067 listed in Storey-Fisher et al. (2024). The Gaia photo-centre is close to component A with a separation of 6.1 mas to the south-west. We reported a different offset with a smaller magnitude in Paper I, where we used the converted J2000 coordinates (see the beginning of Sect. 6). The difference is due to some proper motion seen between the J2000 and J2016 epochs, which may suggest some moving jet component. The direction of the offset vector is close to the jet direction as estimated from the separation vector between the two components A and B, which further hints at the source being a core-jet system. In addition, we could discard this system as a milli-lens due to its high SBR at X-band of 8.65. The fractional variability of the source is 17%/yr, while both components show flat spectra. The linear distance between components is ∼41 pc. This source might be classified as either a core-jet or a CSO, and more observations are needed to study the variability and spectra.

6.2.15. J0502+1626 (0459+163)

This source (Fig. A.13) is unclassified in the literature, and we could not find an optical counterpart. The source has been discarded as a milli-lens due to its high SBR at X-band of 8.26 and the significant difference in the spectra of components. No strong variability is observed, but there is an indication of an increase in the distance between the components. The source could be a core-jet or a CSO.

6.2.16. J0527+1743 (0524+176)

This source (Fig. A.14) is unclassified in the literature. It has an optical counterpart in SDSS classified as a galaxy with 22.47 mag at r-band and photometric redshift z = 0.46  ±  0.08. We discarded this source as a milli-lens because of its SBR being < 0.24 at X-band. The linear distance between components is ∼107 pc. The fractional variability is in slight excess of 20%/yr, but we still judge this source as a good CSO candidate. Since both spectra are flat, the source could also qualify as SMBHB candidate.

6.2.17. J0616-1957 (0613-199)

This source (Fig. 3) is unclassified in the literature. We discarded this source as a milli-lens because of the high SBR at X-band of 8.49, the different spectra between components, and the clear extension of both emission regions towards the centre. The steep spectra and the edge-brightening of the two hotspots are typical for CSOs. There is no available optical counterpart and redshift, but at a fiducial redshift of 1, the source would have a separation between components of ∼200 pc. The fractional variability of the source is < 2(1 + z) %/yr. We deem this source a CSO.

6.2.18. J0732+6023 (0728+605)

This source (Fig. A.15) is classified as a flat-spectrum radio source measured around gigahertz frequencies (Massaro et al. 2014). Our image shows an extended and bright region A, that is model-fitted with three components, and only component B comprising region B. We discarded the source as a milli-lens based on the very different spectra, where region A is very steep and component B has a very flat spectrum. It also did not pass our SBR test at X-band, being with > 23 much in excess of our defined threshold. The maximum fractional variability of the source that we found is ∼7%/yr. This source is a CSO candidate, but could also still be a core-jet system. The linear distance between the components is ∼163 pc, based on the photometric redshift of the optical counterpart of 1.064  ±  0.162 found in DES DR8.

6.2.19. J0923-3435 (0921-343)

This source is unclassified in the literature, and we did not obtain new data due to the source’s low declination which is inaccessible to the EVN. The integrated radio spectrum on VLA-scales peaks at around 1 GHz (Vollmer et al. 2010), which is consistent with the data from Astrogeo (Table A.18). The source could be classified as a CSO, but more observations are needed to study the possible proper motion and variability of the components. Without additional data, we cannot reject the source as a milli-lens at this stage. Analysing more archival data as well as possible new observations will shed light on the nature of this source.

6.2.20. J1132+5100 (1130+512)

We discarded this source (Fig. A.16) as a milli-lens due to its SBR being 0.17 and 0.12 at two different epochs, the different spectra, and the extension of emission region A towards B, which is revealed in the new EVN C-band image and encompasses a significant fraction of the distance between the two regions. With a photometric redshift of z = 1.152  ±  0.303 (DES DR8), the linear distance between components is ∼116 pc, consistent with a CSO. This is also supported by the rather low variability of ∼7%/yr.

6.2.21. J1143+1843 (1140+188)

This source (Fig. A.17) has previously been reported to be a CSO (Sokolovsky et al. 2011) and an SMBHB candidate (Tremblay et al. 2016), which we confirm with our analysis. The rich amount of archival available datasets together with newly obtained data for this source will be discussed in an upcoming, dedicated publication (Álvarez-Ortega et al., in prep.). At this stage, we are not rejecting this source as a lens candidate.

6.2.22. J1218-2159 (1216-217)

This source (Fig. A.18) is unclassified in the literature. There is no optical counterpart. This source is still a viable milli-lens candidate, passing all our current tests. With the overall variability being maximally ∼8(1 + z) %/yr, the source is also a good CSO candidate. This source will be studied in more detail concerning spectra and variability with future, more sensitive observations.

6.2.23. J1306+0341 (1303+039)

This source (Fig. A.19) has no optical counterpart. We discarded the source as milli-lens because of the high SBR of 8.38 and 7.42 in two epochs. The source morphology, spectra and the ∼8(1 + z) %/yr fractional variability are very indicative of the CSO nature.

6.2.24. J1340-0335 (1337-033)

This source (Fig. A.20) has an optical counterpart and associated photometric redshift of z = 0.977  ±  0.489 in DES DR8. The source has been identified to peak at low radio frequencies below 1 GHz (Callingham et al. 2017). We discarded the source as milli-lens because of the high SBR of 9.14 and 12.1 in two X-band epochs and the new EVN epoch having a high SBR of 12.09. The source morphology, showing three emitting regions that are clearly connected, also does not fit a gravitational lens system. The spectra support this source being a CSO, while the variability is only in slight excess of 20%/yr, so we retain is as a CSO candidate. The linear distance between components, which is ∼58 pc, also fits the CSO scenario.

6.2.25. J1344-1739 (1341-174)

This source (Fig. A.21) has also been identified to peak at low radio frequencies below 1 GHz (Callingham et al. 2017), and has no optical counterpart. This source failed our milli-lens criteria for the spectra, whose spectral indices differ by 0.53, and the fact that region A shows a clear eastward extension towards region B in both the EVN C-band image and different Astrogeo images. With the variability being less than 10(1 + z) %/yr, this source is a good CSO candidate.

6.2.26. J1632+3547 (1630+358)

Tremblay et al. (2016) discarded this source (Table A.25) as a CSO candidate, but they still propose it as SMBHB candidate. The source has a photometric redshift of z = 1.081  ±  0.082 (DES DR8) and an r-band magnitude of 23.70 in the BROS catalogue. The linear distance between components is ∼91 pc. Since we did not obtain new data for this source, we cannot make any more definitive statements at this point. Because of the lack of multi-epoch observations, we can also not rule out the source as a milli-lens yet. Analysing more archival data and proposing for new observations should provide more clarity.

6.2.27. J1653+3503 (1652+351)

This source (Table A.25) has been discarded as both CSO and SMBHB candidate (Tremblay et al. 2016). It has an optical counterpart in SDSS with 22.85 r-band magnitude and photometric redshift z = 0.81  ±  0.10. The linear distance between components is ∼78 pc. As for J1632+3547, we did not obtain any data from follow-up observations and cannot discard this source as a milli-lens. A deeper analysis of archival data and new observations will constrain the nature of this source.

6.2.28. J1721+5207 (1720+521)

This source (Fig. A.22) has an optical counterpart in DES DR8 with photometric redshift z = 1.302  ±  0.229. It has a steep integrated spectrum at frequencies around 1 GHz (Massaro et al. 2014), which continues to steepen at higher frequencies, as is confirmed by our data. We rejected this source as milli-lens because of its high SBR of 12.47 for our new EVN data, as well as the spectral indices differing by 0.64 in the 2014 epoch, above our threshold. At the given redshift, the linear distance between components is ∼408 pc. The variability is within 7%/yr, making this source a good CSO candidate.

6.2.29. J1805-0438 (1802-046)

This source (Fig. A.23) has an IR and optical counterpart with 24.4 mag at r-band listed in the OCARS catalogue (Malkin 2018). It is listed as BL Lac object in D’Abrusco et al. (2019) according to its IR colours in WISE, which were selected to match those of known γ-ray emitting blazars. Our VLBI images show three components at lower frequencies, but for some epochs from Astrogeo at X-band, only a subset of any of the components A, B, and C is detected in the image, making the interpretation more difficult. Our new EVN image at 4.9 GHz shows all three components detected, while the 22.2 GHz image shows only component C. Component B seems undetected, which is not surprising given the steep spectrum and the partial non-detection of component B in earlier epochs. The separation vectors between components A and B and components A and C are almost perpendicular to each other, giving the source a peculiar morphology. Because of the uncertain component identification across epochs, the source passes the FR test, although some epochs at S-band deviate more than the uncertainties allow. We attribute this to uncertainties introduced by blending effects due to the low resolution at S-band. The discriminating factors to discard this source as a lens candidate would be the clearly different spectra of components A and B. As an additional test, we have considered the possibility that components A and B comprise one elongated image, and component C the other, less elongated image of a milli-lens system. In this case, the morphology would match the expectations, since the elongations would be along the critical curve. Then, the corresponding values for SBR, FR, and the spectrum are all consistent with our criteria within the uncertainties. For the spectrum, only the 2012 epoch has been considered. Hence, we keep the source in the list of lens candidates at this stage, and we aim at more sensitive X-band observations to shed light on this peculiar source. The small separation between components makes the source likely to be a CSO but, for confirmation, we would need a redshift. The fractional variability in our data is ∼10(1 + z) %/yr.

6.2.30. J2010+1513 (2007+150)

The source (Fig. A.24) does not have an optical counterpart. We observed this source in two sessions with the EVN. While the regions seem blended for the March 2021 C-band epoch due to the large elliptical beam, this is resolved in the June epoch, showing clearly two components for each region and their extension towards each other. The source passes all our criteria but the spectrum, which is significantly different between the regions, and the morphology. Only component A is detected at 22.2 GHz, indicating that the spectrum of region B must steepen even further than what is indicated between 4.3 and 7.6 GHz. The variability of the total flux between our EVN epochs is ≲ 20(1 + z) %/yr, barely fitting the CSO definition. Since there is no evidence of proper motion, this source is still likely to be a CSO.

6.2.31. J2044+6649 (2044+666)

This source (Fig. A.25) is unclassified and has no optical counterpart in the literature. We discarded the source based on the spectral indices differing by more than 0.5. The FR is remarkably similar between the components across epochs, and the variability is < 4(1 + z) %/yr. This source is likely to be a CSO.

6.2.32. J2114+4036 (2112+404)

This source (Fig. A.26) has an optical counterpart identified as a galaxy with photometric redshift z = 0.117  ±  0.015 (Bilicki et al. 2014). We discarded this source as a milli-lens system because of the different spectra between the regions in our new 4.9 and 22.2 GHz data and the appearance of central emission between regions A and B, connecting the two. Region A steepens considerably, while region B stays rather flat, and could be the core of a core-jet system. The flux variability of component A exceeds 30%/yr, making this scenario more likely, although the total flux density varies maximally by ∼14%/yr. The Gaia counterpart with 21.39 G-band magnitude is found very close to component B, with an offset of 25 mas from the VLBI position, which is centred on the brighter component A. For a milli-lens, we should expect the photo-centre in between the lensed images, but here the Gaia position indicates the optical emission to rather come from a jet component, or the accretion disc if component B is identified as the core. The linear distance between the regions is ∼38 pc.

6.2.33. J2209+6442 (2208+644)

This source (Fig. 2) has no optical counterpart. It was rejected as milli-lens because of the low SBR of 0.05 at X-band, as well as the different spectra. The flat spectrum of component B up to 22.2 GHz could indicate it being the core of a core-jet system. The low variability (maximal 4(1 + z) %/yr) and the lack of detected proper motion also suggest a possible CSO scenario.

6.2.34. J2214-2521 (2211-256)

This source (Fig. A.27) has an optical counterpart classified as an extended galaxy in Bilicki et al. (2014) with a spectroscopic redshift z = 0.0868  ±  0.0002 (Colless et al. 2003). There is an associated Gaia counterpart with 20.99 mag in G-band, the photo-centre being 152 mas offset to the south-east from component A. We already discussed the optical properties of this source briefly in Paper I. Finally, we discarded it as milli-lens because of the different spectra in the 2017 epoch and the low SBR of < 0.06 for the 2017 X-band data. The new EVN image shows a new component C compared to the Astrogeo data, and a much lower separation between components A and B than in previous epochs. We deem this an uncertainty in the imaging and model-fitting process due to bad data quality in this case and do not discard the source because of the separation test. In the last two X-band epochs, the source shows fractional variability in slight excess of 20%/yr, making it less likely to be a CSO and rather a core-jet, but we do not exclude the CSO possibility at this stage. More sensitive observations are required to determine the source morphology and variability more clearly.

6.2.35. J2225+0841 (2223+084)

This source (Fig. A.28) has an optical counterpart classified as a galaxy in SDSS with 23.16 mag at r-band and photometric redshift z = 0.73  ±  0.10. The linear distance between components is thus ∼33 pc. We rejected this source as a milli-lens based on the components’ different spectra, and the high SBR ≫ 4 observed in both the EVN C-band and K-band data. The spectral indices for each of the individual components are remarkably similar to what is found from the Astrogeo data between 4.3 and 7.6 GHz, and the EVN data between 4.9 and 22.2 GHz. The source is likely a CSO. The maximum fractional variability of ∼17%/yr also fits this scenario.

6.2.36. J2259+4037 (2256+403)

This source (Fig. A.29) has no optical counterpart. It shows a flat spectrum around 1 GHz for the integrated flux (Massaro et al. 2014). We discarded this source as a milli-lens because of its very different spectral indices (Δα > 1) in our new 5 and 22 GHz data, very low SBR of 0.01 at K-band, and the source morphology which suggests a connecting bridge between the regions in both new C-band images. The source varies maximally by ∼16(1 + z) %/yr and might thus qualify as a CSO.

6.2.37. J2312+0919 (2309+090)

This source (Fig. A.30) has an optical counterpart classified as a galaxy with 18.22 r-band magnitude and photometric redshift z = 0.23 ± 0.08 (SDSS). It may be associated with the Fanaroff-Riley II (FRII) high excitation radio galaxy, 3C 456 (e.g. Macconi et al. 2020), which is located at redshift z = 0.23 and about 2.5″ (8 kpc at the given redshift) away from J2312+0919. The photo-centre of the Gaia counterpart (19.4 G-band magnitude) is significantly offset to the north of the source by 24.7 mas, not supporting a milli-lens scenario if it is associated with J2312+0919. Ultimately, we discarded this source as a milli-lens because the FR is too variable at C-band, the difference being 0.81 while the 2-σ combined error is 0.80. The maximum fractional variability is just within limits for the CSO classification (∼19 %/yr), and the linear distance between components is ∼21 pc, fitting this scenario. Both components are dimming over time, and there is an indication of an increase in separation. More monitoring of this source would prove valuable to constrain the trend in flux density, the separation, and if a steepening of the spectra with time could be observed.

6.2.38. J2324-0058 (2321-012)

This source (Fig. A.31) has an optical counterpart in SDSS with 21.82 r-band magnitude and photometric redshift z = 0.37  ±  0.09 and is classified as a galaxy. The linear distance between components is ∼26 pc. The source passes only the separation test, since it is too variable in FR (the maximum difference of 1.24 compared to the combined error of 0.64), and shows an SBR of 7.11 in the new EVN data. The spectra are also very different, with component A having a flat spectrum compared to the steep-spectrum component B. The spectra and the fractional variability of ∼30%/yr indicate that this is a core-jet system.

6.2.39. J2337-0622 (2334-066)

This source (Fig. A.32) has a relatively bright galaxy (18.35 mag at r-band) as optical counterpart in SDSS with photometric redshift z = 0.18  ±  0.02. There is an associated Gaia source with 19.4 G-band magnitude and a photo-centre 36.9 mas offset from the VLBI position. The linear distance between components is ∼20 pc. We discarded the source as a milli-lens due to the different (and variable) spectra, and the high SBR observed in some X-band epochs (41.2 in 2018). The Gaia position also does not support the lensing scenario. Interpreting the source as a core-jet system, with A being the bright, flat-spectrum core, and B a knot in the jet, the Gaia position could be tracing either an optical jet, as it lies broadly in the direction of the radio jet, or the host galaxy itself. The strong variability (≫20 %/yr) of both components also suggest the source to be most likely a core-jet system instead.

6.2.40. J2347-1856 (2344-192)

This source (Fig. 4) was previously classified as a CSO (Taylor & Peck 2003; Sokolovsky et al. 2011). While it is also claimed as a bona fide CSO in Kiehlmann et al. (2024a), it remains the only one out of their 79 bona fide CSOs that does not show compelling morphological evidence of being one, lacking both a clearly identified nucleus or extensions of the two-sided jet. There is an indication of a jet connecting the two components at S-band for some epochs found in Astrogeo, but this still needs to be confirmed by more sensitive observations. We keep the source among the milli-lens candidates for now. We found an optical counterpart in SDSS classified as galaxy with 22.68 mag at r-band and photometric redshift z = 0.22  ±  0.13. Then, the linear distance between components is ∼125 pc. Excluding the first two epochs with uncertain calibration, the fractional variability is within a 7%/yr, so we can confirm this source as a possible CSO. However, at this stage, it also passes all our criteria to be a milli-lens. Studying the spectra in more detail with new observations could constrain this source better. Taylor & Peck (2003) found no polarised flux stronger than 1.5 mJy in this source, which would make it difficult to use polarisation observations to study potential differences between the two components. We note, however, that the lack of polarisation makes it less likely to be two lensed images of a radio nucleus, and more likely to be a CSO.

6.3. Further considerations

In total, we found nine remaining milli-lens candidates that fulfil all our criteria. For three of them, we did not yet obtain any follow-up observations. We deem it likely that we will reject most or all of the remaining candidates with sufficient follow-up, especially with high-cadence observations, which will reveal proper motion and variability. Even in the unlikely case, that all nine sources were to be confirmed as milli-lenses, we cannot make statistical statements about the population of SMCOs in this study. This is because the parent sample from Astrogeo, that we are using here, while largely unbiased, does not have redshifts available for all the sources. The lensing probability depends crucially on the redshift distribution of the parent sample studied (Loudas et al. 2022). For the SMILE sample, we are able to make predictions based on the redshift distribution and expect to find a number of lenses ranging from zero up to the order of a few, depending on the DM model studied (Loudas et al. 2022). In that sense, some milli-lenses could indeed be expected to be found in our study here, but no statistical conclusions can be drawn from that.

So far, we have not made predictions about other SMCOs acting as lenses, like PBHs or the centres of globular clusters. This study, as well as a calculation of the expected number of lensing events derived from simulations, is currently being conducted.

As was to be expected, many of the objects in our sample are either CSOs or CSO candidates. In total, we found 34 viable CSO candidates whose data are consistent with the CSO definition at this point. We note that confirming these objects as CSOs will require analysis of more VLBI data and monitoring of their variability, which we leave to a future work (Kumar et al., in prep.). It is interesting to note that none of these objects, which have a spectroscopic or photometric redshift, has a size larger than 500 pc, in good agreement with the findings of Kiehlmann et al. (2024b) that there is a cut-off in CSO sizes at around 500 pc. Since our sample selected from Astrogeo is not complete it cannot be combined with that of Kiehlmann et al. (2024b) to make an improved estimate of the statistical significance of this size cut-off. As pointed out by de la Parra et al. (2024), quantitative comparisons between samples should only be made using complete samples that do not use a spectral index filter, for example. Figures B.1 and B.2 show the spectra for all 40 initial milli-lens candidates using the total fluxes from all available VLBI images used in this work. It is important to note that these data come from different VLBI arrays, and no correction has been made for different resolution or (u, v) range in these observations. Thus, care needs to be taken when interpreting these spectra, and they serve illustrative purposes only. We see primarily steep spectra, with two sources (J1143+1834 and J0052+1633) showing a peak around gigahertz frequencies, as expected for GPS sources, that are likely to be CSOs as well. While J1143+1834 has already been classified as a CSO in previous works, for J0052+1633, it serves as additional new evidence for the CSO nature. A detailed investigation of the spectra of all sources will be dedicated to future work (Kumar et al., in prep.).

7. Summary and conclusions

We have presented an investigation of the nature of 40 milli-lens candidates based on a number of VLBI images from the Astrogeo database and sensitive follow-up observations with the EVN at 5 and 22 GHz. Based on a set of criteria, we ruled out 31 of these candidates as milli-lens systems. Hence, we showed the feasibility of ruling out most milli-lens candidates with observational follow-up. The rejected sources may still be interesting science targets, especially with regard to adding more objects to the class of CSOs, and our preliminary classification here will help future studies.

The remaining nine objects that could not be ruled out will be further scrutinised. We applied for high-frequency observations with the Korean VLBI network (Project ID: KVN-24A-392; PI: Jae-Young Kim) to test the detectability of 20 out of the 40 original candidates and will compute their spectra if detected. This and further monitoring at several frequencies with high-sensitivity observations will prove crucial to ruling out or confirming the remainder of the candidates as milli-lenses. Out of the nine remaining candidates, only one (J1143+1834) was detected in our 22.2 GHz observations with the EVN. A multi-epoch, multi-frequency observing proposal is planned for the VLBA, for which all nine sources will be visible, in contrast to the EVN. This will put tight constraints on spectra and flux variability, distinguishing core-jet systems, CSOs, and possible milli-lenses. Our classification of sources as CSOs or SMBHBs is necessarily incomplete at this stage since, for example, we only ruled out CSOs if their fractional variability is much over 20 %/yr in our data. In addition to new VLBI observations, we will investigate archival data from single-dish and connected-interferometer observations to obtain integrated spectra and constrain the flux variability of sources in much more detail (Álvarez-Ortega et al., in prep.). Future VLBI observatories, such as the ngVLA Long Baseline Array (ngVLA-LONG) and Square Kilometer Array (SKA)-VLBI, are expected to help in distinguishing milli-lens systems by providing very sensitive observations that can also constrain the emission at compact and large scales at the same time. This will also enable the study of even larger samples when looking for potential milli-lens systems.

This pilot study is a crucial step in the SMILE project, which aims to discriminate DM models based on an investigation of ∼5000 sources from the VLA CLASS survey. This flux-complete sample of flat-spectrum sources will largely have redshifts available, thus allowing for accurate predictions of the DM model dependent number of observed milli-lenses, as previously presented in Loudas et al. (2022). It is important to note here that our search for milli-lenses does not require a complete sample as a basis. Rather, what is needed is a sample with known redshifts that gives as many lines of sight as possible to constrain the number density of milli-lenses. We refer to such a sample as a ‘sufficient sample’, meaning a sample that is both clearly defined and unbiased for the purpose of the particular study in question, which in our case is the search for milli-lens systems. (See de la Parra et al. 2024 for more details on the choosing of appropriate samples.)

Data availability

The appendices containing VLBI images and tables of the remaining sources, as well as radio spectra and optical images can be found on Zenodo.

Acknowledgments

We thank the anonymous referee for the very useful suggestions and comments on the manuscript. F.P., C.C., D.A., A.K., V.M. & D.B. acknowledge support from the European Research Council (ERC) under the Horizon ERC Grants 2021 programme under grant agreement No. 101040021. We thank Rocco Lico for the useful discussions. We would also like to thank Leonid Petrov for maintaining the Astrogeo VLBI FITS image database and A. Bertarini, L. Vega Garcia, N. Corey, Y. Cui, L. Gurvits, X. He, T. Koryukova, Y. Y. Kovalev, S.-S. Lee, R. Lico, E. Liuzzo, A. Marscher, S. Jorstad, C. Marvin, D. Homan, J.-L. Gomez, M. Lister, A. Pushkarev, E. Ros, T. Savolainen, L. Petrov, A. Popkov, K. Sokolovski, A. Tao, G. Taylor, A. de Witt, M. Xu, and B. Zhang for making VLBI images they produced publicly available. We found that the practice of uploading VLBI images to a publicly available database brings great benefits to the scientific community. The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code(s): EC071. 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). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, the National Science Foundation, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. SDSS telescopes are located at Apache Point Observatory, funded by the Astrophysical Research Consortium and operated by New Mexico State University, and at Las Campanas Observatory, operated by the Carnegie Institution for Science. The SDSS web site is http://www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including Caltech, The Carnegie Institution for Science, Chilean National Time Allocation Committee (CNTAC) ratified researchers, The Flatiron Institute, the Gotham Participation Group, Harvard University, Heidelberg University, The Johns Hopkins University, L’Ecole polytechnique fédérale de Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania State University, Smithsonian Astrophysical Observatory, Space Telescope Science Institute (STScI), the Stellar Astrophysics Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Illinois at Urbana-Champaign, University of Toronto, University of Utah, University of Virginia, Yale University, and Yunnan University.

References

All Tables

All Figures

	Fig. 1. Two examples of (u, v) coverages for the EVN observations at C-band. Left: J2209+6442 (high declination target). Right: J2347-1856 (low declination target).
In the text

Fig. 2. J2209+6442, a source that has both components detected at 22 GHz and has been rejected as a milli-lens but is retained as a CSO candidate. Left: EVN 4.9 GHz data. The contours start at four times the image rms noise of 0.061 mJy/beam and increase by factors of two. The restoring beam size is 1.65 × 4.24 mas at −31.7° PA. Right: EVN 22.2 GHz data. The contours start at four times the image rms noise of 0.118 mJy/beam and increase by factors of two. The restoring beam size is 0.23 × 0.81 mas at 55.2° PA.

In the text

Fig. 3. J0616-1957 EVN 4.9 GHz data. The contours start at four times the image rms noise of 0.142 mJy/beam and increase by factors of two. The restoring beam size is 1.57 × 3.74 mas at 66.9° PA. The source is not detected at 22.2 GHz. Despite being rejected as a milli-lens candidate, it shows clear morphological evidence of being a CSO with its extended radio lobes, which were only revealed through our new EVN observations.

In the text

	Fig. 4. J2347-1856 EVN 4.9 GHz data. The source is not yet rejected as a milli-lens but has strong indications to be a CSO. The contours start at four times the image rms noise of 0.296 mJy/beam and increase by factors of two. The restoring beam size is 1.55 × 4.44 mas at 68.91° PA. The source is not detected at 22.2 GHz.
In the text

	Fig. 5. Summary charts for milli-lens rejection criteria. Left: Venn diagram showing which combination of criteria lead to source rejection. Right: Bar diagram showing the total number of sources rejected for each criterion alone.
In the text