© The Authors 2026

1. Introduction

According to the Λ Cold Dark Matter (ΛCDM) cosmological paradigm, structure forms hierarchically, with smaller systems such as galaxies merging to build more massive ones (White & Rees 1978; White & Frenk 1991). The most massive galaxies host a central supermassive black hole (SMBH), which grows through gas accretion and mergers and ultimately regulates star formation and shapes the structure of the host galaxy through feedback processes (Di Matteo et al. 2005; Croton et al. 2006).

During galaxy mergers, the SMBHs involved undergo orbital decay via dynamical friction and can form bound pairs (Begelman et al. 1980; Volonteri et al. 2003; Colpi 2014). When the two SMBHs are actively accreting material, the system becomes a dual active galactic nucleus (AGN; Koss et al. 2012; Rosas-Guevara et al. 2019; Volonteri et al. 2022), with projected separations from tens of kiloparsecs (kpc) down to sub-kpc scales (De Rosa et al. 2019; Pfeifle et al. 2025). Dual AGN represent the direct precursors of gravitationally bound binary SMBHs at parsec separations (e.g., Kelley et al. 2017; Volonteri et al. 2022), which may eventually coalesce and emit gravitational waves (GWs). The GW signals from SMBH binaries are expected to be detectable with pulsar timing arrays (Arzoumanian et al. 2018; Agazie et al. 2023; EPTA Collaboration et al. 2023) and, at higher frequencies, with the upcoming Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2023; Colpi et al. 2024).

Dual AGN therefore offer a unique laboratory to investigate gas accretion and feedback in merging galaxies and the connection between black hole growth and hierarchical structure formation. Characterizing their observed separations, accretion rates, BH masses, and host-galaxy properties provide key constraints for models of SMBH coalescence and for the prediction of the GW detection rate of future missions such as LISA, which will be sensitive to merging SMBHs with masses of 10^4 − 8 M_⊙ up to z ∼ 20 (e.g., Amaro-Seoane et al. 2023). Quantifying the abundance and demographics of dual AGN is therefore essential to linking galaxy evolution with the emerging GW landscape (Sedda et al. 2023; Perna et al. 2025).

It is important to clarify the terminology used in this work. In the literature, the term “dual AGN” generally refers to systems hosting two actively accreting SMBHs over a wide range of projected separations, from ∼100 kpc down to sub-kpc scales (e.g., Liu et al. 2011; De Rosa et al. 2019; Pfeifle et al. 2025). The Cosmic Duets systems targeted here represent a specific subset of close dual AGN selected to probe the small projected separation regime (< 7 kpc), where the merging galaxies are expected to be in an advanced stage of interaction and the two SMBHs reside within a common merger remnant. While wider dual AGN systems trace interacting galaxies that are still distinct, Cosmic Duets is focused on the regime most directly connected to the late stages of SMBH orbital decay, merger-driven accretion, and the onset of feedback within a single galactic potential.

The direct identification of dual AGN at approximately kpc separations remains challenging due to their intrinsic rarity, the high angular resolution required to resolve two active nuclei, and the potential increase in obscuration associated with late-stage mergers (Koss et al. 2018; De Rosa et al. 2019). Simulations and observations indicate that column densities and obscuring material tend to increase as mergers progress (e.g., Capelo et al. 2017; Blecha et al. 2018; Chen et al. 2023a; Ricci et al. 2017). Several discoveries have in fact been serendipitous (e.g., Junkkarinen et al. 2001; Glikman et al. 2023), including cases revealed with spatially resolved Integral Field Unit (IFU) spectroscopy using JWST (Übler et al. 2024; Matsuoka et al. 2024; Perna et al. 2025) and candidates identified from complex broad-line structures in spatially integrated spectra (Maiolino et al. 2024), highlighting the difficulty of systematically uncovering such systems.

Initial targeted searches relied on indirect or resource-intensive techniques, such as the selection of double-peaked emission lines (Zhou et al. 2004; Comerford et al. 2009; Rubinur et al. 2019). Other systematic approaches have targeted dual AGN through spectroscopic galaxy pairs (Liu et al. 2011; Jing et al. 2025), radio-selected samples (Fu et al. 2015a,b, 2018), mid-infrared preselection (Satyapal et al. 2017; Pfeifle et al. 2019), quasar clustering studies (Hennawi et al. 2006, 2010), probabilistic mid-IR pairs (Barrows 2023), and serendipitous discoveries in lens surveys (Inada et al. 2008, 2012). However, many of these surveys primarily probe systems at projected separations of several tens of kiloparsec, corresponding to earlier merger stages or physically associated pairs rather than the compact (∼few kpc) regime where both SMBHs are embedded within a common merger remnant. While a subset of studies (e.g., Fu et al. 2015a, 2018; Satyapal et al. 2017; Pfeifle et al. 2019) have focused on AGN pairs with projected separations < 10 kpc, the majority of candidates identified across these surveys exhibit projected separations on scales of tens of kiloparsecs, largely driven by large low-redshift samples such as Liu et al. (2011). In addition, these works predominantly probe redshift regimes that are distinct from those targeted here and are based on selection techniques sensitive to different AGN populations (e.g., radio-selected or narrow-line AGN in Fu et al. 2015a, 2018 and mid-infrared-selected AGN in Satyapal et al. 2017; Pfeifle et al. 2019). As a result, they often explore a different region of parameter space in separation, dynamical state, redshift, and gas obscuration compared to the sub-arcsecond dual AGN targeted here.

A major breakthrough came with the Gaia mission (Gaia Collaboration 2016) of the European Space Agency (ESA), which enabled systematic all-sky searches for closely separated AGN through multiple complementary techniques. Dual AGN candidates can be identified as multiple Gaia detections within a few arcseconds (Lemon et al. 2017; Chen et al. 2022; Shen et al. 2023; Lemon et al. 2024), via characteristic variability (Kochanek et al. 2006; Krone-Martins et al. 2019), or through astrometric jitter in the photocenters of Gaia quasars (“vastrometry”; Shen et al. 2019; Hwang et al. 2020; Chen et al. 2022; Schwartzman et al. 2024, 2025; Gross et al. 2025). High-resolution imaging from ESA/Euclid provides a powerful complement, enabling the identification of closely separated AGN candidates at sub-arcsec scales. Machine-learning methods applied to Euclid images allow for systematic searches for dual photometric AGN in order to uncover sources that are often too faint for Gaia detection (Ulivi et al. 2025).

These developments have enabled recent large-scale efforts to expand the census of close dual AGN, particularly at z > 0.5 and projected separations below 1″ (a few kiloparsecs at cosmic noon). An example is the compilation presented by Pfeifle et al. (2025), which brings together candidate and confirmed systems identified with heterogeneous selection techniques, although the number of spectroscopically confirmed dual AGN at these separations remains limited. This parameter space, in which simulations place SMBHs at separations consistent with a common merger remnant (Volonteri et al. 2022), has been targeted by innovative Gaia-based techniques. In particular, the innovative Gaia multipeak (GMP) method identifies sub-arcsecond multiple AGN candidates from Gaia light-profile peaks (Mannucci et al. 2022). Typically, GMP candidates are brighter than Euclid-selected sources, which reach a depth of I_AB ∼ 24.5 in the VIS filter, whereas Gaia detects sources down to G ∼ 20.5 mag (Vega-like), enabling efficient IFU follow-up to spectroscopically distinguish dual AGN, lensed quasars, or AGN–star projections (Mannucci et al. 2022, 2023; Ciurlo et al. 2023; Scialpi et al. 2024). However, as an optical selection technique, GMP is primarily sensitive to unobscured, optically bright AGN and is therefore likely biased against heavily obscured dual AGN, particularly at higher redshift, a limitation shared by most current high-z dual AGN searches.

A fraction of GMP-selected sources are expected to be gravitationally lensed quasars by foreground galaxies, producing multiple images at sub-arcsec separations. Although rare, these compact lensed systems are powerful probes of the central regions of lens galaxies, where stellar mass dominates (Treu 2010). In these regions, they enable precise measurements of the mass enclosed within the Einstein radius (total mass), stellar mass-to-light ratios (M_★/L), and the stellar initial mass function (Treu et al. 2010; Shajib et al. 2024), while also constraining the relative contributions of baryons and dark matter and the overall mass profile of the lens galaxy (Sonnenfeld et al. 2015; Massey et al. 2010; Newman et al. 2017). Moreover, strongly lensed quasars provide independent constraints on cosmological parameters, such as the Hubble constant (Wong et al. 2020; Lemon et al. 2024), and have been extensively used in the literature to test models of structure formation and lens galaxy evolution (e.g., Kochanek et al. 2006; Treu & Ellis 2015).

Beyond their role in probing lens galaxies and cosmology, strongly lensed AGN are also powerful tools to study the quasars themselves. The magnification and amplification produced by lensing acts as a natural “boost” of the AGN signal, allowing detailed studies of high-redshift quasars that would otherwise be too faint to observe. Lensing can also induce microlensing variability, which probes the internal structure of the AGN, including the accretion disk and broad-line region (e.g., Sluse et al. 2012; Mediavilla & Jiménez-Vicente 2021). By combining lens modeling with high-resolution imaging and spectroscopy, lensed AGN provide a unique laboratory to investigate the physical properties of quasars across cosmic time.

Several large surveys have systematically targeted lensed AGN to build statistical samples, providing insights into both lens galaxies and background quasars. The Sloan Digital Sky Survey (SDSS) Quasar Lens Search (Inada et al. 2005, 2012) produced a total of about 62 strongly lensed quasars in the SDSS Data Release 7 (DR7) quasar catalog, 26 in the well-defined statistical sample, and 36 additional lenses identified with various techniques, providing one of the largest homogeneous samples for statistical studies of strong lensing. The Sloan Lens ACS Survey (Bolton et al. 2008) used spectroscopic selection from SDSS to find 131 strong galaxy-galaxy lenses at z ∼ 0.05 − 0.5, providing detailed constraints on the mass, light, and kinematics of massive early-type galaxies. The SL2S survey (Sonnenfeld et al. 2013, 2015) and the COSMOS HST survey (Faure et al. 2008) identified tens of strong lens candidates, allowing the study of stellar-to-dark matter fractions, radial mass profiles, and magnification properties in a larger but still moderate redshift range. The Grism Lens-Amplified Survey from Space (Treu & Ellis 2015) extended lens spectroscopy to cluster fields at z ∼ 0.3 − 0.7, delivering near-IR spectra and emission-line redshifts for both lenses and background sources, with a focus on faint systems.

