© The Authors 2024

1. Introduction

Supermassive black hole (SMBH) pairs at sub-kiloparsec scales (∼0.5″) are expected to form at the center of galaxies during the hierarchical assembly of structures(e.g., Begelman et al. 1980; Mayer et al. 2007). During galaxy mergers, SMBH pairs can undergo orbital decay via dynamical friction, form a gravitationally bound binary [e.g.,][among many others]Kelley17,Volonteri21, and eventually coalesce with the associated emission of gravitational waves (GWs), as was predicted by Albert Einstein’s general theory of relativity (Einstein 1936). Dual active galactic nuclei (AGN) are systems containing two active nuclei powered by accretion onto two different SMBHs that are nested inside their host, but that are not yet mutually gravitationally bound. These systems are the natural precursors of coalescing binary SMBHs, which are strong emitters of low-frequency GWs (Colpi 2014). Quantifying their overall abundance would help us to estimate the GW event rate in the future ESA mission Laser Interferometer Space Antenna (LISA), which will unveil the rich population of SMBHs of ∼10^4 − 7 M_⊙ forming in collisions out to redshifts as large as z ∼ 20 (e.g., Amaro-Seoane et al. 2023; Arzoumanian et al. 2018). Moreover, dual AGN are relevant to understanding the theories of hierarchical structure formation, growth, and demography of SMBHs, AGN fuelling and feedback, and the role of merging in AGN evolution (Satyapal et al. 2014; Perna et al. 2023; Sedda et al. 2023). The discovery of two nearby AGN at the same redshift can reveal the presence of a dual system (two different SMBHs), but can also be due to the lensing of a single quasar (QSO) by an intervening foreground galaxy. Lensed systems at small projected separations are a powerful tool with which to address key questions in astrophysics, like studying intrinsically faint and distant lensed sources, investigating dark matter properties of the deflector (Massey et al. 2010), and providing constraints on fundamental cosmological quantities, like the Hubble constant, the cosmological constant, and the density parameter of the Universe (Wong et al. 2020; Lemon et al. 2024). These studies are possible since many properties of the lensed systems depend on the age, scale, and overall geometry of the Universe. Therefore, both dual and lensed AGN are very powerful tools with which to answer a wide range of open questions.

The direct observation of dual and lensed AGN at such small separations is challenging due to the combination of their rarity and the limited spatial resolutions of current wide-field facilities. The ESA Gaia satellite (Gaia Collaboration 2016) has been used in various ways to select these elusive systems. For example, dual and lensed systems have been discovered by looking for multiple Gaia sources associated, within a few arcsecond, with known AGN (Lemon et al. 2017; Chen et al. 2022; Shen et al. 2023). In addition, source variability can also be exploited to find multiple AGN, as was proposed by Kochanek et al. (2006) and employed by Krone-Martins et al. (2019), among others. More compact systems have been identified by looking for the astrometric jitter associated with AGN variability in unresolved dual systems (“varstrometry”; Hwang et al. 2020; Shen et al. 2019; Gross et al. 2023). Despite all these efforts, until recently, very few dual AGN with separations below ∼7 kpc have been confirmed at z > 0.5 (e.g., Lemon et al. 2019; Chen 2021; Chen et al. 2023b; Glikman et al. 2023, see Mannucci et al. 2023 for a full list), making any quantitative testing of cosmological predictions impossible, especially at high redshift.

Recently, a novel technique dubbed “Gaia multi-peak” (GMP; Mannucci et al. 2022) was developed by our group to obtain large and reliable samples of dual and lensed AGN candidates with sub-arcsecond separations. This method looks for AGN that present multiple peaks in the 1D light profiles observed by Gaia. The GMP method consists of selecting AGN with large values of the parameter ipd_frac_multi_peak (FMP) provided in the Gaia catalog since the Early Data Release 3 (EDR3, Gaia Collaboration 2021). This parameter gives the fraction of differently oriented scans in which the object appears to have multiple light peaks. The statistical properties of the GMP selection have been described in Mannucci et al. (2022) and Mannucci et al. (2023).

GMP-selected objects could be physically associated dual AGN, gravitationally lensed systems, or an AGN projected close to a foreground, Galactic star. Therefore, accurate, spatially resolved spectroscopy is necessary to unveil the nature of each system. Integral field spectroscopy is particularly well-suited to extract spatially resolved spectra of each component of these systems because, in most cases, the positions of the various components are not known a priori. Moreover, either space observations or the use of adaptive optics (AO) from the ground are needed due to the small angular separations of the components. The first spatially resolved spectra of GMP-selected systems, obtained with Keck/OSIRIS (Larkin et al. 2006) and HST/STIS, were presented by Mannucci et al. (2022, 2023) and Ciurlo et al. (2023).

Our aim is to discover and classify more GMP dual-AGN candidates. Here, we present AO-assisted observations with MUSE (Multi Unit Spectroscopic Explorer, Bacon et al. 2010) narrow field mode (NFM) at ESO VLT (Very Large Telescope), of 12 GMP-selected systems at z > 0.3, the largest such sample to date. We also present the spectral-decomposition analysis that we applied to spatially unresolved spectra to identify the AGN plus star systems, and to remove them from the target list of the AO observations.

The manuscript is structured as follows. In Sect. 2, we present the target selection and how an accurate analysis of the integrated spectra can help removing the AGN-star systems. In Sects. 3 and 4, we describe the observations and the steps taken to analyze the spectra and identify the nature of each component. In Sect. 5, we analyze and classify each system. Finally, in Sect. 6 we summarize our findings and discuss their implications in the larger context of observed dual AGN at small separations, comparing them with other results from the literature.

2. Target selection

We selected 12 targets observable from VLT from a starting sample of AGN with FMP > 8, to ensure the presence of a double source (Mannucci et al. 2022).

Among these systems, we selected targets with the following requirements: a declination accessible from Paranal (Dec < +15°); a far-enough distance from the galactic plane (|b|> 15°) to reduce stellar contamination; and a sufficient brightness in the near-IR to be used as tip-tilt (TT) star for the MUSE AO system (J < 17.5). This threshold was chosen to be one magnitude brighter than the limit given in the MUSE manual for P109, to allow for AO observations in turbulence class up to 50%. J band magnitudes were either obtained from 2MASS, or estimated from the observed Gaia Bp and Rp magnitudes with an uncertainty of ∼0.27, as described in the appendix. Moreover, we chose targets far from large galaxies in the Local Group (mainly M31, Large and Small Magellanic Clouds) in order to decrease stellar contamination. Finally, we preferentially selected objects with z > 0.5 to probe the redshift range at which most of the merging activity was expected (Chen et al. 2023a), and to avoid the presence of bright host galaxies that could affect the GMP selection (Mannucci et al. 2022).

Nine of the targets were selected among spectroscopically confirmed AGN, ; that is, objects for which previous ground-based, seeing-limited spectroscopy is available (Sect. 2.1).

We complemented this sample with 3 additional photometric GMP-selected quasars (see Sect. 2.2) from the AGN Gaia catalog for which only low-resolution spectroscopy (R ∼ 20–100) from Gaia is available.

The observations were carried out from July 2022 to March 2023. The main properties of these 12 targets are reported in Table 1.

