Optical Spectroscopy of Dual Quasar Candidates from the Subaru HSC-SSP program

The Sloan Digital Sky Survey (SDSS) uses a 2.5 m wide-angle optical telescope at Apache Point Observatory (APO) in New Mexico. Since the start of operations in 2000, its spectral database has grown to four million objects covering over 14,055 square degrees of sky (up to DR14; Abolfathi et al. 2018; Pâris et al. 2018). DR14 is the second data release of the fourth phase of the Sloan Digital Sky Survey (SDSS-IV). It provides a quasar catalog from the extended Baryon Oscillation Spectroscopic Survey (eBOSS), that also includes previously spectroscopically confirmed quasars from SDSS-I, II (Schneider et al. 2010), and III (Eisenstein et al. 2011; Pâris et al. 2012, 2017). The SDSS quasar selection includes magnitude- and color-selected samples and objects identified through radio, infrared, and X-ray surveys including FIRST, GALEX, 2MASS, and ROSAT. In total, the catalog contains 526,356 quasars up to z ∼ 5. The spectral data are taken with fibers, each subtending 3'' (SDSS-I/II) or 2'' (BOSS/eBOSS) diameter on the sky, connected to plates that each have a field of view of approximately 7 deg² and connecting 640 (SDSS-I/II) or 1000 (BOSS/eBOSS) fibers. The earlier spectroscopic data (SDSS-I/II) covers wavelengths from 3800 to 9200 Å at a spectral resolution of ∼2000 (Schneider et al. 2010). The later data (BOSS/eBOSS) cover 3600–10,400 Å with resolution varying from ∼1300 at 3600 Å to ∼2500 at 10,000 Å (Pâris et al. 2018).

Hyper Suprime-Cam (HSC; Miyazaki et al. 2018) is a wide-field (17 diameter) optical imager installed at the prime focus of the Subaru Telescope. The Subaru Strategic Program (SSP; Aihara et al. 2018) is a survey using this instrument that began in 2014 March with the second public data release in 2019 (PDR2; Aihara et al. 2019). The survey includes three layers of imaging data: Wide, Deep, and UltraDeep. The imaging data used in this study is from the Wide layer with 796 deg² of i-band imaging having at least one 200 s exposure, with median seeing of 058 and depth 26.2 mag for 5σ detection. Our dual quasar search is first based on the i band. We then include the images taken in the other four broadband filters (i.e., g, r, z, and y; Kawanomoto et al. 2018) to provide color information and stellar mass estimation for objects at z ≲ 1.

Since HSC does not yet have a sample of spectroscopically confirmed quasars based on optical selection, we chose to start our search for dual quasars using the catalog of spectroscopically confirmed SDSS quasars (Pâris et al. 2018). This catalog provides additional information such as emission line properties that can be used to estimate black hole masses. To search for dual quasars, the SDSS DR14 quasar catalog ensures that at least one has already been identified. Our idea is to identify those that have one (or more) close companions that lack spectroscopic data because of the limitation of the SDSS spatial resolution (∼15), fiber collisions (∼55''), fiber diameter (2'' and 3''), or faintness of the optical companion (22.2 mag for i-band 5σ depth).

We match the SDSS quasar catalog with the HSC photometric catalog to make use of the higher resolution and deeper optical imaging available through the HSC SSP. There are 34,476 SDSS quasars imaged by HSC that are unsaturated (i_AB ≥ 18), have no bad pixels, and have model magnitudes available in HSC. HSC image cutouts of size 60^'' ×60^'' are generated for each object in each band, together with the variance image and model point-spread function (PSF). As described in Silverman et al. (2020), we run a Lenstronomy-based (Birrer & Amara 2018) modeling algorithm on the images with a composition of point sources and Sérsic profiles. We then pick out candidates that appear to have two or more point sources within the field, with separation ranging from 06 to 3^''. After filtering out the artifacts such as satellites and known lenses (Inada et al. 2008, 2010, 2012), there are 421 candidates left, of which 401 have photometry in all five HSC bands. The PSF positions estimated by the algorithm are given in Columns 2–5 in Table 1, with index 1 referring to the known SDSS quasars and index 2 referring to the unknown sources. The algorithm also measures the HSC PSF magnitudes of each sources (Columns 8–13 in Table 1), required to determine which companion objects have colors consist with quasars, thus removing stellar contamination. Figure 2 shows the distribution of bolometric luminosity and redshift for the 401 dual quasar candidates and the parent population of SDSS quasars having HSC imaging. The measurement of bolometric luminosity is described in Section 6.1. We further separate our candidates into blue (HSC g − r ≤ 1) and red (HSC g − r > 1) companions as shown in the figure. We give higher priority to the blue sources in the color–color diagram to remove objects that are likely to be stars (see Section 5.1 for details) when designing our observations.

Zoom In Zoom Out Reset image size

We observed twenty candidates on 18.09.2019 and 19.09.2019 using the Faint Object Camera and Spectrograph (FOCAS; Kashikawa et al. 2002) on the 8.2 m Subaru telescope (S19B-079, PI: J. Silverman). Nine targets were observed on the first night with the grism VPH850+O58, which covers the wavelength range from 5800 Å to 10350 Å and has a spectral dispersion of 1.17 Å pixel⁻¹. During the second night, we observed eleven targets with the 300B+Y47 grating, which covers wavelengths from 4700 to 9100 Å with a dispersion of 1.34 Å pixel⁻¹. The slit width was 08 for both nights. Detailed information is listed in Table 2.