While these surveys provide valuable statistical samples, they typically target lenses with image separations ≥1″ and redshifts below two, limiting the study of very compact and higher-redshift systems. In contrast, the GMP method (Mannucci et al. 2022) identifies subarcsec multiple AGN candidates, including some of the most tightly separated lensed systems known, with quasar (QSO) redshifts in the range 0.5 < z < 3.5 and projected separations between 0.15″ and 0.8″. These compact lenses enable high-resolution imaging and spectroscopic follow-up, allowing us to extend detailed analyses to a population of lensed AGN previously inaccessible to systematic statistical studies, while simultaneously probing both the structure of the lens galaxy and the properties of the AGN. Moreover, our technique targets Einstein radii smaller than the effective radius of the lens galaxies, enabling us to probe their inner regions even at high redshift (Sonnenfeld et al. 2019; D’Amato et al. 2026).

Here we present the European Southern Observatory (ESO) Large Program 114.27BY (PI: M. Scialpi) consisting of adaptive-optics (AO) assisted Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010) narrow-field mode (NFM) observations at ESO’s Very Large Telescope (VLT) targeting 150 GMP-selected systems at sub-arcsec separations and z > 0.5. In this paper, we describe the dataset of 30 systems collected during the first year (2025) as well as the observing strategy, data reduction, and calibration procedures, and we outline the analysis used to classify each system as dual AGN, lensed AGN, or foreground stellar alignment.

The paper is structured as follows. In Sect. 2 we describe the GMP-based target selection and the use of integrated spectra to filter out AGN+star contaminants. Section 3 presents the MUSE observations, the data reduction, and the resulting dataset of images and spectra. Section 4 details the spectral modeling and analysis procedures used to classify the systems and study their properties, while Sect. 5 reports the spectroscopic properties and final classification of each observed system. Section 6 examines the distributions of dual and lensed AGN as a function of separation, magnitude, and redshift while taking into account selection effects. Finally, in Sect. 7 we summarize the results, discuss the implications for the population of dual and lensed AGN at subarcsec separations, and outline the prospects for the full MUSE Large Program. All magnitudes we report are in the Vega system. We adopted a flat ΛCDM cosmology with parameters from Planck Collaboration VI (2020), with H₀ = 67.7 km s⁻¹ Mpc⁻¹ and Ω_m, 0 = 0.31.

2. Target selection

We selected our AGN candidates using the GMP technique (Mannucci et al. 2022, 2023) applied to Gaia DR3 light profiles (Gaia Collaboration 2021), as briefly summarized in the following, while full details can be found in Mannucci et al. (2022). The fraction of Gaia scans displaying a double peak is quantified by the FMP parameter (Gaia Collaboration 2021, 2023), ranging from 0 (never detected) to 100 (always detected). We adopted a selection threshold of FMP > 8, shown by Mannucci et al. (2022, 2023) to be a reliable criterion for identifying multiple systems.

The primary source must be brighter than G < 20.5, corresponding to the efficiency limit of the GMP method (Mannucci et al. 2023). The effective selection window for GMP dual-AGN candidates is therefore ∼0.15″ − 0.8″, as detectability depends on both angular separation and luminosity ratio. Systems at separations smaller than the Gaia point spread function (PSF) (∼0.15″) or with large luminosity contrasts (L_A/L_B ≳ 15) may fail to produce a measurable multipeak signature, while sources with separations ≳0.8″ are typically resolved as distinct detections in the Gaia archive. This upper limit naturally selects dual AGN that are physically close enough to reside within the same host galaxy, while also ensuring that the GMP method is still sensitive to multiplicity. In Table 1 we list the properties of all systems, and for those resolved by Gaia we report the individual magnitudes of both components (G1 for the primary and G2 for the secondary).

To minimize stellar contamination, we restricted the selection to regions at Galactic latitudes |b|> 20° and excluded areas surrounding large Local Group galaxies (M31, LMC, SMC). We also excluded over-crowded regions of the sky, defined as having more than 12 Gaia detections per square arcminute within 3 arcmin of the target. Priority was given to systems at z > 0.5, both because galaxy merger activity peaks in this redshift range (Chen et al. 2023a) and because at lower redshift the brightness of the host galaxy can significantly reduce the GMP detection efficiency (Mannucci et al. 2022). This ensures an optimal balance between astrophysical relevance and selection completeness.

The GMP selection was applied to a parent sample of AGN drawn from two complementary datasets: a catalog of spectroscopically confirmed AGN and a catalog of photometric AGN candidates, which are described below.

(i) Spectroscopic AGN. This sample includes sources with available spectra and secure redshift measurements, drawn from large spectroscopic surveys such as the Sloan Digital Sky Survey (SDSS; Lyke et al. 2020) and the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2025). For these sources, we adopt version 8 of the Milliquas compilation (Flesch 2023).

(ii) Photometric AGN. This sample consists of AGN candidates selected from recent photometric catalogs (Delchambre et al. 2023; Storey-Fisher et al. 2024; Flesch 2021; Kunsági-Máté et al. 2022; Fu et al. 2024; Salvato et al. 2025), which lack spectroscopic confirmation but significantly extend the accessible parameter space.

The spectroscopic sample provides robust AGN classifications and accurate redshift measurements for the primary component up to z ∼ 4. In contrast, the photometric catalogs extend the GMP selection to regions not yet covered by large spectroscopic surveys, increasing the sample size at the expense of lower purity.

Together, these datasets enable an efficient identification of multiple-peak AGN candidates across a broad range of magnitudes and redshifts, approaching full-sky coverage.

2.1. Spectral confirmation of photometric candidates

To confirm the AGN nature of photometric candidates and obtain robust spectroscopic redshifts, we carried out follow-up observations through several programs using multiple instruments. Specifically, we employed the ESO Faint Object Spectrograph and Camera (v.2, EFOSC2) mounted at the Nasmyth B focus of the 3.58 m ESO New Technology Telescope (NTT; Buzzoni et al. 1984; Program IDs: 109.22W4, 112.25CT, PI: Mannucci; 113.26TB, PI: Scialpi), the Device Optimized for the LOw RESolution (DOLORES) on the Telescopio Nazionale Galileo (TNG; Program IDs: A47TAC_28, A48TAC_67, PI: Scialpi), and the FOcal Reducer/low dispersion Spectrograph 2 (FORS2) on the VLT in Long-Slit Spectroscopy (LSS) mode (Program ID: 115.28CQ, PI: Scialpi). These observations allowed us to confirm the AGN classification of the primary sources and to determine their spectroscopic redshifts. The systems in which the presence of an AGN as the primary source is confirmed are further analyzed to reliably exclude the presence of secondary stellar components, thereby removing systems that correspond to AGN–star pairs through the analysis described in Sect. 2.2. Observations were conducted under varying atmospheric conditions, with a typical seeing better than 1.1″. A 1″ slit was used by default and widened when necessary to ensure that both components remained within the slit. Exposure times were optimized using the ESO and TNG Exposure Time Calculators to achieve a continuum S/N of ∼10 − 15 per spectral resolution element for targets with G ≃ 20. Each exposure sequence included additional overheads for acquisition and calibration frames, resulting in total integrations of 20−60 min per target depending on the instrument and conditions.

The redshifts of all targets observed with MUSE are listed in Table 1, while the full catalog will be presented in a future paper (Scialpi & Marconcini, in prep.).

2.2. Spectral decomposition

To remove stellar contaminants as secondary components in our sample, we built a Python routine that quantitatively estimates the likelihood of contamination via spectral decomposition of the unresolved spectra (see Sect. 2.1). The spectral decomposition used here is un updated version of the one described in Scialpi et al. (2024) and will be fully presented in Marconcini & Scialpi (in prep.). We applied our custom spectral decomposition to all targets with available integrated spectra, including AGN from DESI and SDSS and photometric candidates confirmed through our follow-up observations, for a total of 650 sources. All these spectra lack the spatial resolution to disentangle the single components, which may correspond to a second AGN (either an independent AGN or multiple images of the same lensed AGN) or a foreground aligned star. This decomposition process has proven crucial to remove a large fraction of AGN+star contaminants prior to IFU follow-up observations (for a similar approach see Shen et al. 2023).

In brief, the observed spatially unresolved spectrum is fit using a Python implementation of the Penalized PiXel-Fitting algorithm (pPXF, Cappellari & Emsellem 2004) with a linear combination of quasar (from a library covering different ionization levels; Temple et al. 2021) and stellar (G, K, M types) spectral templates from the MaStar library (Yan et al. 2019), covering the wavelength range 3622 − 10354 Å at R ∼ 1800.

First, we convolve the quasar and stellar templates with the data spectral resolution. Then, we create a single template grid composed of 26 stellar spectra (G0-G9, K0–K7, M0–M8) and 66 quasar spectra. Quasar spectra are created using the qsosed package, with intrinsic visual extinction in the range [0, 1] with spacing 0.1 and spanning different levels of ionization.

Then, we used pPXF to infer the best-fit combination of stellar and quasar templates to reproduce the data. In doing so, the free parameters are the stellar spectral type, the QSO extinction (A_V; E(B − V)), the redshift (z_QSO), and the flux normalization between templates (F_star, F_QSO). Moreover, we allowed for a spectral shift of the stellar component along the line of sight to account for its proper motion of ±300 km/s.

The routine first degrades the spectral resolution of the templates to match the observed data, then finds the linear combination of stellar and QSO templates that best fits the spectrum, accounting for dust extinction, redshift, and radial velocity shifts. While the formal statistical uncertainties on the redshift derived from the pPXF fit are typically smaller than the spectral resolution, we adopt a conservative estimate based on the spectral resolution of the data. Given the resolving power of our observing setups, this corresponds to a redshift uncertainty of Δz ∼ 0.001 over the redshift range of our targets (z ∼ 0.7 − 3.2). This accuracy is sufficient for secure line identification and for planning IFU follow-up observations.

To quantify the probability of stellar contamination in our spectra, we defined seven of the brightest stellar absorption features for each stellar type and evaluated the χ² in a spectral window of ±300 km/s around each feature to estimate the improvement relative to a pure quasar template. Based on this, we assign a probability of stellar contamination when the improvement in χ_R² exceeds 10%, or alternatively when ΔBIC > 6. This procedure removes most AGN+star contaminants prior to IFU follow-up. Mannucci et al. (2022) found that 30% of GMP systems at 0.5″ are expected to be AGN+star, increasing to 60% at 0.8″ (scaling roughly with the square of the separation at small separations; see the discussion in Sect. 5.1).

For this first year, we also targeted some objects that had a non-zero probability of containing a star based on the spectroscopic decomposition, in order to test the method. The preselection procedure removed roughly 75% of the expected stellar contaminants, implying that about 30 AGN+star systems were filtered out prior to the MUSE follow-up. Among the 30 objects observed in year one, only 11 were confirmed to be AGN+star systems (see Table 1), confirming the method’s effectiveness.

2.3. Targets for MUSE follow-up