2.1. Spectroscopic targets

We chose nine sources among ∼300 GMP-selected targets with spectroscopic redshift from the Milliquas v7.2 catalog (Flesch 2021) and from our own spectroscopic follow-up (Program ID:109.22W4, PI: Mannucci) obtained with the ESO Faint Object Spectrograph and Camera (v.2, EFOSC2) at the New Technology Telescope (NTT). The results of this spectroscopic campaign will be presented in a future publication.

The spectra of part of the sample were examined to exclude objects showing the presence of a star in their ground-based unresolved spectrum. To achieve this, we developed the novel procedure described in the next section.

2.1.1. Decomposition of ground-based spectra

To identify the systems resulting from chance alignments between AGN and Galactic stars, we developed a procedure to examine the total combined spectra of unresolved multiple GMP-selected sources before the MUSE follow ups. This method enabled us to eliminate from the target list those systems showing significant stellar features at velocity ∼0 km/s with respect to Earth, indication that the companion of the AGN is a star (see also Shen et al. 2023).

Ground-based spectroscopic data (see Table 1) covering both components of the pairs were primarily used to detect the presence of an AGN and to measure its redshift. In GMP-selected systems, the luminosity ratio between primary and secondary component is usually less than a factor of six for all the objects with G > 18.5 (Mannucci et al. 2023). Both components are expected to contribute significantly to the observed spectrum. As a consequence, we proceeded by decomposing the total ground-based spectrum into a combination of a quasar and a star with an innovative routine (see Appendix A in Mannucci et al. 2023 for the first application of this routine).

For the dust-obscured quasar component, we adopted the parametric spectral energy distribution (SED) model from Temple et al. (2021) to reproduce the average observed quasar spectra across the redshift range 0 < z < 5.

For the stellar component, we used a library consisting of a combination of the MaStar (MaNGA Stellar Library) survey from SDSS-IV (see Yan et al. 2019), covering the spectral range 3622 − 10 354 Å at a resolving power of R ∼ 1800. We only considered stars of spectral types G, K, and M, as hotter stars are very rare at high Galactic altitudes. To fit the observed spectra, we considered all the possible combinations of 30 stellar spectra and 6 emission-line QSO templates (see Appendix B in Temple et al. 2021 for details). Different QSO templates span from the minimum to the maximum of the range of observed QSO SEDs, accounting for different levels of ionization of emission-line regions. The decomposition routine degrades the spectral resolution of the template spectra to be the same of the observed data and finds the linear combination of stellar and QSO templates that provides the best-fit to the observed spectrum.

The free parameters of our decomposition procedure include the stellar spectral class, the amount of dust extinction in the QSO spectrum, A_V, E(B − V), and the redshift of the QSO, z_QSO, if not previously known. Other free parameters include the flux normalization of the two components, F_star and F_QSO, and the radial motion of the stellar companion, allowing for a maximum shift of the stellar spectrum of ±300 km/s¹ For each set of parameters we test every possible combination of stellar plus QSO templates, performing a loop over the parameter space until obtaining the set of free parameters that minimizes the residuals. Before comparing the modeled and observed spectra, we accounted for the spectral resolution of ground-based observations by degrading the templates resolution. This is necessary to properly compare emission and absorption features which play a key role in this routine. Since in many cases we observe the co-presence of strong QSO emission lines and weak stellar absorption lines, we determine whether to include a stellar component by estimating the χ² over the main stellar absorption features for each spectral type.

With all these ingredients, we ran the spectral decomposition and classified a system as an alignment of an AGN and a star if the addition of a stellar component decreases the χ² computed over the main absorption lines by at least 30%. While fitting the data with a combination of QSO and stellar templates, and in order to reduce degeneracies between different spectra combination, we also imposed that the magnitudes of the primary and secondary systems obtained from the fit be within the sensitivity limits of the GMP method (see Fig. 2 in Mannucci et al. 2023). This decomposition requires a S/N and a resolution of the spectra sufficient to identify and characterize the features of the star. We expect this decomposition routine to be very effective in revealing the presence of a star if it has similar apparent luminosity of the primary QSO. If the luminosity ratio is close to the largest value allowed by the GMP selection (∼2 mag for targets with G < 18.5), this procedure could be unable to detect the presence of a star.

We developed this optimisation routine in November 2022, and used it since then to examine spectroscopic data of GMP sources for the presence of stellar features and define the MUSE target sample. At that date, six spectroscopic targets had been already observed (see Observation date in Table 1). Consequently, we applied this new decomposition method to the targets J0015–2526, J1059+0622, and J0052–3045. As discussed in Sect. 6, this procedure turned out to be very effective in identifying and removing AGN–star systems before spatially resolved observations. As a consequence, the chance to observe dual and lensed systems with IFU observations drastically increases.

2.2. Non-spectroscopic Gaia AGN candidates

We complemented the nine spectroscopic targets described above with 3 systems from the Gaia AGN catalog selected according to the procedure described below. We selected objects with values of FMP above the adopted threshold of 8 (see Sect. 2.1) among those classified as AGN candidates based on the Discrete Source Classifier (DSC). DSC is a module that classifies Gaia sources probabilistically into 5 classes: quasar, galaxy, star, white dwarf, and physical binary star. QSO candidates are identified by high values of the probability parameter classprob_dsc_combmod_quasar, having values between 0 and 1. We selected objects classified as very probable QSO (classprob_dsc_combmod_quasar > 0.97). This threshold value was chosen as a compromise between reliability of the classification and the number of observable targets. We also required high values of two additional parameters provided by the Quasar Classifier (QSOC, Delchambre et al. 2023) module. The first parameter, Zqsoc, indicates the reliability of the redshift estimated by Gaia, while the processing binary flag, flags_qsoc, reports the potential errors that can occur during the prediction process from QSOC. We considered sources with Zqsoc above the threshold of 0.99 (out of a maximum of 1) in order to obtain optimal redshift estimates and with the processing binary flag flags_qsoc = 0, indicating the absence of any warning flags.

3. Observations and data reduction

All the objects were observed with VLT/MUSE, which provides optical, integral-field unit (IFU) spectra (programs ID: 109.22W5 and 110.23SM, PI: Mannucci). This paper includes data acquired until March, 2023. The main properties of the 12 observed targets are reported in Table 1. We used the 7.5″ × 7.5″ field of view (FoV) NFM of the instrument, assisted by the AO system (GALACSI, Ströbele et al. 2012) employing 4 laser guide stars (LGS, Lewis et al. 2013). The spatially resolved spectra obtained have a spatial sampling of 0.025″ × 0.025″ and cover a spectral range from 4750 Å to 9350 Å with a resolving power ranging from 1700 in the blue to 3400 in the red. In NFM-AO mode the range 5780–6050 Å is blocked to avoid contamination by sodium light from the LGS (Na Notch filter). We employed slightly different observing strategies in P109 and P110. For the former semester, we used 4 MUSE position angles (PAs) rotated by 90 deg from each other, for a total exposure time of 35 min per object. Having verified that 3 PAs are enough to obtain a good sky subtraction and remove the instrument artifacts, in P110 we reduced the overheads by decreasing the PAs to 3, for a total exposure time between 19 min and 27 min per target, depending on the magnitude of each of them. To ensure the highest quality observations and data accuracy, we require seeing full width at half maximum (FWHM) better than 0.8″ and a fraction of lunar illumination less than 50%.