Data reduction was carried out with the focasred package, following the Subaru Data Reduction CookBook. ²⁶ We started with bias subtraction and flat-field correction. We did not correct for the spatial distortion, to avoid introducing additional noise. Then we applied L.A. COSMIC, which is a cosmic ray removal tool based on a variation of Laplacian edge detection developed by Van Dokkum (2001) to remove cosmic ray events. Subsequently, we implemented the IRAF apall task in apextract package to extract the trace of the 1D spectrum. Wavelength calibration was performed using ThAr lamp exposures. Finally, we observed the standard star G191B2B with the same setup as the science exposures for flux calibration.

We observed nine candidates with the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2004) on the 8.1 m Gemini-North telescope (GN-2019B-Q-128, PI: J Silverman) from 2019 September to 2019 December in the queue mode. Depending on the redshift of the dual quasar candidate, a different setup including the choice of grating was configured to cover a wavelength range suitable to detect broad emission lines (e.g., C iv, Mg ii, Hα) characteristic of quasar activity (Table 1).

The data were taken with the GMOS-N Hamamatsu detectors. CuAr lines were used for wavelength calibration. The observed standard stars were EG131 and Feige 66. The reduction was done using the semi-automated Python Spectroscopic Data Reduction Pipeline (PypeIt; Prochaska et al. 2020). Parameters were tuned according to the documents, to fit the usage of each object.

Considering the different algorithms used in IRAF and PypeIt, we normalize the results by rescaling the flux calibrated spectra to match the PSF magnitudes of the sources measured from our image decomposition. These HSC photometric data points are plotted as star symbols in 1D spectra shown in Section 5.2, Appendix B.1, and B.2.

To confirm a dual quasar system, we classify the companion based on its optical spectrum. Generally, when broad emission lines such as C iv, CIII], Mg II, Hβ, or Hα are observed, we classify the companion as a type 1 quasar, while diagnostic diagrams ("Baldwin, Phillips, & Terlevich" (BPT) diagrams (Baldwin et al. 1981) and WHAN diagrams (Cid Fernandes et al. 2011)) are used to classify type 2–like sources. Broad absorption lines are also common among quasars (e.g., Hewett & Foltz 2003); for example, C iv may form a so-called "P-Cygni" profile, with a broad absorption feature impacting the blue wing of the core emission. There is one such case that we classify as a "BAL quasar." Table 3 lists the six dual quasars that we have confirmed to date. Classifications of the rest of the targets are reported in Appendix B. The details of the emission-line measurements are given in Section 4.