From the subset of GMP-selected sources with available spectroscopic data (both fiber and single-slit spectra), and after applying all the selection cuts described above, we retained only those targets whose spectral decomposition showed no evidence of a stellar secondary component. We then applied additional constraints to ensure that the selected systems were observable from Paranal (Dec < +15°) and sufficiently bright in the near-infrared to serve as tip-tilt reference stars for the MUSE AO system (J < 17.5). J magnitudes were obtained from 2MASS (Skrutskie et al. 2006), VHS (McMahon et al. 2013), and other surveys. In cases where the objects were too faint to be detected, we estimated their magnitudes from a linear combination of the observed R_p and B_p magnitudes from Gaia (Gaia Collaboration 2021), with a typical uncertainty of ∼0.27 mag (see the appendix of Scialpi et al. 2024):

$$\begin{array}{c}\hfill J=R_p-0.82(B_p-R_p)-0.26,\sigma =0.27.\end{array}$$(1)

Applying this procedure to the full GMP-selected AGN sample, we identified the best 150 targets suitable for MUSE-NFM follow-up observations. These sources span the redshift range 0.5 ≲ z ≲ 3.5, as higher-redshift AGN are rarely detected by Gaia, and exhibit angular separations between ∼0.15″ and ∼0.8″ (corresponding to projected physical separations of approximately 1 − 7 kpc across this redshift range), consistent with close AGN pairs in interacting or merging systems. Simulations suggest that at such separations the SMBHs may often reside within a common merger remnant (e.g., Volonteri et al. 2022; Chen et al. 2023a).

These sources span the redshift range 0.5 ≲ z ≲ 3.5, as higher-redshift AGN are rarely detected by Gaia, and exhibit angular separations between ∼0.15″ and ∼0.8″ (corresponding to projected physical separations of approximately 1 − 7 kpc across this redshift range), consistent with close physical AGN pairs in merging systems (e.g., Volonteri et al. 2022; Chen et al. 2023a).

3. The Cosmic Duets large program

In this work, we discuss the observing strategy, data reduction, and analysis procedures of our Cosmic Duets ESO Large Program 114.27BY (PI: Scialpi), whose ultimate goal is to obtain AO assisted MUSE-NFM observations to characterize dual and lensed AGN at subarcsec separations in the redshift range 0.5 < z < 3.5. Here we present the results for the objects observed during the first year of the program, from 1 October 2024 to 30 September 2025. The program, motivated by the identification of 150 GMP-selected dual and lensed AGN candidates observable with MUSE/NFM (Sect. 2), has a total allocation time of 137 hours. The same observing strategy and analysis procedures described here will be applied to the rest of the sample, which will be observed in cycles P116 and P117.

3.1. Observations

The observations were carried out with MUSE in NFM assisted by the GALACSI AO system at VLT (Bacon et al. 2010; Ströbele et al. 2012; Arsenault et al. 2008). GALACSI operates in closed-loop at 1000 Hz, delivering near-diffraction-limited images. The instrument field of view (FoV) is about 7″ × 7″. The pixel size is 0.025″ × 0.025″, covering the observed-frame wavelength range 4750−9350 Å with a spectral resolving power R ∼ 1700 − 3400, corresponding to a velocity resolution of Δv ∼ 180 − 90 km s⁻¹. The wavelength range 5780−6050 Å is blocked by a sodium notch filter to avoid contamination from the laser guide stars (LGS). It uses four sodium LGS to correct for high-order turbulence and a natural guide star (NGS) for tip-tilt correction. Data were collected over multiple nights from October 2024 to September 2025 under excellent seeing conditions (FWHM < 0.85″) and with lunar illumination < 50%. Each target was observed with three dithered exposures offset by 0.2″ to improve PSF sampling, for a total exposure time of 30−45 min depending on target brightness. This strategy ensures sufficient signal-to-noise (S/N), effective sky subtraction, and mitigation of instrumental artifacts.

In all cases, the target itself was used as the AO reference source. In some observations, closing the AO loop on these double systems proved challenging, as the seeing conditions varied and only the total system magnitude (typically brighter than J ∼ 18) was known, while the individual magnitudes of the two components were not. Apart from these occasional issues, the AO system consistently delivered excellent correction, with measured PSF FWHMs at ∼7000 Å between 0.14″ and 0.17″. The multiple components of all targets are well separated (see Fig. 2) and enabled detailed analysis of each source.

In contrast to the common practice of rotating the field between exposures (e.g. by 90°), we maintained a fixed position angle (PA) for all observations. This choice is driven by the nature of our targets: pairs of point-like sources separated by ≲0.5″ and of comparable brightness, with the brightest used as the NGS for AO correction. Although the AO loop remains closed during individual exposures, rotating the field between them increases the risk that the NGS temporarily falls outside the control region or that the deformable secondary mirror returns to the wrong configuration when the PA is reset, potentially leading to aborted exposures and increased overheads. Moreover, field rotation is not required for our science case. Unlike deep observations of extended, low–surface-brightness objects, where rotation mitigates systematic background structures, our compact point sources allow robust background estimation directly from multiple regions of the FoV. Rotating the field therefore provides no practical advantage while introducing operational risks, while a fixed PA ensures stable AO performance and observational efficiency.

3.2. Data reduction

We reduced our data using the MUSE pipeline (Weilbacher et al. 2020) within the EsoRex (ESO Recipe Execution Tool) environment. The standard reduction involved bias and dark subtraction to remove instrumental noise, flat-fielding to correct for detector response variations, wavelength calibration using arc lamp exposures, characterization of the line spread function to account for spectral resolution, sky subtraction and flux scaling. The sky was estimated from the faintest 50% of the field, avoiding regions containing the science target, which occupies only a small portion of the FoV. Flux calibration was performed using a standard star observed during the same night with the same instrumental setup. Circular aperture extraction ensured that the star’s total flux was captured across the full wavelength range. Finally, the individual exposures were astrometrically aligned to correct for coordinate offsets and combined by averaging, resulting in the final calibrated datacube.

3.3. Datasets

We present the first 30 systems observed during the first year of the Cosmic Duets Large Program (ESO 114.27BY; PI: Scialpi). Table 1 summarizes the main information of these targets, including exposure times (T_exp), number of exposures (N_exp), projected separations among the components (Sep), and preliminary classification. The following subsections describe the extraction of images and spectra from the calibrated MUSE-NFM datacubes.

3.3.1. Image extraction

For each target, we extracted a white-light image by summing up all spectral channels. Figure 2 shows cut-outs of the 16 double systems (lensed or dual AGN) with logarithmic scaling to enhance faint sources. The separation between the components is indicated in each panel. Figure 3 shows the three quadruply imaged systems, with the radius of the smallest circle encompassing AGN images (Einstein radius, R_E) indicated. All 30 observed systems are resolved into multiple point sources with angular separations ranging from 0.2″ to 0.8″, confirming the success of the GMP selection method. Figure 1 shows the distribution of dual AGN and lensed systems from this sample and from Scialpi et al. (2024), all observed with MUSE–NFM, as a function of projected separation and redshift.

Fig. 1. Distribution of dual (green) and lensed (blue) AGN as a function of projected separation (top panel) and redshift (bottom panel). The sample, observed with MUSE–NFM, includes sources presented in this work and by Scialpi et al. (2024). For quads, the Einstein radius was used to characterize the separation.

3.3.2. Spectral extraction

All spectra were extracted using a fixed 5 pixel-radius (0.125″) circular aperture, chosen to match the typical FWHM of the PSF across all wavelengths. For J2344–4259 and J0530–3730, a smaller aperture of 4 pix was used because the components are closer than 0.25″. This choice optimizes the trade-off between flux recovery and contamination from the companion sources. The extracted spectra are independent, with flux uncertainties of ∼10%.

To assess whether aperture extraction is affected by cross-contamination (Husemann et al. 2018; Pfeifle et al. 2023), we performed an independent test using profile-fitting extraction with Gaussian and Moffat models. A white-light image was created by summing all spectral channels, and fit with two Gaussians to obtain initial estimates for the centroids and widths (σ_x, σ_y). Because the sources are very close, we fixed their relative separation and imposed common widths, as PSF variations at this scale are negligible compared to flux uncertainties.

Using these initial parameters, we then fit two Gaussians independently in each spectral channel. The resulting fit parameters, the widths (σ_x, fit, σ_y, fit) and amplitude (A_fit), were smoothed as a function of wavelength by fitting a polynomial of degree 3, in order to remove outliers caused by noisy or bad channels. The flux of each component was then computed as

$$\begin{array}{c}\hfill F_A=A_{\mathrm{fit}}2\pi {\sigma }_{x,\mathrm{fit}}{\sigma }_{y,\mathrm{fit}},\end{array}$$(2)

where A_fit is the fit amplitude of the Gaussian for component A, and σ_x, fit and σ_y, fit are the fit Gaussian widths along the x and y directions, respectively. This procedure ensures that the extracted flux accounts for the PSF shape and minimizes contamination from the nearby source.

This test allowed us to evaluate whether contamination or wavelength–dependent PSF variations systematically alter the flux recovered by fixed apertures. At blue wavelengths, however, the PSF broadens significantly due to the reduced efficiency of AO correction, causing the fit widths to increase and the profile–fitting technique to overestimate the true flux. Tests using Moffat profiles or mixed models (two Gaussians or two Moffats per component) did not improve the situation: despite their greater flexibility, they proved highly sensitive to noise and channel–to–channel fluctuations, leading to spurious variations in the extracted flux at low S/N.

By comparison, circular-aperture extraction yields stable measurements at all wavelengths, with flux variations below 5% and negligible cross-contamination (< 1%) thanks to the small separation and the excellent MUSE–NFM spatial resolution. Because AO assisted optical observations have modest Strehl ratios, the extracted spectra contain only a fraction of the total flux; however, this fraction is identical for both components in each pair. We therefore normalized the spectra to their Gaia magnitudes (Gaia Collaboration 2021). When Gaia resolves both components, each spectrum was normalized to its own magnitude; when only a combined magnitude is available, we first measured the flux ratio between components and normalized the sum of their spectra to the system magnitude.

These normalization factors reflect not only the AO Strehl ratio, but also the absolute flux–calibration accuracy of MUSE, which contributes significantly to the overall uncertainty. Assuming similar Strehl ratios and PSF for all components, we applied these normalization factors to extract the final spectra from each datacube. The resulting normalized spectra are shown in Figs. A.1 (dual AGN) and A.2 (double lensed AGN). Figure 2 shows the corresponding circular apertures used for extraction. In all the figures, the brightest component is shown in purple and the faintest in blue. For the three quadruple-lensed AGN (Figs. 3 and A.3), the third and fourth images are shown in orange and green, respectively, while the lens galaxy appears in red when detected.

Fig. 2. White-light images of the double AGN systems (lensed or dual) in units of 10−15 erg/s/cm−2. The target name and the separation between the two components are indicated in each panel. The spectra are extracted from the circular apertures shown on the maps, with purple and blue apertures corresponding to the A (brightest) and B (faintest) components, respectively. Dual systems are marked with the flag D, lensed systems with L. For systems identified as gravitational lenses, if an additional lensed image is detected, it is shown inside the red aperture and labeled as component G.