We reduced our data by using the processing pipeline for the MUSE Instrument (Weilbacher et al. 2020) running under the EsoReflex (ESO Recipe Flexible Execution Workbench, Freudling et al. 2013) graphical environment. We used the standard workflow, subtracting the sky by extracting it from 50% faintest parts of the field of view, away from the target, which takes up only a small portion of the FoV.

In all cases, the target itself is used as reference star for the AO system. We didn’t find problems in closing the AO loop on these double targets if the total magnitude is significantly brighter than the nominal MUSE limit of J ∼ 18, depending on seeing conditions. The AO system always provided an excellent correction, and we measure FWHMs of the point-spread-function (PSF) at ∼ 7000 Å between 0.14″ and 0.17″. J0052–3045 exhibits the worst PSF, with a FWHM ≃ 0.17″ in the direction along the separation. Despite this, since the separation among the components is 0.44″ (see Tab. 2), the two sources are well distinguished even in this case.

4. Data analysis

MUSE data confirm that the GMP selection is able to successfully recover multiple systems: all the 12 observed systems are resolved into two point-sources with separations between 0.30″ and 0.86″, within the expected range according to Mannucci et al. (2022). In Figure 1, we show the cut-outs extracted from the data cubes around the position of the targets, whilst in Figure 2 we present the spatially resolved spectra extracted from the squared apertures shown on the images. In all cases, the separations are much larger than the PSF FWHM, therefore all the spectra are independent and can be used to separately study the nature of each source and measure its properties. The dimensions of the squared aperture could influence the data by introducing corrections related to the PSF. However, these effects remain consistent for both components. Hence, this does not impact our analysis. Note that, given the limited Strehl ratios (SR) of our observations, our spectra only contain part of the total flux of the objects, the same fraction for the two components of each system. Our results are not affected by this effect because they are all based on relative quantities.

Fig. 1. White light images (target name and separation of the two objects are reported on the figure). The spectra are extracted from the square apertures plotted over the maps. Orange squares enclose the targets later identified as stars, while the confirmed AGN are displayed with violet and light blue apertures (A and B components).

We modeled each spectrum with either an AGN or a stellar spectrum. For AGN, we used the templates created by Vanden Berk et al. (2001) and by Temple et al. (2021). To account for the effects of interstellar dust extinction, we considered 3 different extinction laws to the AGN templates: Calzetti et al. (2000), Cardelli et al. (1989), and an empirical quasar extinction law presented in Temple et al. (2021). We fit all the observed spectra with these AGN templates, by varying E(B − V) and the source redshift. For each dual and lensed AGN, the best-fit values and classification are presented in Table 2. Specifically, according to the quality of the fit, we used the Cardelli et al. (1989) extinction law to correct both spectra of J0402-3227, while for all other spectra we applied the empirical extinction law established by Temple et al. (2021). On the other hand, Galactic extinction has a limited impact due to the high Galactic latitudes of the targets (|b|> 21°), and therefore is negligible and ignored.

For stars, we used the same library as detailed in Sect. 2.1, considering stars of G, K, and M spectral types. These stars have much redder colors than AGN and, in most cases, their spectral features dominate the red part of the spectrum. Stars hotter than spectral type G are exceedingly rare far from the galactic plane; thus we exclude them. Our fitting procedure tests all the possible stars with different normalizations and selects the best one with the χ² method. We report all the results in Table 3.

For the systems classified as composed by two AGN, our goal is to assess similarities and differences among the spectral features of the two images and determine whether they are associated to the same AGN lensed by a foreground galaxy, or to two distinct quasars at small separations. In principle, the spectra of the two images of a lensed system are expected to be identical, except for a normalization factor, since gravitational lensing is an achromatic phenomenon. In practice, intrinsic variability, microlensing, and differential dust extinction can introduce differences in the observed spectra and mimic the presence of a dual system.

1. Variability can introduce differences between the spectra when coupled with the different time delays between the two images. In our case, this is not expected to be a problem. In fact, the time delay between two images at sub-arcsecond separations is a few days at most (Wambsganss 1998; Lemon et al. 2024), while at the luminosities probed in this work (L_bol = 10^{46 − 47} erg s⁻¹), the BLR is expected to extend for hundreds of light days (Bentz et al. 2013; Grier et al. 2019; Prince et al. 2023), and the narrow-lines to be emitted over at least hundreds of parsecs. As a result, no significant variations in line profile or line ratios is expected due to this effect. However, small differences are expected in the continuum and, thus, in the emission broad line equivalent width (EW).

2. While lensing is acromatic, microlensing, caused by compact objects (for instance, stars) within the lens galaxy, can produce chromatic effects. Microlenses differentially magnify parts of the quasar emission regions, leading to time- and wavelength-dependent changes in the flux ratios of the images (Wambsganss 2006). Microlenses generally act on source regions at scales that are comparable to or smaller than a few micro-arcsec; that is, the size of their angular Einstein radius. Consequently, it is possible that both the quasar’s continuum region, such as the accretion disc, and the broad line region will experience micro-lensing effects. Given that magnification varies in relation to source size, the power-law continuum emission of quasars will exhibit greater magnification at bluer wavelengths, emitted by more compact regions (e.g., Mosquera et al. 2013; Jiménez-Vicente et al. 2015). This also applies to the emission lines with higher ionization, which are expected to be subject to microlensing effects larger than those of the lower-ionization ones. Line shapes can also be affected, since different parts of the BLR could be differentially magnified/de-magnified. In contrast, the narrow lines, emitted over regions of hundreds or thousands of pc, are not subject to any microlensing effect. Therefore, we have to define if the differences in the continuum shape, the discrepancies in line equivalent width and the deformations in the line profiles can be due to the possible combination of differential extinction and/or microlensing within the lensing galaxy or not. The analysis of this phenomenon is reported in Sect. 4.4.

3. The light paths corresponding to the different images of a lensed system pass through different regions of the lens, and can intercept different amounts of gas and dust. As a consequence, the observed spectra could be characterized by different amount of dust extinction and different EWs of the absorption lines associated with the lens. This is addressed in Sect.4.1.

Absorption lines also provide information on the nature of the system. If one lensed image exhibits one or more narrow absorption lines (NALs) at the quasar’s redshift (intrinsic NALs), then the other lensed image should also display the same feature. This is because the intrinsic NALs reflect the gas composition within the quasar. However, differences in the EW of NALs may arise in the case of lensing, either due to microlensing or due to the geometry of the lens, with the line of sight of each image crossing different regions of the lens (e.g., Chartas 2007; Green 2006). In contrast, dual AGN can present different lines shapes, both for narrow and broad lines, and can easily have different absorption features, since they are separated by few kpcs.

In summary, to assess the nature of all the systems in our study, we jointly analyze the emitted spectra of both AGN, focusing on the following properties:

• Redshift

• EWs of the detected lines, and line profile.

• Emission line flux ratios

• Presence of narrow absorption lines

• Continuum level and shape