We show the completeness of this work in Figure 3 with stacked histograms binned by redshift and i-band magnitude difference. The colors refer to the status of the targets, either observed (gray for the failures and green for the successes) or not observed (separated by the companions' g − r color). The upper panel gives a statistical view of Figure 2 on the x-axis. Our follow-up observations covered redshift from 0 to ∼3.2, including 20 blue targets and 12 red targets. The i-band magnitude differences of the pairs correspond to flux ratios from 0.03 to 1. Therefore, our selection algorithm is able to cover pairs that differ in flux by up to a factor of ∼30.

Zoom In Zoom Out Reset image size

Out of the 32 targets we followed up spectroscopically, we confirm six quasar pairs, giving a success rate of ∼19% for finding physically associated quasar pairs. This success rate is significantly lowered by the inclusion of candidates with colors in agreement with the stellar locus, which are necessary to cover those red quasars. Below, we discuss means of limiting our rate of failure. In addition, we identified three AGN–galaxy pairs (only one SMBH in the pair is active) and two project quasars pairs in which the quasars are at significantly different redshift with velocity offset larger than 2000 km s⁻¹ (Hennawi et al. 2006).

Figure 4 provides a view of these results in color–color space, with only the companions plotted. In the left panel, the blue curve is a stellar locus adopted from Covey et al. (2007), which shows the average distribution of stars as a reference. The gray dots are our 401 candidates. The magnitudes used here are the PSF magnitudes estimated from our decomposition methods based on HSC imaging (Ding et al. 2020; Li et al. 2021). We claim that the measurements of AGN magnitudes are precise since the light distribution of the point sources are quite unique; the uncertainties are under 0.02 mag within a certain model (e.g., 2 Point Sources + 1 Sérsic Profile). In some cases, the spectrum is too noisy or too featureless to tell which class they belong to; these unclassified targets are labeled with a purple cross.

Zoom In Zoom Out Reset image size

For comparison, we show the kernel density estimation (KDE) plot for normal HSC stars, galaxies, and quasars, respectively. The galaxies are selected with r_extendness = 1 and i_extendness = 1, which means the HSC Pipeline (Bosch et al. 2018) classified them as extended sources in both the r and i bands. Stars are selected as point sources in these two bands with r_extendness = 0 and i_extendness = 0. While this sample will include some quasars, stars outnumber quasars by several orders of magnitude, and it is safe to consider it as a statistical star sample. Quasars are selected from our parent sample pool, i.e., the HSC-imaged SDSS quasars (Section 2.3). We reject objects with g − r outside the 1%–99% range of the distribution for the stars, galaxies, and quasars, to remove objects with poor photometry. Finally, ∼4000 targets are randomly selected from each of the categories as the comparison sample for Figure 4.

We find that the classification of our candidates based on spectroscopy is consistent with the KDE distribution of the comparison categories. Since the comparison quasar sample is from SDSS, they are mostly type 1 quasars, and thus are concentrated at the lower left region of the r − i versus g − r plot. As expected, this is also the region in color space where we have a high success rate. However, contamination from galaxies and stars is still significant. We represent the figure as a distance–color space (Figure 4 right panel), where distance is defined as the offset of our targets from the stellar locus in the r − i versus g − r space (left panel). This separates the different populations, to a certain extent. The quasars are systematically further away from the stellar locus than the stars, and are bluer than the galaxies. We draw a vertical dashed line at g − r = 1.0 and a horizontal line at "distance from stellar locus=0.06," dividing the figure into four regions. We have a high success rate of finding quasars in the upper left region (7 out of 9). One of the contaminants is a blue star, while another one is a narrow-line galaxy with emission-line ratios that suggest it is a composite AGN and star-forming object.

The lower region mostly contains stars, since they indicate a very close distance to the stellar locus. However, one true quasar pair at z ∼ 3.1 lies in the lower region. Its host galaxy is undetected, making the source star-like, thus it lies very close to the stellar locus.

Galaxies are more widely distributed in this diagram, as shown by the comparison samples. The diffused galaxies are filtered out by our selection algorithm, leaving those with luminous bulges or star formation regions that might appear as fake point sources.

Hennawi et al. (2010) studied the colors of dual quasars to z ∼ 5.5 (see their Figure 4). High-redshift quasars (z > 3) are more scattered in the color–color space than the low-redshift quasars. And due to their relatively low number density, they are not recognizable in the contours of Figure 4. One has to consider redshift as well as the color when pursuing high-redshift quasars.

We now focus on the color of the quasar pairs. We plot both sources of each pair on Figure 5. In the upper panel, red dots are the same data set as in Figure 4, and black dots are the parent sample of these systems, i.e., the known SDSS quasars. Our eight quasar pairs are shown with star symbols, where the color indicates the two sources in a given system. For example, the blue and cyan star symbols denote the two projected quasar pairs (SDSS J0213-0221 and J2251+0016).

Zoom In Zoom Out Reset image size

Richards et al. (2001) show that the locus of quasar colors as a function of redshift is very tight. The lower panels plot this locus; the shadowed region contains 95% of the sample at any given redshift. We overlay our physically associated and projected quasar pairs. Eight of the 12 quasars in physically associated pairs fall in the 95% region of g − r, 9 out of 12 fall in the 95% region of r − i, but most of our objects lie above the mean color. Thus, dual quasars have colors that are not drawn exclusively from the color distribution of SDSS-selected quasars. The fact that they tend to be redder suggests that they are affected by dust in their hosts, due to the interaction of the two objects (Richards et al. 2003).

The two projected quasar pairs are not interacting with one another, suggesting that they should lie in the 95% region. This is true for J2251+0016, but not for J0213-0211. The redshift separation between the two objects in the latter source is 6000 km s⁻¹, so it is conceivable that they are physically associated with one another.

It is also worth mentioning that the redshifts of the sources are themselves not precisely determined in all cases: the emission lines are broad and often asymmetric, so this can affect the measurement of the apparent line-of-sight velocity offset. The typical shift of C iv is ∼800 km s⁻¹ compared to Mg ii, while the latter usually only have a small blueshift compared to [O iii] (Richards et al. 2002).

In the following section, we report the details of three newly identified, physically associated quasar pairs. The projected pairs are presented in B.1, and the quasar–galaxy pairs in B.2.

We report on the properties of three newly discovered dual quasars that are physically associated, which is defined as a quasar pair with line-of-sight velocity offset smaller than 2000 km s⁻¹ (Hennawi et al. 2006). The lens hypothesis is ruled out in each case by the fact that the spectra of the two objects in each pair are not identical.

5.2.1. SDSS J220906.91 + 004543.9

The SDSS survey reported this object as a broad-line quasar at z = 0.4466. The Subaru HSC images are shown in Figure 6(a). The color image (RGB), shown in the left panel, is generated using the g, r, and i bands (Lupton et al. 2004). The magenta circle shows the position of the SDSS fiber, whose diameter is 3'' for this target. The green dashed rectangular region marks the position of the long slit used for spectroscopy. The red and blue rectangular regions are the extraction apertures on the 2D spectrum for the main SDSS quasar (source A) and the companion (source B), respectively. We keep this format of indices for all the sources reported later in this paper. The four figures on the right-hand side show the image decomposition results using the i-band image and Lenstronomy-based algorithm (Birrer & Amara 2018). From left to right, the panels are the following: (1) original HSC data, (2) model reconstructed with point sources + host galaxies, (3) host galaxies of the system (data-point source models), and (4) residual of data-model normalized by the pixel errors. According to the decomposition results, the i-band magnitudes of the point sources are 19.92 and 20.25, respectively. We carry out a similar analysis for all five bands and find that the southern source is redder (g − r = 0.7) than the typical SDSS quasar (Figure 5). The same is true for the primary quasar. We also measured the separation of the two point sources. The point sources are separated by 167, which corresponds to 9.55 kpc at z = 0.4466.

Zoom In Zoom Out Reset image size

The optical spectra are presented in the middle panels of Figure 6(b). The 1D spectra are extracted from the rectangular regions shown in the HSC color image. We corrected the flux-calibrated spectrum by matching to the HSC photometry, which are plotted as stars, dark red for source A, and dark blue for source B. The position of emission lines are also labeled with dashed vertical lines having the same coding. The regions affected by telluric oxygen and water vapor absorption are shaded in gray (oxygen bands: 6280–6340 Å, 6860–6950 Å, 7590–7720 Å; water vapor bands: 9270–9380 Å). We show a zoom-in view of the 2D spectrum with the corresponding sources (A and B) labeled on the left. The traces are well-separated. The FOCAS spectrum of this target covers the Hβ and Hα regions.

The bottom panel (Figure 6(c)) provides supplemental and detailed information of the pair. Line-fitting results of the two sources are given with the first (second) row corresponding to source A (B). We fit the Hβ and Hα regions using the methods described in Section 4. Vertical dashed lines indicate the centroid of each emission line. The tables on the right-hand side record the widths of narrow and broad lines with errors. The offsets of the line centroids between the two sources or between the broad and narrow components of the same line are listed. For example, "An_Ab" stands for the offset between the narrow line and the broad line of source A. "Ab_Bb" stands for the offset between the broad line of source A and the same broad line of source B. The velocity is positive when the former line is redshifted relative to the latter line. For SDSS J2209+0045, the measurement of the narrow lines indicate very small offsets, and both sources have broad emission lines > 2000 km s⁻¹. The difference of the line ratios tells us that this is not a lens system, and we thus classify it as a quasar pair.

In panel (a), the decomposition of the i-band image resolves the host galaxies of this system. The model indicates a Sérsic index of 0.8 for the main quasar and 3.8 for the companion. After subtracting the two Sérsic profiles and point-source components from the science image, the residual image seems to indicate an asymmetry of the galaxies' morphology, possibly resulting from the ongoing merger event.

5.2.2. SDSS J233713.66 + 005610.8

SDSS classified this object as a broad-line quasar at z = 0.7080. The fiber is centered on the southern source as shown in Figure 7(a). The color of the two sources are very similar, with g − r typical of quasars and r − i being redder by 0.3 mag (Figure 5 bottom panels). The HSC image also shows a third source to the north of the pair, thus we fit 2 PSFs + 3 Sérsic profiles for this system.

Zoom In Zoom Out Reset image size

The GMOS spectrum of this target was taken with the R831+RG610 grism, and the Hβ region was placed at the center of the spectrum (Figure 7(b)). The traces of the two objects in the 2D spectrum are well-separated and we can clearly see the [O iii] and Hβ features in both spectra. Notably, the [O iii] doublet appears to have blue wings in both sources. To support this statement, we fit the lines with two different models (panel c). In the left column, no outflow components are used, thus each [O iii] line is modeled with a single Gaussian. This leaves a significant residual at the peak of [O iii] in source A and Hβ in source B. Note that the [O iii] and the narrow component of Hβ are constrained to have the same profile. When we include additional outflow components for the [O iii] doublet in the two sources (orange curves in the right column of Figure 7 panel (c)), the fit is significantly improved. Based on [O iii], we measured a redshift of 0.7089 for source A and 0.7083 for source B. The velocity offset of the narrow lines between two sources is 102 ± 7 km s⁻¹. This, together with the different line ratios, cause us to classify the target as a dual quasar system.

J2337 is the only target for which we have to fit an iron template to reconstruct the spectra. The relative strength of the iron emission in source B is higher than that in source A, while the equivalent width (EW) of the [O iii] core component (magenta components) and FWHM of Hβ are narrower in source B. This follows the trend expected by Eigenvector 1 (e.g., Shen & Ho 2014).

5.2.3. SDSS J123821.66 + 010518.6

The SDSS BOSS program reports this target as a broad-line quasar at z = 3.1327. Both sources are red in g − r color and lie close to the stellar locus (Figure 5). The primary source is two magnitudes brighter than the secondary. The host galaxies are undetected in this system in our image decomposition, which is unsurprising given the redshift of the pair. Only two PSFs are applied to the model.

The GMOS spectrum of this target was taken with the B600+CG455 grism. Lyα and C iv λ1549 were both observed (Figure 8(b)), and the latter was used to measure the redshifts of the sources. Considering that the asymmetry of the C iv λ1549 line might affect our measurement of the redshift, we supply the results via matching SDSS templates ²⁷ to our data (Figure 8(c) middle panels). We estimated the reduced χ² values from z = 0 to z = 5, which yields the best match of source A at z = 3.133 (template No.30) and z = 3.135 for source B (template No. 32). With a slight difference from these two measurements, the classification of a physically associated pair holds for both.

Zoom In Zoom Out Reset image size

Furthermore, we found that the companion source has the "P Cygni" profile typical of BAL quasars, while the primary source does not. We fit the profile of the companion with two Gaussian functions (Figure 8(c), left panels), and find that the BAL component is blueshifted from the emission line by 8303 ± 625 km s⁻¹ with an FWHM of 7742 ± 411 km s⁻¹, compared to an FWHM of 10026 ± 870 km s⁻¹ for the emission-line component. According to Weymann et al. (1991), a quasar with contiguous absorption extended over more than 2000 km s⁻¹ and with a blueward velocity of at least 5000 km s⁻¹ can be considered as a BAL quasar. Therefore, we classify this as a BAL quasar. Hamann et al. (1993) studied line profiles of 40 BAL quasars, and found that most BAL quasars have stronger absorption than the emission, which suggests that the BAL regions do not completely envelop the continuum source. For J1238B, the equivalent width of the C iv λ1549 emission line is 55 Å, smaller than 60 Å of the BAL. Also, there are multiple narrow absorption lines in the brighter source, while they are absent in the fainter source, which tells us that the absorption clouds are in front of the brighter source only.

This is the highest-redshift quasar pair at close projected separation (7.67 kpc) yet found. Fortuitously, the BAL nature of the companion firmly rules out the possibility of it being a lensed system, thus mitigating such complications (Shen et al.2021).

For type 1 quasar pairs, we estimate the BH masses for each source using the viral method (Vestergaard & Peterson 2006):

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}={10}^{\alpha }{L}^{\beta }{\left(\displaystyle \frac{\mathrm{FWHM}}{1000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{\gamma }{M}_{\odot },\end{eqnarray} \tag{ 1 }$$

where FWHM is the full width at half maximum of the specific broad emission lines used for the calculation, L is the corresponding monochromatic luminosity, and α, β, and γ are the indices calibrated from reverberation-mapping methods. Our measurements are based on C iv, Mg ii, Hβ, and Hα. The parameters used for each are listed in Table 4. For each observed line, the monochromatic luminosity is measured from the power-law fit to the continuum, and the FWHM is measured from Gaussian fit described in Section 4. We further note that the normalization of the spectra are set by the HSC photometry after running our 2D decomposition. Thus, the quasar luminosity is free of host galaxy contamination.

In Figure 9, we plot the black hole mass for the objects in our sample as a function of the spectroscopic redshift, and label both quasars in each dual quasar system with the same colored symbol. We show the median value for cases with multiple estimations of BH mass from different emission lines for a single object. We also provide the value of the mass ratio of the two, connected by the dashed lines in the figure. The measurement errors on the black hole masses are propagated from the uncertainties on the parameters as output from the line-fitting routine. Also, the virial estimators have an intrinsic scatter of 0.3 dex. The errors on the redshift measurements are negligible. For comparison, we use the values for the full SDSS quasar sample as provided by Rakshit et al. (2020) as indicated by the background sources with contours indicative of the number density of quasars. The comparison sample includes 525,103 quasars from SDSS DR14 with black hole mass measurements. The full quasar sample covers a redshift range from 0 to 4.0 and black hole masses between 10⁷ and 10¹⁰ M_⊙. Our dual quasars have masses comparable to single quasars and have black hole mass ratios within a factor of five (Figure 9).

Zoom In Zoom Out Reset image size

Figure 10 displays the bolometric luminosity and Eddington ratio of our dual quasars as a function of black hole mass. The log bolometric luminosity ranges from 44.2 to 47.5, while the log Eddington ratio ranges from −2.2 to 0. The bolometric luminosity L_bol in the upper panel of Figure 10 is calculated using monochromatic luminosity and corresponding bolometric correction factors provided by Richards et al. (2006): BC5100 = 9.26, BC3000 = 5.15, and BC1350 = 3.81. The Eddington luminosity L_Edd is estimated from the BH mass as:

$$\begin{eqnarray}&&{L}_{\mathrm{Edd}}\cong 1.3\times {10}^{38}\left({M}_{\mathrm{BH}}/{M}_{\odot }\right)\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 2 }$$

The Eddington ratio is then given by L_bol/L_Edd. We separate our quasar pairs into low-redshift (z < 1) and high-redshift (z > 3) pairs and compare them with the corresponding SDSS population in Figure 10. Panels (a) and (c) show the low-redshift cases, while (b) and (d) show the high-redshift objects.

Zoom In Zoom Out Reset image size

We find that the bolometric luminosity and Eddington ratio of both objects in the low-redshift quasar pairs do not have significant difference from single quasars at the same BH mass. This agrees with the results of Stemo et al. (2020). Figure 3 shows that we do not have a bias toward equally luminous pairs in follow-up observations. However, the selection algorithm itself could be biased toward candidates with stronger point-source features, and thus these candidates would be more likely to have properties similar to those of the single quasars when successfully confirmed. For the high-redshift pair (SDSS J1238+0105), the main source is luminous with a high Eddington ratio, while the fainter source has large uncertainties in the measurements, making it harder to draw a strong conclusion.

To highlight, we have so far not found any quasar pairs in the redshift range between 1 and 3 even with suitable candidates as shown in Figure 3. In this redshift range, our spectra are only able to cover one emission line, either C iv λ1549 or Mg ii λ2798, hence the identification of a weak emission line is difficult, especially in cases of type 2 quasars.

Our decomposition algorithm returns five-band photometric magnitudes of the host galaxies of our dual quasar systems. Based on this photometry, we carried out a stellar mass estimation for systems at z < 1 using CIGALE (Code Investigating GALaxy Emission; Boquien et al. 2019). Li et al. (2021) used CIGALE models for the study of quasar hosts with HSC. The computation of stellar mass includes an assumed star formation history (SFH), nebular emission, attenuation, dust emission (not included in our cases due to the lack of IR data), and AGN emission (not included in our cases because we have already subtracted the point sources) to generate a model SED based on a fit to the HSC photometry.

In Figure 11, we show the rest-frame U − V color of the host galaxies of our dual quasars as a function of the stellar mass. Recall that we modeled either one or two galaxy components for the pairs. In Figure 11, the star symbols represent the hosts of the main quasars, while the circular markers of the same color represent the companions if two separate galaxy components are identified. The U- and V-band magnitudes are measured from the best-fit SED model. A typical color uncertainty of 0.28 mag is applied to all the sources, and the uncertainties of the stellar masses are given by CIGALE. For comparison, we provide a sample set of 709,006 normal galaxies over the redshift range 0.2 and 1 with 10⁷ < M_* < 10¹³ (Kawinwanichakij et al., submitted). We find that the host galaxies of our quasar pairs are mostly massive, but are bluer than red-sequence galaxies. Most of them are located below the boundary line $U-V=0.16\times \mathrm{log}\ {M}_{* }+0.16$, which separates the star-forming and quiescent galaxies. Capelo et al. (2015) suggested in their simulation that the interaction of galaxies may trigger star formation in the merger stage (see their Figures 1 and 2). This also agrees with the observational findings of host galaxies of single quasars (Matsuoka et al. 2015; Ishino et al. 2020; Li et al. 2021).

Zoom In Zoom Out Reset image size

The correlation between the BH mass and mass of the host galaxy has been quantified by many studies, but the reason for this connection remains uncertain. Silk & Rees (1998), Fabian (1999), and King (2003), among others, suggest that mergers work as a global process to build this correlation. If this is true, we may expect to see that systems prior to merging are offset from the local correlation, with the merger bringing them onto the correlation. However, it is difficult to trace the full life cycle of a merger with observational snapshots. Therefore, we will first test this hypothesis from the view of the Horizon-AGN simulation.

Horizon-AGN is a cosmological hydrodynamical simulation with the size of the box L_box = 100h¹ Mpc. Dubois et al. (2014) measured the properties of their mock galaxies and studied the correlation between the spin directions and the large-scale cosmic web structure. Volonteri et al. (2016) included BHs and AGN in their simulation. BHs are created in gas and stellar overdensity regions that exceed the threshold of star formation, n₀ = 0.1 H cm⁻³, where the seed mass is assumed to be 10⁵ M_⊙. The accretion onto BHs follows the Bondi–Hoyle–Lyttleton (BHL) accretion rate:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}=\displaystyle \frac{{G}^{2}{M}_{\mathrm{BH}}^{2}{\rho }_{\mathrm{gas}}}{{\left({c}_{s}^{2}+{v}_{\mathrm{rel}}^{2}\right)}^{3/2}},\end{eqnarray} \tag{ 3 }$$

where M_BH is the BH mass, ρ_gas is the average gas density, c_s is the average sound speed, and v_rel is the average gas velocity relative to the BH velocity.

A BH is assigned to a halo when it is located within 0.1 × R_vir of the center (central BH) or within R_vir (off-center BH). They found that a fraction of galaxies can host both a central BH and an off-center BH, thus forming a binary system. When both BHs have L_bol > 10⁴³ erg s⁻¹, they form an AGN pair.

We successfully traced the evolution of five such AGN pairs from Horizon-AGN simulation with primary quasar luminosity > 10^45.3, and the secondary > 10^44.3 (M. Volonteri, 2021 private communication). The real space separations of these sources are between 5 and 30 kpc. These criteria are the same as those we used to select our observational candidates. This left us with five pairs for which we are able to trace the full merger sequence (Figure 12). Star symbols represent the primary BH (central BH), and the dots represent the secondary BHs (off-center BH); each pair is marked with a separate color. The BH mass and galaxy stellar mass are recorded in four stages in the merger scenario: (1) initial stage when both sources are labeled as AGN and are well-separated. (2) The host galaxies merge first. (3) The BHs form a close SMBH binary separated by 1 kpc or less, where they become unresolved at this stage, thus the M_BH is the total mass of the pair. (4) They coalesce into one single source. Figure 12 shows the evolution trajectories of these five AGN pairs. In some cases, the final coalescence stage is very close to the SMBH binary stage, thus only three stages are visible. The thick black line is a fit to the M_BH − M_* correlation based on the observational data set of local AGNs (from Häring & Rix 2004):

$$\begin{eqnarray}&&\mathrm{log}\left({M}_{\mathrm{BH}}/{10}^{7}{M}_{\odot }\right)=0.27+0.98\mathrm{log}\left({M}_{* }/{10}^{10}{M}_{\odot }\right).\end{eqnarray} \tag{ 4 }$$

The 1σ confidence interval is indicated in gray.

Zoom In Zoom Out Reset image size

To compare our data with the evolutionary traces in Horizon-AGN, we plot black hole mass and stellar mass based on the image decomposition of our dual quasars in Figure 13. In total, we are able to measure BH masses and stellar masses for five quasar pairs: J0847, J1214, J1416, J2209, and J2337. Here, our quasar pairs are plotted as colored points, and circles with the same color plot the two sources in a pair. Unlike the simulation, here we can only capture one phase in the merger process. As a highly simplified exercise, we sum the values of the two components in each pair as indicated by the star symbols, mimicking their possible final states after coalescence. Note that, in two cases, J0847 and J1416, their host galaxies cannot be separated; the model prefers to fit the image with one Sérsic profile, thus we have only one measurement for the stellar mass. The gray dots are 3315 single quasars at 0.2 < z < 1.0 selected from SDSS DR14 (see details in Li et al. 2021 for comparison). We found that our dual quasars mostly have overmassive BHs compared to the local correlation, and this excess persists after the merger. This is partly due to the bias of the SDSS survey that only the most luminous and massive quasars are detected.

Zoom In Zoom Out Reset image size

We found similar tendencies in the left panel of Figure 12 and in our results (Figure 13). Either before or after the merger, most of the dual AGNs in the simulation are located above the local correlation. This indicates that the BHs have already grown to a massive state before the merger. This agrees with the distribution of our data in Figure 9, showing the dual quasars possess BHs with masses comparable to single quasars. In fact, the mergers do not bring them any closer to the correlation; rather, the merger products lie even further above the correlation. we trace the evolution of the Horizon AGNs with redshift in Figure 12 (middle and right panels). We find that the mergers do lead to a boost in the BH mass (see where the dashed lines and solid lines merge), especially the blue case. However, the gain in stellar mass is not as dramatic. In one case (the green one), we even see a decrease in stellar mass during the merger. We suggest that, during the merger, the systems lost matter due to the interaction, stellar feedback caused by starburst, or AGN-driven effects, while the central black holes of course do not lose mass. This is why we see the mergers bring the AGN pairs above the local correlation.

However, biases may still be present in these results: the less massive (i.e., fainter) BHs could be missed by our selection algorithm, which prefers point-like sources. When the central source is not bright enough to outshine its host, we may fail in identifying such candidates. Also, with our spectroscopic observations, we may fail to identify a low-mass BH when the spectrum is too faint to recognize the emission lines. Therefore, the pairs presented in this work might not represent the bulk population of mergers.

On the other hand, as suggested by Jahnke & Macciò (2011), statistical convergence by mergers in hierarchical assembly could be an intrinsic mechanism for establishing the mass relations. In this case, no causal connection is required to build up the relation. An estimation of merger fractions of quasars would examine this hypothesis via indicating the statistical importance of mergers to BH growth.

Therefore, a larger observational sample is required to make more definitive statements about the relative evolution of black holes and their hosts in galaxy mergers with dual black hole growth. In a future study, we will expand our sample size both from candidate selection and additional follow-up spectroscopic observations. The HSC/SSP is continuing to release new data, thus we plan to update our parent samples from PDR2 to S20A (DR3). Together with a more refined selection algorithm, we expect to grow our parent sample by a factor of 2. We have also achieved new data from the subsequent Gemini/GMOS queue and the 3.58 m ESO New Technology Telescope (NTT). We will report our findings in follow-up works.

Differing from physically associated pairs, there is no inevitability for projected quasar pairs to form, no more than by chance alignment. However, they can still provide useful information. For example, Prochaska et al. (2013) studied the H i environment of the foreground quasars using the absorption features of background quasars based on a sample of 650 projected quasar pairs with proper separations between 30 kpc and 1 Mpc. Their follow-up study (Prochaska et al. 2014) shows that a similar method can also be applied to study the environment of C ii λ1334 and C iv λ1549. Since the separations of our candidates are smaller than 30 kpc, these candidates would be helpful to study the very central region of the quasars' environment, where an enhanced absorption is expected. In this work, we have discovered two new projected pair quasar systems at kpc-scale separation.

B.1.1. SDSS J021352.67−021129.4

The SDSS BOSS program reports this target as a broad-line quasar at z = 2.7794. Both of the sources are very blue, below the stellar locus, and in the quasar region (Figure 5). The host galaxies are not resolved at this high redshift, thus two PSFs are enough to model their optical emission (Figure 16a).

Zoom In Zoom Out Reset image size

The FOCAS spectrum covers C iv λ1549 and C iii] λ1909. Unfortunately, the observing condition was not ideal, with a seeing around 2^'', plus the target was very close to the moon. This resulted in a noisy profile and a cutoff at the blue end (Figure 16b). The redshift of source A is solid because Lyα, C iv λ1549, and C iii] λ1909 are detected by SDSS and our own spectrum. However, the solution of source B is ambiguous. We cannot confidently say that the strong emission line from the companion is C iv, because there is no obvious C iii] λ1909 feature at the corresponding position. We fit the line with a Gaussian function (Figure 16c), which resulted in a velocity offset from component A of 5779 ± 197 km s⁻¹. If we consider the line to be C iv λ1549, then it is shifted by 5932 ± 119 km s⁻¹ relative to source A, hence we classify it as a projected quasar pair here. On the other hand, if the emission line is Mg ii λ2798 or C iii] λ1909, the redshift of source B would be 1.3021 and 2.3742, respectively. The C iii] assumption is less likely to be true, because in that case, C iv λ1549 is expected to appear at the blue end of the spectrum, which is actually absent.