Fig. 3. White-light images of the three quadruply lensed AGN systems. The target name and the value of the Einstein radius (Re) is indicated in each panel. The spectra are extracted from the colored circular apertures. The colors of the apertures follow the order of brightness: A (purple, brightest), B (blue), C (orange), and D (green, faintest). For J0957, we also detect the emission-line lensing source (E) in red.

4. Data analysis

To analyze the data, we followed the methodology applied in a previous MUSE-NFM pilot program adopting the same observational strategy, and target selection described in Scialpi et al. (2024). Details on the main steps of the analysis are provided in the following.

4.1. Spectral modeling

For systems classified as AGN+AGN, the extracted spectra of each component were modeled independently using AGN templates from Vanden Berk et al. (2001) and Temple et al. (2021), including interstellar dust extinction following the empirical quasar attenuation law of Temple et al. (2021). This fitting aims to derive the reddening (E(B − V)) and redshift of each AGN individually. Because the spectra were extracted from circular apertures and the components are well separated, we treat each spectrum as originating from a single source; potential differences in emission-line profiles between the two AGN are thus captured within the respective template fits. Galactic extinction is neglected, as its effect is minimal at the high Galactic latitudes of our targets (|b|> 20°).

For systems composed of an AGN and a foreground star, we modeled the AGN emission using the same template-fitting procedure described above, while the stellar component was fit with templates from the Yan et al. (2019) library, including G-, K-, and M-type stars to determine its spectral type.

4.2. Narrow absorption lines

In addition to continuum and emission-line diagnostics, we identified and analyzed narrow absorption lines (NALs) in our MUSE spectra, which in some cases provide complementary information to distinguish lensed from dual AGN systems. Intrinsic absorbers are close in redshift to the quasar, typically with velocity shifts of < 3000 km s⁻¹ relative to the systemic redshift. These absorbers, often showing blueshifted velocities, consist of gas associated with the quasar or its host galaxy, and intrinsic NALs trace the quasar environment (e.g., Hamann et al. 2011; Misawa et al. 2007; Lu & Lin 2019). In lensed systems, if one image exhibits intrinsic NALs, the other image should generally show the same features. Differences in equivalent width (EW) can arise due to microlensing (see Sect. 4.3) or the geometry of the lens, where each light path intersects different regions of the lens galaxy (e.g., Chartas et al. 2009; Green 2006). Conversely, if intrinsic absorbers differ significantly between components in redshift, depth, or profile, this indicates that the two sight lines probe physically distinct AGN, favoring a dual AGN scenario.

Intervening absorbers consist of gas clouds or galaxies along the line of sight, unrelated to the quasar itself. These features appear at redshifts different from that of the quasar, but can be present at the same redshift in both components of a system if both sight lines intersect the same intervening structure (Martin et al. 2010). While identical intervening NALs in both spectra are consistent with lensing, this is not by itself conclusive, and must be interpreted alongside other diagnostics. Even simple measurements of absorber redshifts provide valuable information on the intervening structures along closely spaced sight lines (separations of a few kiloparsecs at the absorber redshift), offering a first step toward probing the small-scale properties of the circumgalactic and intergalactic medium (e.g., Dutta et al. 2024). To search for NALs in our sample, each spectrum is normalized by a locally fit continuum and inspected for narrow features. All absorber redshifts have been measured, providing a comprehensive catalog of both intrinsic and intervening NALs in our sample. For intervening absorbers, a line is considered common (from the same absorber) if it is present in both components with a velocity difference Δv < 50 km s⁻¹ and consistent depth ΔEW/EW < 20%. For intrinsic absorbers, differences between components can indicate dual AGN or quasar-driven outflows.

4.3. Classification of dual and lensed systems

To discriminate between lensed QSOs and dual AGN, we compared the two components in terms of redshift, EWs, emission-line flux ratios, line profiles, continuum shape and NALs. Before classifying a system as a dual AGN, it is important to verify that observed differences in continuum or broad-line EWs cannot be explained by gravitational microlensing. Microlensing occurs when compact objects (e.g., stars) in the lens galaxy preferentially magnify the most compact emission regions, such as the accretion disk and the high-ionization portion of the broad-line region (e.g., Wambsganss 2006; Mosquera et al. 2013; Jiménez-Vicente et al. 2015). Macro- and microlensing magnifications (M and μ) were estimated via the Macro–Micro Decomposition (MmD) method (Sluse et al. 2007). Deviations caused by microlensing are recognized by their wavelength-dependent signatures and by the absence of corresponding changes in narrow lines, ensuring they are not misinterpreted as evidence of dual AGN.

Below we summarize the criteria adopted for classification:

• Redshift (Δz). Since lensed images originate from the same background AGN, they must share the same intrinsic redshift. Gravitational lensing cannot produce differential line-of-sight velocity shifts larger than the small, apparent differences due to microlensing on the broad line shape. Therefore any measurable offset (Δz ≥ 0.001, i.e., Δv ≳ 300 km s⁻¹) cannot be produced by lensing and indicates two physically distinct AGN.

• Equivalent widths. For each rest-frame spectrum extracted from circular apertures, line fluxes were measured via trapezoidal integration after estimating the local continuum from the edges of a velocity window. Uncertainties were propagated from individual flux errors and continuum estimates. The EWs were then compared between components. Significant differences not compatible with microlensing, variability, or differential extinction, support a dual AGN interpretation.

• Emission-line ratios. Line flux ratios, especially for narrow emission lines, are a robust diagnostic. Significant deviations between components indicate physically distinct AGN, since narrow-line regions are spatially extended and largely unaffected by microlensing.

• Continuum shape and line profiles. Differences in continuum or line profiles were quantified via a cross-correlation procedure, progressively modifying one spectrum to match the other and optimizing parameters via χ² minimization. Initial fits considered only a redshift offset and flux scaling; further refinements allowed independent variations in continuum and emission lines, including wavelength-dependent extinction modeled with a second-degree polynomial. Minor deviations consistent with these effects do not indicate a dual AGN.

• Narrow absorption lines. Narrow absorption lines provide complementary information for classification (see Sect. 4.2). We identify both intrinsic and intervening absorbers. Each spectrum is normalized by a locally fit continuum, and absorber redshifts are measured. For intervening absorbers, a line is considered common if present in both components with Δv < 50 km s⁻¹. For intrinsic absorbers, we consider differences between the two components to be significant if the absorber centroids differ by Δv > 300 km s⁻¹. In these cases, gravitational lensing (including differential potential and time-delay effects) cannot produce systematic centroid shifts and they are naturally explained by two distinct sight lines intersecting different absorbing regions (outflows/host-galaxy kinematics).

Overall, a system is classified as a lensed AGN if all diagnostics are consistent with the same source (allowing for microlensing and extinction effects), and as a dual AGN if one or more diagnostics show significant differences.

5. Results

Table 1 summarizes the classifications obtained for the observed systems. The cross-correlation analysis described in Sect. 4 was used as a quantitative diagnostic to assess spectral similarity between components, particularly for identifying lensed AGN candidates. However, the final classification of each system as dual AGN, lensed AGN, or AGN+star projection is based on a combination of spectral, kinematic, and morphological criteria, including emission-line profiles and ratios, velocity offsets, flux ratios, and the presence of lensing features in the datacube. Applying this multicriteria approach, we identified six systems as dual AGN, ten as doubly lensed AGN candidates, three as quadruply lensed systems, and 11 AGN+star. In the following, we describe the observed systems in detail, highlighting their spectral properties and the procedures adopted for their final classification.

5.1. AGN+star systems

The spectral decomposition used during the target selection stage (see Sect. 2.2) is designed to remove stellar contaminants prior to IFU follow-up and succeeds in excluding the majority of such systems; see Sect. 4. Nevertheless, in our sample of 30 observed systems, 11 are revealed by the MUSE data to be AGN+star projections that were not rejected during the initial selection. These cases correspond to configurations in which the stellar component is significantly fainter than the AGN, so that its spectral features are too faint to be detected with sufficient confidence in the integrated spectrum used for preselection.

In addition, several spatially unresolved observations were obtained under atmospheric conditions that reduced the effective S/N (seeing of 0.8 − 1.2 arcsec), further limiting the detectability of stellar absorption features. Once the spatial and spectral information of the IFU is taken into account, the nature of the secondary point source becomes evident, and the stellar component can be isolated and fit, allowing a robust AGN+star classification. The spectra of these systems are reported in Fig. C.1 for completeness.

Considering all GMP systems observed with MUSE (this work and Scialpi et al. 2024), we find that the occurrence of AGN+star projections strongly depends on the projected separation between the two point sources. Figure 4 shows the fraction of AGN+star systems, N_{AGN + star}/N_tot, in bins of separation.

Fig. 4. Fraction of AGN+star projections as a function of projected separation. The dotted red line shows the expected contamination probability function (Eq. 3 with D = 0.65″), representing the separation distribution of random stellar contaminants. The shaded area shows the separation range excluded by the GMP selection.

The probability of stellar contamination can be described by

$$\begin{array}{c}\hfill P(r)=1-exp(-\frac{\pi r^2}{4D^2})=1-exp(-\pi r^2\rho ),\end{array}$$(3)

where D is the mean distance to the nearest contaminant star and ρ = 1/4D² is the surface density of contaminants. For r ≪ D, P(r) behaves approximately as a parabola (Banerjee & Abel 2021), while it asymptotically approaches one for r ≫ D. For illustration, the curve with D = 0.65″ is overplotted in Fig. 4 in red. The contamination fraction rises significantly for separations ≳0.6″, where AGN+star projections become more common than double AGN candidates. At separations > 0.8″, contamination exceeds 80%, and the curve’s asymptotic approach to 1 indicating that dual AGN searches at these scales are expected to have low efficiency. These pairs are therefore excluded from our GMP sample (dashed area). These levels agree with estimates by Mannucci et al. (2022, P(0.5)∼0.3) and Ulivi et al. (2025), who found a lower critical separation (0.4″) due to deeper sensitivity, resulting in a smaller value of D and therefore higher contamination.

5.2. Dual AGN systems

Six of the AGN+AGN systems in our sample have been classified as dual AGN following our spectroscopic and spatially resolved analysis. Only one of these sources, J0032–1053, was also selected by an independent method and previously identified as a dual AGN candidate by Hwang et al. (2020), who flagged the object due to significant non-zero proper motions in Gaia DR2. In the following, we describe the analysis and properties of each individual system.

5.2.1. J0032–1053

J0032–1053 is confirmed as a dual AGN system at redshift z = 2.439 ± 0.001, with the two nuclear components separated by a projected distance of 5.5 ± 0.1 kpc and a relative velocity of ∼200 ± 130 km s⁻¹. Spectroscopic analysis reveals notable differences between the two nuclei, labeled AGN A and AGN B (purple and blue, respectively, in Fig. A.1). AGN A is more luminous, with its observed continuum flux requiring a downward scaling factor of ∼3.2 for direct spectral comparison. The line ratios also differ: the lower-ionization species, such as C II] λ2326, are proportionally stronger in the AGN B spectrum relative to the high-ionization C IVλ1549 emission.