If all these quantities are consistent within the uncertainty in the two spectra, we classify the system as a lensed QSO candidate. Otherwise, significant differences lead us to identify them as dual AGN. The number of quantities that show significant differences between the sources is directly proportional to our confidence in the dual nature of the system.

4.1. Cross-correlation

In order to automatize the classification of the system (lensed or dual) through the comparison of the two spectra, we developed a cross-correlation approach. This method employs 3 distinct sets of free parameters, adapted to the characteristics of the examined system through χ² minimization.

The basic concept of our approach is to test if two AGN spectra can be produced by a lensed AGN. We thus modify one of the two AGN spectra in all the possible ways that are compatible with the effects produced by lensing, with the aim of matching the spectrum of the other AGN, applying our modifications step-by-step as follows:

1. Initially, we cross-correlate the two spectra using two free parameters only, namely, a shift in redshift and a magnification factor in flux. In the case where the first parameter significantly deviates from zero, i.e., the two components show different redshifts, we classify the system as dual. In the opposite case, if the shift is compatible with zero, if a perfect match (within the errors) is obtained, the systems is classified as lensed candidate, otherwise we proceed applying other modifications.

2. As a first refinement to the previous approach, we introduce an additional free parameter that accounts for the different normalization between the continuum and the lines that the two spectra display. This is possible in a lensed system, because the continuum can vary on shorter timescale as a consequence of the AGN intrinsic variability or of microlensing effects, while the narrow lines are expected to be stable. This is the case for J0015–2526 for example, where we find that the spectra of the two components can be matched by using two different normalization factors for the continuum and the emission lines, differing by 67%. Consequently, this system is classified as a lens.

3. As a second refinement, we multiply the observed spectra for a second degree polynomial in order to account for possible differences in the extinction law between the two light paths. This can happen as the images that we observe come from the same source but pass through different regions of the lensing galaxy. We proceed with this further correction in the case of J0525–4153, where the two sources have different values of E(B − V), 0.155 and 0.014 respectively.

4.2. Line-fitting procedure

Another crucial step to classify the systems, distinguishing dual from lensed AGN, entails the analysis of the spectral lines. To analyze the emission lines, we fit multiple Gaussian functions with different centers, amplitudes, and widths to each line in all observed sources. Gaussian fitting with two components (i.e., one narrow and one broad) was employed for broad lines, in addition to a power law for the continuum, while only one Gaussian was used for the narrow lines. For all the lines of each object, we measure the center and we quantify their EW. Central wavelengths and EW values are reported in Table 4.

4.3. Presence of the lens

For all systems constituted by two AGN, we further explored the possibility of identifying the lensing galaxy in our datacubes, as expected in case of a lensed system.

First, we examined the broad-band image obtained by collapsing the cube along the wavelength axis to look for the presence of a foreground galaxy. To increase the image sensitivity, we also produced separate images in 4 different wavelength ranges (5300, 6500, 7100, 8700)Å devoid of bright sky. No lens galaxy is directly detected in any of these images of the systems analyzed in this work.

Secondly, to identify potential faint or dust-obscured galaxies hidden by the presence of bright images of the lensed QSO, we subtracted the contribution of the AGN and examined the residuals. We fitted the PSFs by summing two 2D elliptical Gaussian functions, and proceed to subtract the contribution of the two AGN components. In all our systems, the observed residuals can be attributed to an imperfect fit to the PSF, and no lensing galaxy is unambiguously detected.

A lens non-detection is not a smoking gun to discard the lensed AGN nature if the lens can be fainter than the detection limit. Thus, we investigate whether the luminosity upper limits of the undetected lens can be used to exclude its presence and therefore secure the dual AGN nature.

We can roughly estimate the expected luminosity of the lensing galaxy as a function of its unknown redshift. We use half of the observed separation as a proxy for the Einstein radius and calculate the mass of the lens galaxy inside it as a function of redshift (M_⋆(z)) assuming a point-like mass distribution (e.g., eq. (21) in Narayan & Bartelmann 1996). We assume that the lenses are early-type galaxies; in other words, the faintest systems for a given mass, as they have the oldest stellar populations. Assuming that stars constitute the majority of the mass inside the Einstein radius and that they follow a Sersic profile with an index of n = 4, by using the mass-size relation in van der Wel et al. (2014) we determine the needed stellar mass as a function of redshift. From this, we use Longhetti & Saracco (2009) to estimate the K- and V-band magnitudes of the possible lensing galaxy as a function of redshift. In every scenario, the possible lenses are brighter if either at low redshift (low mass, but nearby galaxy) or high redshift (close to the source, distant but very massive galaxies), and fainter for intermediate redshifts.

The resulting galaxies are then convolved with the PSF and added to the datacubes and residual maps.

In all cases except J0402–3237, the lenses are expected to be below the detection threshold at any redshift. In J0402–3237, characterized by the largest separation of 0.67″ and therefore requiring the largest masses, the lens should be detected for any redshift, and its absence can aid us in classifying the system as dual rather than lensed.

In conclusion, for all the systems except J0402–3237, the absence of a visible lens galaxy in MUSE data does not preclude the presence of a lensing effect. The sensitivity of all spectroscopic data may not be sufficient to reveal extended objects with low surface brightness and no bright emission lines, such as the lensing galaxy. This limitation persists even with the subtraction of the AGN PSF. The intricate structure of the PSF, influenced by atmospheric conditions, remains the primary factor hindering a clear delineation of the lens.

4.4. Gravitational microlensing

Several techniques can be used to unveil microlensing in lensed quasars. In this work, we applied the so-called Macro-micro decomposition (MmD) from Sluse et al. (2007), which is based on comparing the differences in the (microlensed) continuum and broad lines to those in the (un-affected) narrow lines, once the macrolensing magnification is already accounted for. This method enables an easy visualization of the differential microlensing affecting the QSO spectrum.

As stated by Sluse et al. (2007) we assume that the spectrum of an observed lensed image, denoted as F_i(λ), can be decomposed into two distinct components. The first component, F_M, exclusively accounts for macro-lensing effects, while the second, F_Mμ, incorporates both macro and micro-lensing influences. We also defined the macro-magnification ratio between image 2 and image 1 as M = M₂/M₁(> 0), the relative micro-lensing factor between the images as μ, and A is the scaling factor between the F₂ and F₁ continua.

The behaviors of μ and M exhibit distinct characteristics. The factor M may vary with wavelength due to differential extinction in the lensing galaxy. On the other hand, μ is presumed to remain constant within the limited wavelength range covered by a line profile, but chromatic microlensing may occur due to the wavelength dependence of the source size; the blue continuum-emitting region may be smaller and more microlensed than the red one.

Hence, we calculated μ and M in close proximity to each emission line for every system.

To definitively assess the presence or absence of the microlensing effect, one needs to consider two lensed images and determine both components, F_M and F_Mμ, using the following equations:

$$\begin{array}{cc}\hfill F_1& =M\times F_M+M\times \mu \times F_{M\mu }\hfill \\ \hfill F_2& =F_M+F_{M\mu }\hfill \end{array}$$