B.1.2. SDSS J225147.82 + 001640.5

SDSS reports the northern source as a broad-line quasar at z = 0.4096. According to the HSC image, the main quasar is hosted by a face-on barred spiral galaxy that is well-resolved after we subtract the point sources (Figure 17(a) data—point source).

Zoom In Zoom Out Reset image size

The FOCAS spectrum of this target covers the Hα and Hβ regions of the main quasar. However, we did not find similar features at similar positions in source B. Instead, they appear to be located at different wavelengths (Figure 17(b)). We find the Hβ line is located near the 7600 Å sky absorption feature, together with the [O iii] doublet, which indicates a redshift of 0.5778. We tried to reconstruct the line profile that is affected by the sky absorption through line fitting in Figure 17(c). The bottom left figure shows a trial with a broad Hβ component that matches the redward part of the Hβ profile well. Then we fit the profile with narrow components only, with the separation being fixed. We can see that, with only narrow Hβ, the line profile is not fit adequately. This indicates the necessity of the broad component in this source (component B). The width of the broad component is 3408 ± 203 km s⁻¹. Therefore, we classify this target as a projected quasar pair. The Hα line of the companion should be located at around 10,250 Å, which is likely buried in the noise at the red end of our spectrum.

In three cases, we found the companion to be a galaxy, without signs of a quasar, at the same redshift as the primary quasar. All of these cases are located at low redshift, and the HSC images are able to resolve the features of both galaxies. These systems show signs of interaction, indicating an ongoing merging event.