Two intervening absorption systems are detected in both components (see Fig. 5). The first system is located at z = 1.0794 ± 0.0004, where the Mg IIλλ2796.35, 2803.53 doublet is observed at ∼6943 Å, along with Fe IIλλ2586.65, 2600.17 transitions around 6448 Å in the observed frame. A second intervening absorber is detected at z = 1.4791 ± 0.0006, showing multiple Fe II transitions (λ2344, λ2374, λ2382, λ2586, λ2600) spanning ∼4800 − 5400 Å, as illustrated in Fig. 5.

Fig. 5. Absorption systems of the dual AGN J0032–1053. The upper panel shows component A, while the bottom panel shows component B. The two intervening absorption systems are indicated: the system at z = 1.079 is shown in dark red and the system at z = 1.479 in teal. Intrinsic C IV absorption at the AGN redshifts is shown in blue.

In both AGN A and AGN B, C IV absorption is detected, but with markedly different velocity shifts relative to their systemic redshifts. For AGN B, adopting a systemic redshift of z_sys, B = 2.442, the C IVλλ1548, 1550 doublet is detected at z_out, B = 2.4388, corresponding to a blueshift of v_B → B ≃ −170 km s⁻¹, well within the typical range of AGN outflows (e.g., a few hundred km s⁻¹; Hamann et al. 2011). The line ratio is consistent with an intrinsic narrow absorption line (NAL), indicating absorbing gas associated with AGN B, but without a fast outflow. For AGN A, with systemic redshift z_sys, A = 2.43967, a weaker doublet-like feature is detected at z_out, A = 2.4175. Interpreted as intrinsic C IV absorption, this implies a much larger blueshift of v_A → A ≃ −1935 km s⁻¹, consistent with a fast quasar outflow. Compared to the systemic redshift of AGN B, the offset becomes v_A → B ≃ −2140 km s⁻¹. The shallow profile and non-unity doublet ratio suggest that the sight line to AGN A intercepts a lower-density, possibly more diffuse region of the same outflow traced in AGN B, but at lower optical depth.

Alternatively, the absorption feature in AGN A may originate from a distinct intervening system unrelated to the intrinsic outflow seen in AGN B. In this case, the absence of C IV absorption at the systemic redshift of AGN A supports the interpretation that only AGN B exhibits an intrinsic outflow, while the two sight lines do not intersect the same absorbing structure, despite their projected 5.5 kpc separation.

5.2.2. J0144–5745 and J2100+0612

These two systems are confirmed dual AGN at redshifts z = 0.791 and z = 0.779, with projected separations of 3.4 ± 0.1 kpc and 3.9 ± 0.1 kpc, respectively. In the MUSE spectra at these redshifts, we detect the Mg IIλλ2796, 2803 doublet at the bluest wavelength, Hβ and [O III] λλ4959, 5007 at the reddest wavelength, and additional lines at intermediate wavelengths. In both systems, the brightest AGN (purple in Fig. A.1) exhibits a broader Hβ line, and in J2100+0612 a broader Mg II is also observed. Additionally, the line ratios between the narrow lines differ between the two nuclei, with measured ratios of [O III] λ5007/Hβ = 4.2 ± 0.2 for AGN A and 7.8 ± 0.2 for AGN B. No narrow absorption line systems are detected in either component of these two dual AGN. These differences suggest varying ionization conditions or metallicities between the nuclei

5.2.3. J0637–4137

The system J0637–4137 at z = 2.743 ± 0.001 is a strong candidate for a dual AGN, with the two components separated by a projected distance of 3.3 ± 0.1 kpc. The shape and EW of the C III] line in the fainter component (light blue in Fig. A.1) differ from those of the brighter nucleus, providing evidence that the system hosts two distinct AGN.

Multiple absorption systems are detected in both spectra (Fig. A.1). All of these are intervening and thus not physically associated with the quasars themselves. We identify three such systems at z_abs1 = 1.2523, z_abs2 = 2.1352, and z_abs3 = 2.3924, each showing consistent metal absorption features (e.g., Mg II, Fe II, C IV, Si II) in both components. The close match in absorption redshifts along the two sight-lines provides independent confirmation that the quasar pair is physically associated, complementing the evidence from the nearly identical emission-line redshifts and reducing the likelihood of a chance line-of-sight alignment.

5.2.4. J1837–6711

The system J1837–6711 is a dual AGN. First, the narrow emission lines differ between the two components, with [O II], [Ne III], Hϵ, Hδ, and Hγ showing distinct flux ratios. Since the NLR is spatially extended, microlensing cannot alter its flux, and these differences therefore require two separate narrow-line regions. Second, the Fe II emission complexes around Mg II also differ. The blue side of Mg II shows stronger absorption near ∼2700 Å in the purple component, and higher bumps at ∼2900 Å and ∼3200 Å. Fe II originates in the outer BLR, which is only weakly and smoothly affected by microlensing. Thus, the observed Fe II variations cannot be produced by microlensing of a single BLR, and instead imply intrinsically distinct BLR conditions. Taken together, these differences in the spectral properties cannot be explained by microlensing of a single quasar. They instead strongly indicate the presence of two active nuclei in the same system, consistent with an ongoing galaxy merger at z = 1.120.

5.2.5. J2204–6530

J2204–6530 is confirmed as a dual AGN system with a projected separation of 3.6 kpc. The two components show clear differences in velocity, line centroids, and line profiles, consistent with the presence of two distinct active nuclei.

Using the spectral decomposition and fitting procedure described in Sect. 4, we measured redshifts of z_A = 1.880 and z_B = 1.892. These values are based on broad emission lines (C III] and Mg II; the latter is partially affected by absorption), as no narrow lines suitable for precise systemic measurements are available in the spectral range (see Fig. A.1). Accounting for these limitations, we estimate a relative velocity of ∼1236 ± 250 km s⁻¹ between the two components.

This velocity shift cannot be explained within a gravitational lensing scenario. In lensed systems, the spectra of multiple images are expected to be nearly identical, as the typical time delays (of order days; Treu 2010) are too short to produce such large differences in emission-line centroids and profiles through intrinsic variability.

Although the measured offset lies at the high end of the velocity distribution observed for close dual AGN, it remains within the range reported in the literature (e.g., Pfeifle et al. 2025), particularly when accounting for the large uncertainties associated.

As shown in Fig. 6, both AGN display intrinsic NALs at different velocities. In both spectra, the Mg IIλλ2796, 2803 doublet is detected in absorption, although at different observed wavelengths. For AGN A, the Mg II NALs are redshifted by 700 ± 45 km s⁻¹ relative to its systemic redshift, while in AGN B the absorption appears at −1500 ± 300 km s⁻¹ relative to AGN B, and at −250 ± 50 km s⁻¹ when referred to the systemic redshift of AGN A.

Fig. 6. Enlargement of the Mg II spectral region of J2204–6530 (AGN A in the top panel and AGN B in the bottom one). Both AGN show the Mg II doublet absorption lines at rest-frame wavelengths of 2796 Å and 2803 Å. In each panel, the redshift of each QSO and the NAL velocity relative to its systemic redshift are reported, assuming that the two outflows are independent.

This configuration can be interpreted in two ways: each AGN hosts its own independent outflow with distinct velocities; a single outflow (likely associated with one AGN) intersects the line of sight of both nuclei, producing absorption at different projected velocities depending on geometry and relative orientation. While these velocities are consistent with outflowing gas, the current data do not allow a precise characterization of the outflow properties. Distinguishing between alternative scenarios requires spatially resolved kinematic analysis of the absorbing gas, which is beyond the scope of this work and will be addressed in future follow-up studies. We note that inflowing gas is unlikely to explain the observed absorption, as the measured velocity offsets are mostly blueshifted relative to the systemic redshifts, and redshifted absorption expected from inflows is not observed. The line profiles and partial coverage signatures also favor an interpretation as AGN-driven outflows rather than accreting material.

5.3. Double lensed AGN

In this section we present the systems classified as doubly lensed AGN from our data and analysis. Out of the sixteen double AGN systems analyzed, ten have consistent redshifts between their components and show minor spectral differences that are fully consistent with microlensing effects (see Fig. A.2). In these cases, variations appear primarily in the continuum flux and, in some instances, in the equivalent widths of broad emission lines, while narrow lines remain unaffected. No intrinsic differences between the images are required to explain the observations. These signatures are expected for gravitational microlensing, where compact objects in the lensing galaxy (e.g., stars) preferentially magnify the most compact emission regions, such as the accretion disk and the high-ionization portion of the broad line region (Wambsganss 2006). Consequently, the blue continuum and high-ionization lines can be more strongly magnified than the red continuum or low-ionization lines (e.g., Mosquera et al. 2013; Jiménez-Vicente et al. 2015). Narrow emission lines, emitted on scales of hundreds to thousands of parsecs, are essentially immune to microlensing effects. Using the macro–micro decomposition (MmD) method (Sluse et al. 2007) on all available emission lines, we obtained detailed estimates of M and μ for nearly all sources in the sample. Systems with values consistent with literature expectations (μ ∼ 1.1 − 1.5 for microlensed continuum regions, while narrow-line regions remain unaffected) are classified as lensed AGN. The observed differences in these ten systems are therefore fully explained by microlensing alone, supporting their classification as lensed AGN candidates.

The systems J0259–0901, J0551–4629, J1447–0501 and J2257–6555 show nearly identical spectra from both images, with no significant differences detected. This strongly supports their classification as lensed AGN, with each pair of components corresponding to multiple images of the same source. In the MUSE datacube of J0551–4629, the lensing galaxy is directly detected and is highlighted in red as component C in Fig. 2, providing a definitive confirmation of the lensing nature of the system.

J0802+1944, J2344–4259, and J2353–0655 show some differences in the continuum shape, with marginally significant effects on the broad or narrow emission lines within the measurement flux errors. In these cases, microlensing can primarily account for the observed variations in the continuum, driving the lensed classification of the systems.

While J0034–1623, J0317–5604, and J0404–2836 exhibit differences in their Mg II emission line profiles, these variations are quantified by the differential microlensing factor μ = A/M, where A is the continuum flux ratio and M the Mg II line flux ratio between the images. The resulting μ values agree with differential amplification expected from the literature (Sluse et al. 2007, 2012). For J0034–1623, A ≈ 5.60 and M ≈ 6.21, yielding μ ≈ 0.90, indicating that Mg II is de-magnified relative to the compact continuum source and explaining the observed differences in line shape and EW. For J0317–5604 (A ≈ 1.57, M ≈ 1.46, μ ≈ 1.07) and J0404–2836 (A ≈ 1.82, M ≈ 1.73, μ ≈ 1.05), the Mg II BLR is only marginally affected. While physically the same in both images, microlensing induces small variations in the EWs due to the different amplification of the continuum versus the more extended BLR. Narrow lines in J0404–2836 remain unaffected (μ ≈ 1), confirming that the observed differences in broad lines are due to microlensing. Although dual AGN could, in principle, produce differences in Mg II, this analysis strongly supports microlensing as the main cause.