Up to a scaling factor, F_M only depends on A, while F_Mμ depends solely on M. The A parameter can be accurately determined as the value for which F_M(A) = 0 in the continuum adjacent to the line profile. Assuming A remains constant along the line profile, F_M exclusively reflects the contribution from the emission profile. If we assume that a part of the observed emission experiences only macro-amplification (without micro-amplification), we can estimate the macro-amplification factor M. Given the relation A = M × μ, we can then estimate the micro-amplification factor μ for the continuum at the wavelength of the line profile. On the other hand, if the emission line is amplified exactly like the continuum, microlensing cannot be distinguished from macrolensing.

5. Results

The details of the resulting classification are listed in Table 3 for the 5 AGN/star associations, and in Table 2 for the seven dual AGN or lensed QSO candidates. Following the cross-correlation procedure described in Sect. 4, we classify 3 objects as lensed AGN candidates and 4 as dual systems.

We classify J0015–2526, J0525–4153, and J1059+0622 as lensed AGN candidates since their components exhibit overlapping spectra, with similar trends in continuum and EW for all the lines.

On the other hand, J0052–3045 is a clear dual system because its spectra show a significant difference in the lines on the NLR and in redshift, estimated to be Δz = 0.010 (∼3000 km/s) through cross correlation analysis. However, this value could be smaller due to the absence of narrow lines in one of the AGN (light blue in Fig. 2) making it difficult to accurately estimate the redshift. As a result, its measurement relies solely on the broad narrow line MgII, which exhibits a distinctly different profile between the two components. In the other cases (J0402–3237, J0840–0529, and J2228–3047), the spectra reveal different EWs and line ratios. In J0402–3237 and J2228–3047, the continuum also differs. Additionally, in the latter object, the largest difference is observed in the EW of the iron lines, specifically in the region around 7000 Å.

Fig. 2. Spectra of the systems observed with MUSE (target name and redshift are reported in the panels). They have been extracted over the square apertures marked on the white light images in Figure 1 (with the same color-coding). Vertical dotted blue lines mark the position of the main expected emission lines. For J0052−3045, emission lines are plotted in light blue and violet as the redshift for the two AGN is not the same. The star in J0007−0053 has been multiplied by ten to optimize the visualization. The gap around 6000 Å is due to the New Generation Controllers (NGC) deployed in adaptive optics.

We note that some spectra reveal the presence of NALs from intervening material. Here, we only discuss the absorption features of J0052-3045, which are used for the classification of the system. The other absorption systems will be described in a future work.

In the following, we describe the properties of each of the observed systems and the procedure followed for their final classification.

5.1. AGN+star

5 of the observed systems are associations of an AGN and a foreground star (see Table 3). Figure 2 presents the MUSE spectra of these 5 systems, where the orange line shows the spectrum of the star. The classification of the secondary system as a star is performed through the technique presented in the previous section; that is, by finding the best-fitting stellar spectrum. Two of these systems, J0007+0053 and J2139–4331, have previously been observed and studied by Lemon et al. (2020). In that work, J2139–4331 was already classified as an AGN+star system, a classification confirmed by our spectra, while no classification was provided for J0007+0053 due to the limited (seeing-limited) spatial resolution.

Moreover, it is noteworthy that among these systems, the separation of J2139–4331 measured by MUSE data (0.82″) differs significantly from the Gaia measurement (0.666″). This disparity can be attributed to the high proper motion of the star measured by Gaia.

5.2. Lensed AGN systems

Both the components of the seven remaining systems show AGN spectra. We start by describing the results of the MUSE spectral analysis for the 3 systems classified as lensed AGN candidates. The values of all the spectral features described below are reported in Table 4.

5.2.1. J0015−2526

The classification of this system as a lensed AGN candidate is based on the striking similarity of the two spectra, and confirmed through the cross-correlation discussed in Sect. 4.1 and shown in Fig. 3. For J0015−2526, we quantified a magnification factor for the emission lines (Fact_lines) and applied a correction to this value to reproduce the continuum (Fact_cont). The use of two different values is needed because the continuum could be affected by microlensing and could show variability on a short timescale (e.g., Fian et al. 2018). The correction factor obtained, 67%, is in agreement with this scenario. The values of these two scaling factors are reported in Fig. 3.

Fig. 3. Cross correlation J0015−2526_B (light blue AGN spectra): multiplied in flux and shifted to the same redshift by the quantities indicated in the labels. Color-coding as in Fig. 2.

In the spectra of J0015−2526 there are six detected lines, of which MgIIλ2798, Hδλ4103, and Hγλ4342 were fitted following the procedure explained in Sect. 4.2. Gaussian fitting with two components (i.e., one narrow and one broad) was employed for MgII, in addition to a power law for the continuum, while only one Gaussian was considered to model the Balmer lines due to the lower S/N. Central wavelength and EW values were computed for each line, and were found to be consistent between the two AGN. Additionally, the two sources have the same emission line flux ratios within the errors, providing strong evidence that this is indeed a lensed system. Notably, the H_δ/H_γ ratio, which is determined from two narrow lines, is perfectly equal in the two components.

Since J0015−2526 is classified as a lensed system, we tried to identify the lens as we explained in Sect. 4.3.

First, we produced 4 images in 4 wavelength ranges devoid of bright sky lines, i.e.: λ = (5300, 6500, 7100, 8700) Å. We summed these images to search for any possible lensing galaxy in the field. The lens galaxy is not detected, in either the summed image or the partial images.

Second, we investigated the possible properties of a lensing galaxy to test if its non-detection in the MUSE datacube is consistent with the lensed AGN nature of J0015−2526. Based on the measured quasar redshift (z = 0.970). We tested this for 3 lens redshifts: 0.25, 0.50, 0.75. For each redshift value, as explained in 4.3, we derived the mass enclosed by the Einstein ring, and used van der Wel et al. (2014) and Longhetti & Saracco (2009) to assess the size and magnitude. We then simulated the 3 galaxies assuming a Sersic profile with n = 4 and we add each galaxy to our data. With these magnitudes (G > 21.8, 22.5, 22.0, respectively), none of the tested lensing galaxies is brighter than the limiting sensitivity of our data (G ∼ 20, estimated by adding compact galaxies of various brightness to the data). Furthermore, we observe some absorption lines in the spectra of J0015 that are indicative of an intervening absorber at z = 0.75. The lens corresponding to this redshift could be a typical elliptical with mass M = 10¹¹ M_⊙, and an effective radius R_eff = 0.4″, with a magnitude of G = 22.4, thus undetectable in our data.

5.2.2. J0525−4153

The AGN J0525−4153, at z = 2.155, is the first observed object selected from the Gaia catalog without a previous spectroscopic confirmation. The cross correlation, obtained by also considering a second degree polynomial for a supplementary difference in the extinction, clearly shows a close match between the spectra (Fig. 4). Continuum and lines are perfectly overlapping, as well as the several narrow absorption lines, strongly suggesting that this is a lensed AGN.

Fig. 4. Cross correlation J0525−4153_B (light blue AGN spectra): multiplied in flux and shifted to the same redshift by the quantities indicated in the labels. Color-coding as in Fig. 2.

To study the emission lines, we primarily focus on the MgII and CIV emission line components, which we analyzed through the Gaussian fitting procedure. For the MgII emission line, we employed two Gaussian profiles and a power-law component, while in the case of the CIV line we made a fit with a Gaussian and a power-law profile only. We estimated the equivalent widths of both lines, which were found to be compatible, yielding mutually consistent line ratios (MgII/CIV). Moreover, we observed that CIV is blue-shifted by the same velocity in both spectra, providing further confirmation of the lensed nature of this system.