B.2.1. SDSS J021930.51−055643.0

SDSS classified this target as a broad-line quasar at z = 0.2917. The fiber is centered on the northern source as shown in the left panel of Figure 18(a). We decompose the emission based on model fitting using two PSFs and two Sérsic functions (Figure 18(a)). The photometric magnitude of the companion is fainter by 1.2 mags on average in all five bands. We see that the host galaxies of these two sources are irregular, given the strong residuals, thus indicating an ongoing interaction between them.

Zoom In Zoom Out Reset image size

We provide a zoom-in view of the FOCAS spectrum of the Hα and Hβ regions in Figure 18(b). We found that the Hα emission of the main quasar is spatially extended and blends with emission from the companion. To disentangle the emission from each component, we perform an extraction of the profiles using astropy and scipy packages in Python. We fit the spatial profile with two Gaussian curves at each wavelength position. The central positions of the Gaussians are constrained strictly to each trace (±2 pixels from the peak) to deblend the signals. The aperture size is set to be a constant value (1.57 pixel in Gaussian width (σ) for this case) for both sources, which is decided before the extraction via a prior fitting to the binned data. Here, we show the performance of this strategy at the position of Hα as an example (Figure 18(c); spatial line profiles). The x-axis is the spatial axis along the vertical direction in the 2D spectrum in Figure 18(b), and the y-axis shows the counts at each pixel. Both sources are fairly well-described with a Gaussian curve, we take the area under each curve as the photon counts at that wavelength position. Besides Hα, we see other emission lines in the main quasar spectrum also tend to be more extended on the side of the companion than the other side. We suggest this to be evidence of an interaction between the two sources.