Several of the systems in our sample show narrow absorption features at redshifts lower than the quasar systemic redshift. The fact that these absorbers lie at lower redshifts than the background AGN indicates that they are produced by foreground structures along the line of sight rather than by intrinsic outflows from the quasar. For instance, J0034–1623 shows a strong absorber at z = 1.6339, while J0317–5604 hosts three intervening systems at z = 1.4204, 1.6146, and 1.6372. Similar features are found in J0404–2836 (z = 1.1231, 1.2510, 1.3038), and in J0551–4629 (z = 0.7903, 1.0151), among others. The presence of the same absorption systems in both images confirms that the light paths of the two images cross the same foreground structures, supporting the gravitational lens interpretation.

In contrast, two systems show absorption features that are not intervening but intrinsic to the quasar. In J1447–0501, the Mg II doublet (λ_rest = 2796, 2803 Å) is observed in absorption, with a velocity of ∼200 ± 50 km s⁻¹ in both images, indicating a moderate quasar-driven outflow. In J2344–4259, several high-ionization lines show stronger blueshifted absorption, including the N V doublet (λ_rest = 1238.82, 1242.80 Å) and the C IV doublet (λ_rest = 1548.20, 1550.77 Å), with velocities of −2748 ± 40 km s⁻¹. These features are also detected in both components, confirming their intrinsic origin and therefore the lensed nature of the system, rather than indicating distinct physical AGN.

5.4. Quadruple lensed AGN

From our GMP selection procedure (Sect. 2), we identified three systems that each exhibit four AGN images (quads; see Fig. 3). For each system, we extracted both the images and the spectra of the four AGN. In one case, we were also able to isolate the spectrum of the lensing galaxy.

We decided to observe even those systems previously known as quads from imaging in order to obtain MUSE spatially resolved spectroscopy. When the lens is sufficiently bright, this allowed us to directly detect and resolve both the lensing galaxy and the multiple AGN components. The high angular resolution provided by AO assisted MUSE observations (∼0.1″) is critical for accurate modeling of such small-scale lens systems. Moreover, this approach ensures that the measurement of lensing prevalence among spectroscopic GMP candidates is independent of prior assumptions.

5.4.1. J0530–3730

This system was identified as a quads candidate by Delchambre et al. (2019) based on a search over Gaia DR2 (Gaia Collaboration 2018). It was also selected using the method described by Ostrovski et al. (2017) and Lemon et al. (2017) as a triple detection around a photometric quasar candidate, and through the GMP technique (Mannucci et al. 2022), because AGN A and C are separated by 0.16″ and AGN A and B by 0.67″ (see Fig. 3). It was confirmed as a quasar at z = 2.838 from NTT spectra (Anguita et al. 2018).

This system was modeled using the automated time-delay cosmography pipeline by Schmidt et al. (2023), based on high-resolution HST multi-band imaging and informative priors. The method provides accurate strong-lens mass modeling, deriving key lens parameters, magnifications, and time-delay predictions.

Using our data, we measured a circle enclosing all four AGN images with a radius of 0.531 ± 0.009″, consistent with previous estimates (Schmidt et al. 2023). Thanks to the high spatial resolution of our observations, we are able, for the first time, to simultaneously spatially and spectrally resolve all four AGN, detecting in all of them the same absorption systems at z = 1.0054 ± 0.0002 and z = 1.2641 ± 0.0002.

5.4.2. J0957–2242

J0957–2242 is a newly discovered quadruple lensed system observed in this Large Program, as it is exclusively GMP-selected. MUSE data reveal the four images of the AGN at z = 2.070 ± 0.005, as well as the foreground lens galaxy, which is detected in emission (component E in Fig. 3). We extract the spectra of the AGN from circular apertures of 5 px (0.125″) radius and that of the lens galaxy from a smaller aperture of 3 px (0.075″) radius (see the bottom panels of Fig. A.3). The smaller aperture for the lens galaxy minimizes contamination from the much brighter AGN emission, while still capturing the central light of the deflector. We estimate an Einstein radius of 0.32″ (∼1.7 kpc), making this one of the most compact quadruply imaged AGN currently known, with an emission-line galaxy acting as the lens (D’Amato, in prep.).

From the MUSE spectra, we estimate the redshift of the lens source to be z = 0.492 ± 0.001, we detect Hβ, the [O III] λλ4959, 5007 doublet, and the [O II] λλ3726, 3729 doublet (Fig. A.3), and we measured the individual line fluxes. These line fluxes can be used to discriminate between photoionization by star formation and AGN activity. However, as shown by Baldwin–Phillips–Terlevich (BPT) emission-line diagnostic diagrams (Baldwin et al. 1981) in Fig. 7 adapted from Lamareille (2010), the available lines are insufficient to confirm or exclude the presence of an AGN. The log₁₀[O III]/Hβ = 0.36 ± 0.2, consistent with both star-forming and AGN-dominated ionization, while the log₁₀[O II]/Hβ = 0.65 ± 0.2, placing the galaxy near the separation line. Accurate classification would require the [N II] λ6584/Hα ratio, which is not covered by the MUSE spectra, following the diagnostic scheme of Lamareille (2010).

Fig. 7. BPT diagram adapted from Lamareille (2010) showing the positions of the local SDSS galaxies. In the left panel, the yellow strip shows the range of [N II]λ6584/Hα values corresponding to the observed [O III]/Hβ ratio (the dashed red area indicates the measurement uncertainty). In the right panel, the yellow star marks the position of the lens galaxy in the [O III]/Hβ versus [O II]/Hβ diagram, which lies very close to the separation line between star-forming galaxies (blue) and AGN (green and cyan).

The ratio of the [O II] lines, O IIλ3726/λ3729 = 0.85 ± 0.10, implies an electron density of n_e ≃ 263 ± 50 cm⁻³, assuming an electron temperature of T_e = 10⁴ K (Osterbrock & Ferland 2006). These relations are valid for ionized gas typical of star-forming galaxies, where collisional excitation dominates. Considering a plausible range of electron temperatures for low-z star-forming galaxies, T_e = 7000 − 18 000 K, the inferred density varies only moderately, from n_e ≃ 223 cm⁻³ to n_e ≃ 313 cm⁻³, remaining within the quoted uncertainties.

We computed the oxygen abundance using the R23 metallicity indicator, defined as R23 = ([O II] λλ3726,3729 + [O III] λ4959, λ5007)/Hβ (Curti et al. 2020). From the measured fluxes we obtain log₁₀(R23) = 0.91 ± 0.20. Using the empirical calibration of Maiolino et al. (2008), this corresponds to two possible oxygen abundances: 12 + log(O/H) = 7.85 ± 0.20 (low–metallicity branch) and 12 + log(O/H) = 8.44 ± 0.20 (high–metallicity branch). To break the degeneracy we computed the R2 and R3 ratios, obtaining log R2 = 0.65 ± 0.15 and log R3 = 0.55 ± 0.15. Both values place the source on the high-metallicity branch according to the trends of Maiolino et al. (2008). The agreement between R2 (dust-sensitive) and R3 (dust-insensitive) also indicates minimal nebular extinction. We therefore adopted the high-metallicity solution 12 + log(O/H) = 8.44 ± 0.20 consistent with moderately metal-rich ionized gas.

Using the observed Hβ luminosity and the Kennicutt (1998) calibration we estimate the star formation rate (SFR) SFR_Hβ = 0.016 ± 0.002 M_⊙ yr⁻¹ (uncorrected for dust). Because R2 and R3 are in agreement, suggesting minimal dust extinction, this value should be regarded as a robust lower limit and is likely close to the true SFR.

Overall, the combination of these emission-line ratios and diagnostics indicates highly ionized gas with moderate electron density, providing meaningful constraints on the excitation mechanism and chemical enrichment even in the absence of Hα and [N II] coverage.

The ionized gas kinematics traced by the [O III]5007 transition reveal a tentative rotation pattern, as shown by moment maps in the upper panel in Fig. 8. To test whether a pure rotating disk is suited to reproduce the observed kinematic features we adopted the MOKA^3D kinematic model (Marconcini et al. 2023). In particular, we assumed a circular thin disk with inner and outer radius of 0−0.15″ (∼0−0.6 kpc) and a disk height comparable with the instrumental PSF (i.e., 0.1″). In doing so, we assumed that the gas kinematics can be described as a purely rotating disk with constant circular velocity, with no radial motions. The free parameters of the fit are the inclination, the circular velocity and the disk PA. Our best fit model returns an inclination of 55° ±3°, an intrinsic circular velocity of 380 ± 11 km s⁻¹ and a PA of 18° ±4°.

Fig. 8. MOKA3D rotating disk best-fit model of J0957−2242 traced by the [O III]5007 transition. Panels show the observed (top), best-fit (center), and residual (bottom) moment maps. The residual maps were obtained by subtracting the model from the observed moment maps. The star marks the kinematic center, inferred as the peak flux emission. Pixels are masked at SN ≤ 4.

As shown in Fig. 8 the inferred best-fit model is able to reproduce the observed ionized gas features with extremely high-accuracy and residuals ≤10 km s⁻¹ both in the line-of-sight velocity and velocity dispersion maps. Then we derived the dynamical mass profile from the deprojected circular velocity while assuming circular motion in a spherically symmetric potential:

$$\begin{array}{c}\hfill M=\frac{V_{\mathrm{circ}}^{2}R}{G},\end{array}$$(4)

where R is the galactocentric radius in kpc, and G is is the gravitational constant. We computed the uncertainties on the dynamical mass propagating the uncertainties of the circular velocity. This procedure provides a direct estimate of the enclosed dynamical mass, under the assumption that the gas trace circular motion in the plane of the disk and that non-circular motions are negligible. As a result, we estimate a dynamical mass within the disk maximum extension (0.15″) of M = 3.1 ± 0.2 × 10¹⁰ M_⊙.

5.4.3. J1116–0657

The quasar J1116–0657 was identified as a quadruple gravitational lens by Blackburne et al. (2008), with a source redshift of z_S = 1.235. They detect a faint lensing galaxy estimated at z_L ∼ 0.7. A simple lens model reproduces the observed image positions but fails to match the flux ratios, possibly due to substructure or microlensing. The estimated time delays between the images are on the order of one day.

Our measurements are in agreement with Blackburne et al. (2008), with an Einstein radius of 0.341 ± 0.005″ (see Fig. 3). Importantly, we were able to extract spectra for all four AGN images individually, providing the first resolved spectral information for this system and enabling the identification of spectral features in each component.