Additionally, we identified the presence of intervening NALs at the same observed wavelength in both spectra, suggesting that they originate from the same gas along the line of sight of the quasar. We identify the observed lines as FeIIλ2382, FeIIλ2585, FeIIλ2600, MgIIλ2796 and MgIIλ2803 at z = 1.597 and z = 1.889. These intervening NALs could be due to the lens galaxy located between the two AGN images, or they could be associated with other objects along the line of sight.. Further investigation of the NALs could yield even more detailed information about the lens galaxy and other intervening objects, their redshift and chemical abundances. However, as in the previous case, our data are not deep enough to detect the lens in a system with this configuration.

In summary, our analysis of the spectra of J0525−4153 has provided strong evidence supporting its identification as a lensed AGN candidate.

5.2.3. J1059+0622

J1059+0622 was previously analyzed by Lemon et al. (2017) who classified it as a lensed candidate based on variability.

The spectra of J1059+0622 exhibit the presence of 3 emission lines, although one of them (CIII]) is partially lost due to the LGS gap (see Fig. 2). Since the spectra seem to be apparently similar as is the case for gravitationally lensed sources, we tried all the cross-correlation types to test this scenario. However, none of the tested spectral modifications produce perfectly overlapping spectra and a clear difference in CIV line is always present.

We fit MgII and CIV lines with a Gaussian component plus a power-law for the continuum (dotted violet and light blue lines in Fig. 5). The estimated EWs and line ratio MgII/CIV for the two components differ significantly. We report the obtained values in Table 4. Since these differences can be caused by microlensing, we tested the presence of this effect with the MmD procedure of Sluse et al. (2007) (described in Sect. 4.4). Applying it to both the observed spectra and the dereddened ones, we identified a plausible scenario involving microlensing of the continuum along with differential extinction, potentially attributable to a lensing galaxy. The derived MmD values for the CIV and MgII in the observed spectra are as follows: A = 1.0,  0.53, M = 0.38,  0.32, and μ = 2.63,  1.66, respectively. The difference of the value of M for the two lines is compatible with the differential extinction Δ(E(B − V)) = 0.73 reported in Table 2, thus these values are perfectly compatible with differential extinction plus microlensing of the continuum (Sluse et al. 2012). The same interpretation holds for the dereddened spectra, for which we obtain the following MmD values for the CIV and MgII lines: A = 0.3,  0.11, M = 0.26,  0.16, and μ = 2.73,  1.62 (Fig. 6). Also in the dereddened spectra, the weak deformation of MgII can be due to residual differential extinction. This might be related to the empirical extinction law applied during the dereddening process Temple et al. (2021), Sect. 4), however, we evaluated the robustness of the decomposition with the “pre-processing” of the spectra by applying the decomposition on the dereddened spectra. The discrepancy in E(B − V) could be ascribed to differential reddening of the two AGN images due to the two distinct lines of sight that probe different regions of the lensing galaxy. Therefore, the spectra are consistent with coming from two images of a single, lensed source.

Fig. 5. Fit to the MgII (top) and CIV (bottom) lines of both components of J1059+0622. Color-coding of Fig. 2. Each component of the fit is shown with dotted violet and light blue line, for components A and B, respectively.

Fig. 6. MmD decomposition following the procedure in Sluse et al. (2007) on the dereddened spectra of J1059+0622. The upper panel shows the CIV line, while the lower panel displays the MgII line. The resulting values of the relevant parameters are reported in the figure

Another indicative feature that prompts consideration of a lensed system is the presence of intervening NALs at a redshift of 1.19 in both spectra. These may potentially be generated by the lens at this particular redshift. However, there is no hint for residual contamination by a potential lensing galaxy. Applying the procedure outlined in section 4.3 similarly to what done for J0015, a typical lensing galaxy at z = 1.19 would not be detected in our data and neither would be if placed at a different redshift, given that the limiting magnitude of the MUSE observation is G ∼ 20.5. Thus, a lens non-detection is not a good enough reason to discard the lensing scenario for this source.

In summary, the most plausible classification for this system is a lensed AGN. Variations in the equivalent widths, preventing a perfect cross-correlation, are likely a result of microlensing and differential extinction.

5.3. Dual AGN systems

In the following, we present the results of our MUSE data analysis of the 4 systems classified as reliable dual AGN.

5.3.1. J0052−3045

J0052−3045 stands out among the analyzed AGN as it exhibits clear differences in redshift, line position, and line profiles, indicating that this is a genuine dual system (see Fig. 2). The MgII line is present in both spectra, with a different shape and a velocity shift of approximately ≃1000 km/s, corresponding to redshifts of z = 1.455 and z = 1.445 for the two components (see Table 4). This shift cannot be attributed to variability in the line profile since the expected time delay between different lensed images is a few days at most. Furthermore, shape and line ratios of narrow emission lines (i.e. [NeV] and [OII]) are very different. These differences cannot be easily explained by the variability of the quasar, hence they provide an additional evidence of the distinct nature of the two AGN.

Moreover, as we can see in Fig. 7, there are large differences in the visible intrinsic absorption lines in the two distinct spectra. The A component (violet) shows the MgII doublet in absorption at rest-frame wavelengths of 2796 Å and 2803 Å, ; in other words, at the same redshift of the QSO. Since these NALs are entirely absent in the second AGN, gas composition of the two AGN is different, providing evidence that they cannot be two images of the same source. The presence of these lines in only one of the two AGN further supports the hypothesis of a dual system, as one AGN has an outflow that is not visible in the other. According to all our findings, J0052−3045 is a dual AGN.

Fig. 7. Enlargement of the MgII doublet absorption lines of J0052−3045. The lines are found at rest-frame wavelengths of 2796 Å and 2803 Å. Color code is the same of Fig. 2.

5.3.2. J0402−3237

The system in question was previously observed by Lemon et al. (2020), who were not able to provide a reliable classification due to the limited spatial resolution of their data. In fact, our MUSE cubes fully resolve the individual components and allow us to classify the system as a dual AGN. We de-reddened the spectra of J0402−3237 by using the MW attenuation, and noted that the values of E(B − V) are markedly different between the two components, namely 0.189 and 0.861. Notably, below 5600 Å the B component (light blue in Fig. 2) is very faint, resulting in a very low signal-to-noise ratio, thus we exclude this portion of the spectrum in the comparison. We tried to improve the match between the two spectra by incorporating a second-degree polynomial to account for potential variations in extinction between the spectra.

Our analysis of the line ratios revealed numerous differences, shown in Fig. 8 and quantified in Tab. 4. Specifically, the EWs of MgIIλ2798, [NeV]λ3426, [NeIII]λ3869 and [NeIII]λ3968 lines are found to be larger in the fainter AGN (i.e., component B, in light blue), while [NeV]λ3346 and [OII]λ3728 are more prominent in the brighter one (i.e., the violet one).

Fig. 8. Details of the spectra of J0402-3237, with the same color-coding of Fig. 5. Left panel: fit to the MgII line. Right panel: part of the spectrum showing several narrow lines, after multiplying the spectrum of the B component by the factor indicated in the label. Note the clear difference in the observed [OII] emission line and in the shape of the continuum.