As a result of line fitting, we did not find significant broad components in both Balmer lines of the companion. We plot the line ratios of the pair in the BPT diagnostic diagram (Baldwin et al. 1981; Kewley et al. 2006, Figure 19), where the green star is the main quasar and the green dot is the companion. We have a very highly constrained measure of the lines from our spatial profile decomposition strategy; the error is smaller than the mark size. We find the companion is located in the AGN/Seyfert region of the diagrams, but considering the close separation (7.74 kpc at this redshift) and the blended 2D spectra, it is possible that the emission lines in the companion are ionized by the main quasar, not another accreting SMBH. Here, we classify it as a quasar–galaxy pair, although it could be an obscured quasar. Further confirmation is required with X-ray observations.

Zoom In Zoom Out Reset image size

B.2.2. SDSS J022105.64−044101.5

SDSS reports the left source (Figure 20(a)) as a broad-line quasar at z = 0.1986. According to the HSC image, the sources have very similar colors, located at the lower right side of the stellar locus (Figure 4), a relatively sparse region of quasars.

Zoom In Zoom Out Reset image size

Hα is observed in our GMOS spectrum. The spectra of each component are well-separated (Figure 20(b)). We see moderate Hα emission in the spectrum of source B, with no broad component required in the fitting (Figure 20(c), left bottom panels). Different from J0219-0556, we do not have Hβ+[O iii] information for J0221-0441, so we can only plot one axis on the BPT diagram (Figure 20(c), right upper panels). This object is located at the AGN and LINER region, due to its strong [N ii] emission. Therefore, our first conclusion was a LINER. As mentioned by Cid Fernandes et al. (2011), the LINER region of BPT diagrams consists of two distinct classes: weak active galactic nucleus (wAGNs) and retired galaxies (RGs). The latter are the result of ionization photons supplied by hot evolved low-mass stars (HOLMES) instead of AGN. One purpose of adding an extra WHAN diagram (Cid Fernandes et al. 2011) is to distinguish these two classes. Source B has a very high [N ii]/Hα ratio, up to 1.7, which places it at the right side of the diagram. The continuum level of this source is around 7.1 × 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹, compared to a flux of 20.3 × 10⁻¹⁷ erg s⁻¹ cm⁻² of Hα. This gives an equivalent width of 2.8 Å, very close to the 3 Å criterion suggested by Cid Fernandes et al. (2011) to distinguish wAGN and RGs. Our data point almost falls on the borderline of these two classes (Figure 20(c) middle bottom figure). Of course, this diagnostic is not a dichotomy, but rather is a continuous evolution from one to another. At the W_Hα borderline, the AGN contributes between ∼1/3 and 2/3 of the ionizing power and decreases as W_Hα decreases. In our case, this fraction is not negligible. We argue that there is an AGN contribution to the lines at some level, but according to the diagnostic diagram, we classify this as a galaxy.