6. Distribution of dual and lensed systems

Building on the classifications derived for each target, we next analyzed the population properties of the GMP-selected systems observed with MUSE–NFM. For this purpose we consider the full sample with available MUSE follow-up, which includes the 19 systems presented in this work together with the seven additional sources from Scialpi et al. (2024). This consolidated sample provides a uniform basis to investigate the redshift, magnitude, and separation distributions of dual and lensed AGN, and understand the selection effects introduced by both the GMP method and the MUSE observations.

The sample is limited to projected separations 0.15″ < sep < 0.8″. The lower limit corresponds to the sensitivity of the GMP method, while the upper limit reflects our focus on systems in the regime of a few kiloparsecs, where simulations predict that the SMBHs are likely embedded within a common merger remnant or in the late stages of galaxy coalescence (Volonteri et al. 2022). At the median redshift of the sample (z ∼ 1.5), these values correspond to projected physical separations of ∼1.7−6.8 kpc.

Only systems with z > 0.5 are included, ensuring that the GMP technique identifies spatially distinct components and not single extended sources.

Figure 9 shows the fraction of lensed systems (blue) and dual AGN (green) as a function of J-band magnitude (top panel) and redshift (bottom panel) for the MUSE sample (this work and Scialpi et al. 2024). In both cases, the fractions are computed per bin so that the sum of lensed and dual systems equals one, providing a normalized view of the relative contributions removing the selection effects. Lensed systems clearly dominate the bright end of the sample (J ≲ 16.5), while dual AGN become increasingly prevalent at fainter magnitudes, highlighting the strong influence of magnification bias. The lensing fraction peaks around z ∼ 2.1, reflecting the combination of survey selection effects (also relevant for dual AGN), the increasing lensing efficiency due to the ratio of angular diameter distances of lenses and sources, and the quasar luminosity function.

Fig. 9. Top panel: Fraction of lensed systems (blue) and dual AGN (green) as a function of J-band magnitude for the MUSE sample (this work and Scialpi et al. 2024). Bottom panel: Same but as a function of redshift. Data are binned in intervals of 0.5 mag (top) and 0.5 (bottom). The number above each column indicates the total number of objects in that bin. The horizontal dashed line marks the 50% level for reference.

These distributions can be interpreted in the broader context of dual AGN searches. Relative to previous GMP studies (Mannucci et al. 2022, 2023; Ciurlo et al. 2023; Scialpi et al. 2024), the present MUSE–NFM sample confirms similar candidates over comparable magnitude, separation, and redshift ranges, but increases the number of spectroscopically classified systems, enabling a more robust characterization of the population properties.

Compared to varstrometry-based selections (e.g., VODKA and VADAR; Shen et al. 2019; Hwang et al. 2020; Schwartzman et al. 2024, 2025; Gross et al. 2025), the GMP method relies on a different observational signature in Gaia data. Varstrometry is sensitive to unresolved systems through astrometric jitter due to flux variability (Hwang et al. 2020), whereas GMP identifies multiple peaks in the Gaia light profile, consistent with multiple components (Mannucci et al. 2022). Both techniques primarily select unobscured, optically bright quasars and, in practice, probe broadly similar regimes in redshift and physical separation. GMP candidates are generally found at z ≳ 0.5 and at kpc-scale separations. Varstrometry, in principle, is sensitive to even smaller separations and lower-redshift systems (Hwang et al. 2020); however, current observational results indicate that confirmed dual AGN from varstrometry largely occupy similar redshift ranges and comparable – or in some cases larger – projected separations than those identified by GMP. Available data suggest that the contamination rates reported in this work and by Gross et al. (2025) are comparable for both techniques, with star–quasar projections affecting each sample at similar levels. Notably, only two sources from the GMP sample also are varstrometry candidates. On the other hand, all currently known varstrometry-selected systems at 0.5 ≲ z ≲ 3.5 are also identified by GMP, provided the secondary AGN has G < 20.5, which is the completeness limit of the GMP selection Mannucci et al. (2023). Finally, the redshift distribution of the dual AGN identified here is broadly consistent with expectations from clustered quasar surveys (e.g., Hennawi et al. 2006, 2010), which find an increasing incidence of close quasar pairs at z ≳ 1, although typically at larger projected separations.

7. Conclusions

We have presented the first-year results of our MUSE Large Program targeting dual and lensed AGN candidates selected via the GMP technique (Mannucci et al. 2022). The first observed dataset comprises 30 targets spanning redshifts 0.7 ≲ z ≲ 3.2 and projected separations of ∼0.24″ − 0.82″, corresponding to ∼1.5 − 7 kpc at typical redshifts. Our observations have led to the discovery of six new dual AGN systems in the redshift range z ∼ 0.9 − 2.5 and with projected separations of < 7 kpc as well as 13 lensed AGN (including three quadruply imaged systems) in the redshift range z ∼ 0.7 − 3.2.

Based on literature compilations (e.g., Mannucci et al. 2022; Chen et al. 2022; Pfeifle et al. 2025), only ∼27 dual AGN are currently confirmed at such close separations across the redshift range 0.5 ≲ z ≲ 3.5, including the six systems presented here. The dual AGN in this work therefore represent ∼22% of the known population in this separation and redshift regime.

Combining this work with previous GMP-based confirmations (Mannucci et al. 2022, 2023; Ciurlo et al. 2023; Scialpi et al. 2024), 15 out of the 27 currently known dual AGN at 0.5 < z < 3.5 with separations < 7 kpc were originally identified through GMP selection (56% of the sample). A larger number of systems (18 out of 27) satisfy the GMP criteria, although three of these were first observed or identified through other selection methods (Junkkarinen et al. 2001; Schechter et al. 2017; Chen et al. 2023b). A similar consideration applies to strongly lensed AGN at subarcsec separations. While several hundred lensed quasars are known overall, only a small fraction have image separations below ∼0.7″ with spatially resolved spectroscopic characterization. The 13 lensed systems identified here therefore represent a substantial addition to the sample of compact-separation lensed AGN at cosmic noon accessible to IFU-based studies.

Using the MUSE data, we systematically classified all targets. The six dual AGN exhibit clear spectral differences between the two nuclei in emission-line profiles and line ratios, with projected separations of 3−6 kpc and relative velocities up to ∼1200 km s⁻¹. Lensed AGN were identified through flux ratios, spectral similarity, and the detection of lensing features in the MUSE datacube. Eleven systems correspond to AGN+star projections, reflecting a low contamination fraction and the reliability of the GMP selection.

We also systematically measured redshifts for both intrinsic and intervening NALs in all spectra. Intervening NALs provide a unique probe of the small-scale structure of the CGM and IGM along closely spaced sight lines (≲10 kpc), tracing the distribution, kinematics, metallicity, and ionization state of cold gas clumps (Mandelker et al. 2019; Dutta et al. 2024). Intrinsic NALs instead offer diagnostics of AGN outflows and dual activity (e.g., Hamann et al. 2012).

The distributions of targets in angular separation, J-band and GaiaG magnitudes, and redshift highlight the effectiveness of the GMP selection and the capabilities of the MUSE AO system. Lensed systems are predominantly detected at the bright end of the J-band distribution, while dual AGN are increasingly identified at fainter magnitudes.

This first-year dataset demonstrates the power of combining GMP preselection with MUSE IFU spectroscopy to systematically identify and characterize dual and lensed AGN at subarcsec separations. As an optical Gaia-based selection, the GMP method is primarily sensitive to unobscured, optically bright AGN. Spatially resolved spectroscopy allows us to disentangle close pairs, measure redshifts and line properties, and detect both intrinsic and intervening absorption systems. The full MUSE Large Program, which will expand the sample to ∼150 targets, will enable statistically robust studies of dual AGN occurrence as well as their physical properties, including black hole masses and luminosities, in relation to separation, redshift, and magnitude. Future analyses will include detailed modeling of outflows, absorber kinematics, and lensing configurations, providing unprecedented insights into the dynamics and environments of AGN at cosmic noon.

In summary, the combination of GMP selection and MUSE IFU observations offers a highly efficient strategy for exploring the sub-arcsec AGN population. This approach bridges the gap between photometric identification and detailed spectroscopic characterization, and it establishes a solid foundation for future large-scale studies of dual and strongly lensed AGN.

Data availability

The observational data used in this work are available from the European Southern Observatory archive: https://archive.eso.org

Acknowledgments

This publication was produced while attending the PhD program 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”. Based on observations collected at the European Southern Observatory (ESO) programs 109.22WA, 110.23SM, 112.25ET, and 114.27BY with VLT/MUSE, and 109.22W4, 112.25CT, 113.26TB, 113.26D7, 114.27FC, and 115.28CQ with NTT/EFOSC2 and VLT/FORS2. The authors sincerely thank the ESO staff in Garching and at Paranal for their availability, patience, and support during both the preparation and execution of the observations. In particular, we are grateful to the ESO instrument experts for their guidance on the MUSE AO system, which was invaluable in planning the observing strategy and ensuring optimal data quality. Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. 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. We thank Elena Pancino for her valuable support regarding the Gaia Archive and properties. We acknowledge financial contribution from INAF Large Grant “Dual AGN and binary supermassive black holes in the multimessenger 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”, from INAF Large Grant “The Quest for dual and binary massive black holes in the gravitational wave era”, (Bando Ricerca Fondamentale INAF 2024), the MIRACLE INAF 2024 GO grant “A JWST/MIRI MIRACLE: Mid-IR Activity of Circumnuclear Line Emission” We also acknowledge financial support by the PRIN-MUR project “PROMETEUS” 202223XDPZM financed by the European Union – Next Generation EU, Mission 4 Component 1, CUP B53D23004750006 and C53D2300080006. AF and EB acknowledge financial support from the Ricerca Fondamentale INAF 2024 under project 1.05.24.07.01 MiniGrant RSN1. MP acknowledges support through the grants PID2021-127718NB-I00, PID2024-159902NA-I00, and RYC2023-044853-I, funded by the Spain Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and El Fondo Social Europeo Plus FSE+. SC and GV acknowledge support from European Union’s HE ERC Starting Grant No. 101040227 – WINGS. We acknowledge financial support from the ASI-INAF agreements n. 2024-36-HH.1-2025 and 2025-29-HH.0. Finally, we thank the referees for their helpful comments and careful review, which improved the clarity of the paper.

References

Appendix A: Spectra of dual and lensed AGN systems

In this section we present the spatially resolved spectra of the full sample of systems composed of at least two AGN. The spectra are extracted from the circular apertures shown in Fig. 2.

Figure A.1 shows the spectra of the dual AGN systems, while Fig. A.2 presents the spectra of the doubly imaged lensed AGN systems. In both cases, the brightest AGN component is shown in purple, while the faintest component is shown in light blue.