Assuming it is a lensed system, we checked for the presence of a lens or the effect of microlensing explained in Sects. 4.3 and 4.4.

J0402−3237 is the only case in which the lens is expected to be visible in the white light image. For any possible lens redshift, we conducted simulations of lensing galaxies, yielding magnitudes below G < 21.5 (at z ∼ 0.7). All of these simulated lenses are detectable under the specific system configuration of z_QSO = 1.28 and sep = 0.666″. However, despite these efforts, we were unable to detect it, further strengthening our initial interpretation as a true dual AGN.

Using the MmD method of Sluse et al. (2007) on MgIIλ2798 and [NeV]λ3426 lines, we find that the microlensing event(s) that could explain the discrepancies between the spectra has a rather strong amplitude, but still not unfeasible. If the observed chromaticity is interpreted as microlensing, we would require M to be the same for MgII and [NeV] (M = 0.4). Anyway, this is not supported by the data because it implies μ = 1.65 (> 1) at the level of [NeV], but μ = 0.37 (< 1) at the level of MgII. This implies a full microlensing of the [NeV] line, which is not possible. Hence, the observed chromaticity of the spectral ratio as (only) chromatic microlensing is not viable. The alternative is considering the possibility that the chromaticity in J0402−3237 results by differential extinction. In this case the quantity that varies chromatically is M (M_MgII = 0.4 and M_NeV = 1.4) while μ remains constant (μ = 0.37). We find that microlensing of the full [NeV] line would be necessary also in this scenario, but it is unfeasible as narrow lines are not affected by microlensing (Sluse et al. 2007). Thanks to all the differences described above, the lack of lensing galaxy and the microlensing analysis, we can classify this system as a dual AGN (Table 4).

5.3.3. J0840−0529

J0840−0529 can also be reliably classified as a dual AGN since, despite having the same redshift, the two spectra exhibit remarkably distinguishable characteristics. For example, the shape of the MgII line is different in the two spectra, as it appears to be narrower in the fainter AGN (component B, light blue one in Figure 9).

Fig. 9. MgII line of J0840−0529’s components, color-coded as in Fig. 2. To optimize the visualization, the B component has been multiplied by the factors in the label.

Also the cross-correlation analysis clearly shows significant differences in the CII and MgII emission lines, i.e.; the line flux ratio between MgII and CII is significantly different. In addition, we calculated the EW of MgII for both the AGN, following the analysis explained in Sect. 4.2 and fitting the line with two gaussian and one power-law. It turned out to be significantly different, 24 ± 3 for component A vs. 41 ± 4 for component B.

As detailed in section 4.4, we employed the MmD method of Sluse et al. (2007) to examine whether the variations in line shape, EWs, and line ratios could be attributed to microlensing. It is not possible to get a plausible decomposition of MgII. The profile deformations are so pronounced that only the case of microlensing acting stronger on the lines than on the continuum could explain the observed properties. Such an event very unlikely, since microlensing is expected to be stronger on the continuum and less evident on the lines. Consequently, the case of a dual scenario with two distinct AGN is more plausible.

5.3.4. J2228−3047

Both AGN of J2228−3047 display prominent emission lines of MgII and CIII]. Several characteristics suggest that this is a dual system. We fit the CIII] line with one Gaussian component plus a power law for the continuum, whereas we added a second Gaussian for the broad-line contribution in the fit of the MgII line (Fig. 10). The fits indicate that the difference in the line central wavelengths in the two AGN is too large to be due to variability of a single source over a few days. We also calculated the EWs of the narrow lines of MgII and CIII] and their line ratio. Values are significantly different between the two components, see Table 4. Furthermore, the cross-correlation analysis also shows intrinsic differences in the EWs of the iron lines.

Fig. 10. Upper panel: Cross correlation J2228−3047_B (light blue AGN spectra): multiplied in flux and shifted to the same redshift by the quantities indicated in the labels. Lower panel: Fit to the MgII line of both components of J2228-3047. Color-coding as was in Figure 8.

Still, this system could be a lensed AGN if microlensing affects the continuum solely through flux magnification and the MgII line through deformation of its shape. We try to systematically apply the MmD (Sluse et al. 2007) of J0840–0529 to see if we can get a plausible decomposition of it. Data are rather noisy, and consequently, the decomposition is uncertain. At the beginning, the data seem to be compatible with a lensing hypothesis implying a microlensing of the continuum μ = 0.7 and a small line deformation of MgII (but not of the CIII] complex). This hypothesis necessitates an additional minor differential extinction or intrinsic variability to account for the observed differences (around 20%) between CIII] (M = 0.4) and MgII (M = 0.32). On the other hand, the absence of a consistent decrease or increase in the flux ratio between the blue and red wings of these lines suggests to us to exclude the lensing hypothesis. This may be harder to explain if μ is effectively achromatic as was suggested by the 2 lines analyzed. Consequently, we cannot completely dismiss the microlensing hypothesis, particularly given the low S/N.

Nonetheless, the results of the analysis of J2228−3047 further suggest the possibility of the duality of this system.

6. Discussion and conclusions

We used spatially resolved, AO-assisted spectroscopy with MUSE-NFM to unveil the nature of 12 multiple AGN candidates selected through the GMP method (Mannucci et al. 2022); that is showing multiple peaks in the Gaia light profile, at redshifts between 0.3 and 2.9. Nine of the targets have been selected from objects with previous spatially unresolved spectroscopy (Sect. 2.1), and 3 from the Gaia QSO candidate catalog (Sect. 2.2). GMP-selected objects could be either, physically associated dual AGN, or, gravitationally lensed systems, or an AGN projected close to a foreground Galactic star. Spatially resolved observations allow us to obtain individual spectra for both components of each target and derive a reliable classification of the systems. A novel decomposition routine, based on reproducing the available, spatially unresolved spectrum of the systems with a combination of an AGN and a stellar spectra (see Sect. 2.1) was used to remove targets with clear hints of the presence of a foreground star.

The spatial resolution of our AO spectra, having PSF FWHM ∼ 0.15″, allows us to resolve all the sources. All the targets in our sample present two components, with separations between 0.30″ and 0.86″ (Fig. 1). This confirms that the GMP method is very effective in successfully recovering multiple systems with separations below 7 kpc at z > 0.3.

5 systems turn out to be associations between an AGN and a foreground star. We examine the remaining seven AGN pairs using multiple steps to determine whether the systems are lensed QSOs or dual AGN (Sect. 4), taking into account all the possible sources of differences between the spectra of the different components of a lensed QSO (intrinsic variability, microlensing effects, and differential dust extinctions). For this, we have compared the two spectra by applying different possible distortions due to these effects. We have also compared the position, fluxes, and line profiles of the emission lines, and looked for the possible presence of lens galaxies. The possible effect of microlensing was considered using the MmD decomposition presented in Sluse et al. (2007).

This accurate analysis of the spectra allowed us to unveil the nature of all the seven AGN pairs (Sect. 5). We classified 4 dual AGN, and 3 lensed QSOs candidate, as was outlined in Table 2. Hence, this work significantly increases the number of sub-arcsec separation lensed and dual AGN to date.