According to our model in Figure 20(a), these two sources share one host galaxy and they are very close to each other, separated by only 123, which projects to 4.04 kpc at redshift 0.199. Therefore, we suggest this to be a late-stage merging system with one of the SMBHs not strongly activated yet or in a post-AGN phase.

B.2.3. SDSS J231152.90–001335.0

SDSS reports the northern source (Figure 21(a)) to be a broad-line quasar at z = 0.3473. The fainter red companion resides to the southeast of the main source. The continuum level of the main source is about seven times that of the companion (Figure 21(b)).

Zoom In Zoom Out Reset image size

The redshift of the companion is supported by both absorption lines and a weak Hα line. We matched the Hβ, Mg, and Na absorption features in the spectrum to a passive galaxy template shifted to z = 0.345, i.e., the same redshift as the main source (Figure 21(c), upper panel). At the corresponding position, we find a bump that is expected to be Hα and fit it with a broad component (Figure 21(c), lower panels). While we matched the peak position of the two sources, the other components in the companion's spectrum are not convincing. It is hard to say whether a broad component is really required. If we believe there is a bump at Hα, then the velocity offset between the narrow components is 366 ± 58 km s⁻¹. Our classification here is mainly based on the absorption lines, which leads to a result of a passive galaxy.

However, according to HSC imaging, the host galaxy of the main quasar overlaps with the companion source, so it is unavoidable that our spectrum of source B suffers from contamination. This issue cannot be well-resolved using the spatial profile fitting strategy (Figure 21(c), upper right panel) due to the low S/N ratio of source B. One result of this contamination is that all the features, including the absorption lines and Hα emission we see in the spectrum of source B, actually come from the host galaxy of the main quasar. In that case, our classification will no longer hold. Considering the morphology of the host galaxy of the main quasar, it is difficult for our decomposition model to reconstruct its features, as shown by the normalized residual image. We suggest that it is likely an irregular galaxy. The separation of these two sources is 175, which projects to 8.6 kpc at redshift 0.346.