Fig. A.1. Final MUSE spectra of the dual AGN systems, with target names and redshifts indicated in each panel. The spectra were extracted from the circular apertures shown in Fig. 2 (same color-coding), corrected for Galactic extinction, normalized to Gaia magnitudes, and shifted to the rest frame. For clarity, the spectra of the primary AGN (component A) have been re-normalized to facilitate comparison with those of component B; the applied renormalization factor is indicated in each panel. Vertical dotted purple lines mark the expected positions of prominent emission lines at the redshift of AGN A. The gap around λobs = 6000 Å (observed frame) corresponds to the New Generation Controllers (NGC) used in AO observations.

Fig. A.2. Final MUSE spectra of the double lensed AGN systems, with target names and redshifts indicated in each panel. The spectra were extracted from the circular apertures shown in the white-light images of Fig. 2 (same color-coding), corrected for Galactic extinction, normalized to Gaia magnitudes, and shifted to the rest frame. Vertical dotted purple lines mark the expected positions of prominent emission lines at the redshift of AGN A. The gap around λobs = 6000 Å corresponds to the NGC used in AO observations.

Finally, Fig. A.3 shows the quadruply imaged lensed AGN systems. The third and fourth images are shown in orange and green, respectively. For the system J0957–2242, the spectrum of the foreground lensing galaxy is also shown in the bottom-right panel.

Fig. A.3. Final MUSE spectra for the three quadruply imaged AGN systems. For each system, the four AGN images are shown in different colors, following the same color-coding adopted in Fig. 3. The top row shows the first two quadruple systems, while the bottom-left panel shows the four AGN spectra of J0957. All spectra are corrected for Galactic extinction, normalized to Gaia photometry, and shifted to the rest frame. Vertical blue dashed lines mark the wavelengths of the main emission lines. Since J0957 data also allowed us to identify the foreground lensing galaxy, its spectrum (component E) is shown separately in red in the bottom-right panel. For reference, the AGN spectrum that contaminates the lens is shown in gray, while the emission lines of the lensing galaxy are indicated by black vertical dotted lines.

Appendix B: Lensed AGN distribution in our sample

To investigate the relative prevalence of lensed and dual AGN as a function of redshift, we define the fractions relative to the total number of AGN in each redshift bin:

$$\begin{array}{c}\hfill F_{\mathrm{lens}}=\frac{N_{\mathrm{lensed}}}{N_{\mathrm{AGN}}},F_{\mathrm{dual}}=\frac{N_{\mathrm{dual}}}{N_{\mathrm{AGN}}},\end{array}$$

where N_AGN = N_lensed + N_dual. These fractions describe the relative contributions within our observed sample and are useful to estimate the expected level of lensing "contamination" in dual-AGN studies. We stress that they should not be interpreted as intrinsic lens fractions suitable for lensing statistics, but only as indicators for selection purposes.

Figure B.1 shows F_lens in broad redshift bins (Δz = 0.5). We fit the redshift dependence of F_lens using both a weighted quadratic relation and a linear model, in order to provide a simple empirical description of the observed trend. The quadratic model is F_lens = az² + bz + c, with best-fit coefficients a = −0.195, b = 0.875, and c = −0.229, obtained using inverse variance weighting (σ⁻²) to account for larger statistical uncertainties at high and low redshifts. For comparison, we also derive a linear fit of the form F_lens = mz + q, with m = 0.259 and q = 0.228. The quadratic relation provides a better description of the data than the linear model in terms of χ², but both parameterizations are shown to guide the eye and to facilitate comparison with future samples. We note that performing this fit is useful for future work, as it allowed us to estimate the expected dual-AGN fraction among GMP-selected systems based on the measured lensing fraction, providing a reference for selection effects in our sample. The resulting quadratic relation, F_lens = −0.195 z² + 0.875 z − 0.229, suggests a broad maximum in the lensing fraction at intermediate redshift, corresponding to the redshift range where survey selection effects and the quasar luminosity function are most favorably aligned.

Fig. B.1. Fraction of lensed AGN within the lensed+dual sample, Flens = Nlensed/(Nlensed + Ndual), as a function of redshift z. Data are binned in intervals of Δz = 0.5. Error bars represent 1σ binomial uncertainties. The solid black line shows the weighted quadratic fit, while the dashed line indicates the linear fit.

Appendix C: AGN+star system

In this appendix we present the spectra of systems identified as fortuitous alignments of an AGN and a star along the same line of sight. We include the 11 systems observed in this Large Program, as well as an additional system, J1649+0812 (16:49:41.30 +08:12:33.5) with a projected separation of 0.59″, observed in the MUSE program (ID: 109.22W5, PI: Mannucci) and previously classified by Chen et al. (2025).

Fig. C.1. Final MUSE spectra of the AGN (purple) and star (orange) systems, with target names and redshifts indicated in each panel. The spectra were extracted from circular apertures of 5 px, corrected for Galactic extinction, normalized to Gaia magnitudes, and shifted to the rest frame. Vertical dotted lines indicate the positions of the expected AGN emission lines. The gap around 6000 Å corresponds to the NGC used in AO observations.

All Tables

All Figures

	Fig. 1. Distribution of dual (green) and lensed (blue) AGN as a function of projected separation (top panel) and redshift (bottom panel). The sample, observed with MUSE–NFM, includes sources presented in this work and by Scialpi et al. (2024). For quads, the Einstein radius was used to characterize the separation.
In the text

Fig. 2. White-light images of the double AGN systems (lensed or dual) in units of 10−15 erg/s/cm−2. The target name and the separation between the two components are indicated in each panel. The spectra are extracted from the circular apertures shown on the maps, with purple and blue apertures corresponding to the A (brightest) and B (faintest) components, respectively. Dual systems are marked with the flag D, lensed systems with L. For systems identified as gravitational lenses, if an additional lensed image is detected, it is shown inside the red aperture and labeled as component G.

In the text

Fig. 3. White-light images of the three quadruply lensed AGN systems. The target name and the value of the Einstein radius (Re) is indicated in each panel. The spectra are extracted from the colored circular apertures. The colors of the apertures follow the order of brightness: A (purple, brightest), B (blue), C (orange), and D (green, faintest). For J0957, we also detect the emission-line lensing source (E) in red.

In the text

	Fig. 4. Fraction of AGN+star projections as a function of projected separation. The dotted red line shows the expected contamination probability function (Eq. 3 with D = 0.65″), representing the separation distribution of random stellar contaminants. The shaded area shows the separation range excluded by the GMP selection.
In the text

	Fig. 5. Absorption systems of the dual AGN J0032–1053. The upper panel shows component A, while the bottom panel shows component B. The two intervening absorption systems are indicated: the system at z = 1.079 is shown in dark red and the system at z = 1.479 in teal. Intrinsic C IV absorption at the AGN redshifts is shown in blue.
In the text

	Fig. 6. Enlargement of the Mg II spectral region of J2204–6530 (AGN A in the top panel and AGN B in the bottom one). Both AGN show the Mg II doublet absorption lines at rest-frame wavelengths of 2796 Å and 2803 Å. In each panel, the redshift of each QSO and the NAL velocity relative to its systemic redshift are reported, assuming that the two outflows are independent.
In the text

Fig. 7. BPT diagram adapted from Lamareille (2010) showing the positions of the local SDSS galaxies. In the left panel, the yellow strip shows the range of [N II]λ6584/Hα values corresponding to the observed [O III]/Hβ ratio (the dashed red area indicates the measurement uncertainty). In the right panel, the yellow star marks the position of the lens galaxy in the [O III]/Hβ versus [O II]/Hβ diagram, which lies very close to the separation line between star-forming galaxies (blue) and AGN (green and cyan).

In the text

	Fig. 8. MOKA3D rotating disk best-fit model of J0957−2242 traced by the [O III]5007 transition. Panels show the observed (top), best-fit (center), and residual (bottom) moment maps. The residual maps were obtained by subtracting the model from the observed moment maps. The star marks the kinematic center, inferred as the peak flux emission. Pixels are masked at SN ≤ 4.
In the text

Fig. 9. Top panel: Fraction of lensed systems (blue) and dual AGN (green) as a function of J-band magnitude for the MUSE sample (this work and Scialpi et al. 2024). Bottom panel: Same but as a function of redshift. Data are binned in intervals of 0.5 mag (top) and 0.5 (bottom). The number above each column indicates the total number of objects in that bin. The horizontal dashed line marks the 50% level for reference.

In the text

Fig. A.1. Final MUSE spectra of the dual AGN systems, with target names and redshifts indicated in each panel. The spectra were extracted from the circular apertures shown in Fig. 2 (same color-coding), corrected for Galactic extinction, normalized to Gaia magnitudes, and shifted to the rest frame. For clarity, the spectra of the primary AGN (component A) have been re-normalized to facilitate comparison with those of component B; the applied renormalization factor is indicated in each panel. Vertical dotted purple lines mark the expected positions of prominent emission lines at the redshift of AGN A. The gap around λobs = 6000 Å (observed frame) corresponds to the New Generation Controllers (NGC) used in AO observations.

In the text

Fig. A.2. Final MUSE spectra of the double lensed AGN systems, with target names and redshifts indicated in each panel. The spectra were extracted from the circular apertures shown in the white-light images of Fig. 2 (same color-coding), corrected for Galactic extinction, normalized to Gaia magnitudes, and shifted to the rest frame. Vertical dotted purple lines mark the expected positions of prominent emission lines at the redshift of AGN A. The gap around λobs = 6000 Å corresponds to the NGC used in AO observations.

In the text

Fig. A.3. Final MUSE spectra for the three quadruply imaged AGN systems. For each system, the four AGN images are shown in different colors, following the same color-coding adopted in Fig. 3. The top row shows the first two quadruple systems, while the bottom-left panel shows the four AGN spectra of J0957. All spectra are corrected for Galactic extinction, normalized to Gaia photometry, and shifted to the rest frame. Vertical blue dashed lines mark the wavelengths of the main emission lines. Since J0957 data also allowed us to identify the foreground lensing galaxy, its spectrum (component E) is shown separately in red in the bottom-right panel. For reference, the AGN spectrum that contaminates the lens is shown in gray, while the emission lines of the lensing galaxy are indicated by black vertical dotted lines.

In the text

	Fig. B.1. Fraction of lensed AGN within the lensed+dual sample, Flens = Nlensed/(Nlensed + Ndual), as a function of redshift z. Data are binned in intervals of Δz = 0.5. Error bars represent 1σ binomial uncertainties. The solid black line shows the weighted quadratic fit, while the dashed line indicates the linear fit.
In the text

Fig. C.1. Final MUSE spectra of the AGN (purple) and star (orange) systems, with target names and redshifts indicated in each panel. The spectra were extracted from circular apertures of 5 px, corrected for Galactic extinction, normalized to Gaia magnitudes, and shifted to the rest frame. Vertical dotted lines indicate the positions of the expected AGN emission lines. The gap around 6000 Å corresponds to the NGC used in AO observations.

In the text