As has been reported by Mannucci et al. (2023), to date there are 15 confirmed dual AGN at z > 0.5 with projected separations below 7 kpc, including the 4 new sources classified in this work. Future observations from the ground (especially with VLT/MUSE, VLT/ERIS, and Keck/OSIRIS) and from space (HST/STIS, JWST) are expected to reveal more of these systems, hence largely increasing the number of confirmed dual AGN, as well as sub-arcsecond lensed systems.

This sample is still too small to derive meaningful statistical information, but it is rapidly growing. In the near future, it will be possible to start testing many of the predictions of the models of SMBH merging, such as the SMBH mass function, distribution of separations, fraction of multiple AGN over the general population of quasars, and their evolution with redshift (e.g., Rosas-Guevara et al. 2019; Volonteri et al. 2022; Chen et al. 2023a). Moreover, the compact lensed systems that we are discovering will allow us to sample the inner part of the lensing galaxies and measure the mass of the baryonic component in regions where the dark matter component is negligible (e.g., Smith & Lucey 2013).

Acknowledgments

Based on observations collected at the European Southern Observatory under ESO programs [MUSE programs ID: 109.22W5 and 110.23SM; NTT program ID: 109.22W4] This publication was produced while attending the PhD program in PhD in Space Science and Technology at the University of Trento, Cycle XXXIX, with the support of a scholarship financed by the Ministerial Decree no. 118 of 2nd march 2023, based on the NRRP - funded by the European Union – NextGenerationEU – Mission 4 “Education and Research”, Component 1 “Enhancement of the offer of educational services: from nurseries to universities” – Investment 4.1 “Extension of the number of research doctorates and innovative doctorates for public administration and cultural heritage”. We thank the anonymous referee for his detailed comments to the original manuscript that contributed in improving the text and leading to stronger conclusions. We are deeply grateful to Dominique Sluse for his assistance in computing the analysis of gravitational microlensing. We thank Alessandro Sonnenfeld for useful discussions. We acknowledge financial contribution from INAF Large Grant “Dual and binary supermassive black holes in the multi-messenger era: from galaxy mergers to gravitational waves” (Bando Ricerca Fondamentale INAF 2022), from the INAF project “VLT-MOONS” CRAM 1.05.03.07, from the French National Research Agency (grant ANR-21-CE31-0026, project MBH_waves), from the INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities”. This work has been partially financed by the European Union with the Next Generation EU plan, through PRIN-MUR project “PROMETEUS” (202223XPZM). G.V. acknowledges support from the ANID program FONDECYT Postdoctorado 3200802 and the European Union’s HE ERC Starting Grant No. 101040227-WINGS. 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 (SDSS) has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington.

References

Appendix A: Estimation of QSO near-IR colors from Gaia photometry

The TT star wavefront sensor used by MUSE-NFM can only work on-axis, i.e, using the science target itself as TT star. For this reason it is necessary to select targets that are bright enough at near-IR wavelengths to obtian a good AO correction. In most cases, the QSO are too faint to have measured J- and H-band magnitudes by 2MASS (Skrutskie et al. 2006) or any other survey. For this reason, we have derived an estimate of these magnitudes from a linear combination of the observed Rp and Bp magnitudes from Gaia. We have considered all the QSO with observed J- and H-band magnitudes, mostly from 2MASS, and computed the combination of the two Gaia magnitudes minimizing the difference between the observed and predicted magnitudes. The derived relations, together with their scatter, are:

$$\begin{array}{c}\hfill J=Rp-0.82\ast (Bp-Rp)-0.26,\sigma =0.27\end{array}$$(A.1)

$$\begin{array}{c}\hfill H=Rp-0.89\ast (Bp-Rp)-0.92,\sigma =0.31\end{array}$$(A.2)

These relations are shown in Fig. A.1. The observed scatter is small enough to obtain a reliable value of the expected near-IR mag to

Fig. A.1. Comparison between the observed J- (left panel) and H- (right panel) magnitudes and the values derived from the Bp- and Rp- Gaia magnitudes using equations A.1 and A.2. Each point is color-coded with the observed (Bp-Rp) color. The solid line shows the one-to-one relation, the dotted lines report the observed scatter around this relation.

All Tables

All Figures

	Fig. 1. White light images (target name and separation of the two objects are reported on the figure). The spectra are extracted from the square apertures plotted over the maps. Orange squares enclose the targets later identified as stars, while the confirmed AGN are displayed with violet and light blue apertures (A and B components).
In the text

Fig. 2. Spectra of the systems observed with MUSE (target name and redshift are reported in the panels). They have been extracted over the square apertures marked on the white light images in Figure 1 (with the same color-coding). Vertical dotted blue lines mark the position of the main expected emission lines. For J0052−3045, emission lines are plotted in light blue and violet as the redshift for the two AGN is not the same. The star in J0007−0053 has been multiplied by ten to optimize the visualization. The gap around 6000 Å is due to the New Generation Controllers (NGC) deployed in adaptive optics.

In the text

	Fig. 3. Cross correlation J0015−2526_B (light blue AGN spectra): multiplied in flux and shifted to the same redshift by the quantities indicated in the labels. Color-coding as in Fig. 2.
In the text

	Fig. 4. Cross correlation J0525−4153_B (light blue AGN spectra): multiplied in flux and shifted to the same redshift by the quantities indicated in the labels. Color-coding as in Fig. 2.
In the text

	Fig. 5. Fit to the MgII (top) and CIV (bottom) lines of both components of J1059+0622. Color-coding of Fig. 2. Each component of the fit is shown with dotted violet and light blue line, for components A and B, respectively.
In the text

	Fig. 6. MmD decomposition following the procedure in Sluse et al. (2007) on the dereddened spectra of J1059+0622. The upper panel shows the CIV line, while the lower panel displays the MgII line. The resulting values of the relevant parameters are reported in the figure
In the text

	Fig. 7. Enlargement of the MgII doublet absorption lines of J0052−3045. The lines are found at rest-frame wavelengths of 2796 Å and 2803 Å. Color code is the same of Fig. 2.
In the text

	Fig. 8. Details of the spectra of J0402-3237, with the same color-coding of Fig. 5. Left panel: fit to the MgII line. Right panel: part of the spectrum showing several narrow lines, after multiplying the spectrum of the B component by the factor indicated in the label. Note the clear difference in the observed [OII] emission line and in the shape of the continuum.
In the text

	Fig. 9. MgII line of J0840−0529’s components, color-coded as in Fig. 2. To optimize the visualization, the B component has been multiplied by the factors in the label.
In the text

	Fig. 10. Upper panel: Cross correlation J2228−3047_B (light blue AGN spectra): multiplied in flux and shifted to the same redshift by the quantities indicated in the labels. Lower panel: Fit to the MgII line of both components of J2228-3047. Color-coding as was in Figure 8.
In the text

	Fig. A.1. Comparison between the observed J- (left panel) and H- (right panel) magnitudes and the values derived from the Bp- and Rp- Gaia magnitudes using equations A.1 and A.2. Each point is color-coded with the observed (Bp-Rp) color. The solid line shows the one-to-one relation, the dotted lines report the observed scatter around this relation.
In the text