A gravitationally lensed quasar discovered in OGLE

Abstract

We report the discovery of a new gravitationally lensed quasar (double) from the Optical Gravitational Lensing Experiment (OGLE) identified inside the ∼670deg² area encompassing the Magellanic Clouds. The source was selected as one of ∼60 ‘red W1 − W2’ mid-infrared objects from WISE and having a significant amount of variability in OGLE for both two (or more) nearby sources. This is the first detection of a gravitational lens, where the discovery is made ‘the other way around’, meaning we first measured the time delay between the two lensed quasar images of −132 < t_AB < −76  d (90 per cent CL), with the median t_AB ≈ −102  d (in the observer frame), and where the fainter image B lags image A. The system consists of the two quasar images separated by 1.5 arcsec on the sky, with I ≈ 20.0 mag and I ≈ 19.6 mag, respectively, and a lensing galaxy that becomes detectable as I ≈ 21.5 mag source, 1.0 arcsec from image A, after subtracting the two lensed images. Both quasar images show clear AGN broad emission lines at z = 2.16 in the New Technology Telescope spectra. The spectral energy distribution (SED) fitting with the fixed source redshift provided the estimate of the lensing galaxy redshift of z ≈ 0.9 ± 0.2 (90 per cent CL), while its type is more likely to be elliptical (the SED-inferred and lens-model stellar mass is more likely present in ellipticals) than spiral (preferred redshift by the lens model).

gravitational lensing: strong, methods: observational, quasars: general

1 INTRODUCTION

Gravitational lensing is a powerful tool; acting as a natural telescope it can potentially allow astronomical observations to reach smaller and fainter objects than would otherwise be observable (Schneider, Ehlers & Falco 1992; Schneider 2006). In strong lensing, a massive foreground galaxy or cluster deflects the light from a background object, resulting in multiple images of the distant source, that is typically magnified by a factor of 10. This provides a wide range of applications of the lensing phenomenon. Lensing analyses support fundamental open questions in astrophysics and cosmology, and have been used to investigate the dark matter density profiles of lensing galaxies (e.g. Schechter & Wambsganss 2002; Pooley et al. 2009; Vegetti et al. 2014; Sonnenfeld et al. 2015), substructures (e.g. Dalal & Kochanek 2002; Kochanek & Dalal 2004; Vegetti et al. 2010), the evolution of galaxies (e.g. Koopmans et al. 2006; Kostrzewa-Rutkowska et al. 2014; Oldham et al. 2017), and the properties of dark energy (e.g. Link & Pierce 1998; Zhu 2000; Suyu et al. 2010; Collett et al. 2012), just to mention a few.

Strong gravitationally lensed quasars with measured time delays between the multiple images provide a straightforward and competitive approach to determining the Hubble constant, entirely independent of the local distance ladder, and allow us to study the expansion history of the Universe. The time delay method was proposed by Refsdal (1964), 15 yr before the discovery of the first strong lensing system, Q0957 + 561, in 1979 (Walsh, Carswell & Weymann 1979). Time delays in gravitationally lensed quasar light curves of different images depend on the gravitational potential of the lensing galaxy and the cosmological distances between the observer, the lens, and the source. In the past few years, there have been many attempts to measure the Hubble constant and other cosmological parameters with increasing precision and accuracy from time delays (Suyu et al. 2014, 2017). The time delay studies require not only detailed models for the mass profile of the main lens but also taking into account any matter structures along the line of sight (Rusu et al. 2017). Also recently, long time monitored systems were used to determine the Hubble constant with a precision below 5 per cent (Bonvin et al. 2017).

Furthermore, in multiple quasar images, we may observe the microlensing effect caused by the stars of the lens (e.g. Wambsganss, Paczynski & Schneider 1990; Wambsganss & Paczynski 1991; Schechter & Wambsganss 2002). This provides stringent constraints on the morphologies of the neighbourhood of black holes in lensed quasars, for example, kinematics and geometric structure in quasar broad emission-line regions (Guerras et al. 2013; Braibant et al. 2014, 2016). Detailed studies of the accretion discs showed that their sizes are robustly larger than expectations from the thin disc theory and the temperature profile seems to be steeper than expected (Bate, Webster & Wyithe 2007; Bate et al. 2008; MacLeod et al. 2015; Jiménez-Vicente et al. 2014, 2015b); however, opposite results have also been reported (Poindexter, Morgan & Kochanek 2008). Because the microlensing effect is also sensitive to the fraction of mass enclosed in stars to dark matter halo mass close to lensed images, at the same time we obtain the stellar composition of the lens galaxies (Bate et al. 2011; Schechter et al. 2014; Jiménez-Vicente et al. 2015a). While microlensing of lensed quasar images has numerous applications to study the AGN physics, on the other hand, it affects the light curves of quasar images, producing problems in the accurate time delay estimation, and, in turn, introducing biases in the characterization of lens systems (Tie & Kochanek 2018). However, recently several algorithms were applied to take into account microlensing in the time delay measurements (e.g. Bonvin et al. 2017).

There are less than two hundred strongly lensed quasars known to date, because (1) they are intrinsically very rare objects and (2) they are discovered by various techniques, each having a number of biases. Considering that we search for lensed sources blended with the light of the lens galaxy, the problem of the deflector's flux contribution impedes the lensed quasar hunt. Systematic searches of gravitationally lensed quasars have been conducted within several surveys, for example, by the Sloan Digital Sky Survey (SDSS) and the Dark Energy Survey (DES). Recently, several techniques were explored to search for lensed quasar systems, for example, identification of spectroscopically confirmed quasars and searching for nearby objects of similar colours (large separation systems) or spectroscopically confirmed quasars with extended morphology (small separation systems) (Oguri et al. 2006; Inada et al. 2008; More et al. 2016), selection of lensed quasars using an image subtraction method as they are time variable sources (Kochanek et al. 2006a; Lacki et al. 2009), data mining and mashine learning techniques with a combination of colour cuts and model-based selection (Agnello et al. 2015; Schechter et al. 2017), Gaussian mixture models of colours (Ostrovski et al. 2017), identifying candidate systems of red galaxies with multiple neighbouring blue objects (Lin et al. 2017), and variability studies of quasars (Berghea et al. 2017).

One of the primary goals of the Optical Gravitational Lensing Experiment (OGLE) – a long-term sky-variability survey (Udalski, Szymański & Szymański 2015) – has been monitoring of selected gravitational lenses. Since 1997, OGLE has performed a continuous, long-term monitoring of the gravitational lens QSO 2237+0305 (Udalski et al. 2006), where the resulting homogeneous OGLE data set is the most extensive photometric coverage of this source, spanning now two decades and showing numerous unique microlensing effects in this extraordinary object. For the second monitored system, HE 1104–1805, the time delay of 157 ± 21  d between the light curves of the two images was determined using the photometric data collected between years 1997 and 2002 (Wyrzykowski et al. 2003).

In this paper, we present the discovery of the first lensed quasar found in the OGLE data set covering the Magellanic System, the area that has been regularly monitored since 2010. Inside the observing area of about 670deg² around the Magellanic Clouds to the photometric depth of I < 21 mag, statistically there should exist about 10 lensed QSO systems, both doubles and quads (Oguri & Marshall 2010). Finding all of them in a search is difficult as the detection efficiency is rarely high, mainly due to the fact that quasar variability amplitude becomes buried in the photometric noise for objects fainter than I ≈ 20 mag. This is true in our search, where we combined the OGLE optical data with the mid-IR (mid-infrared) Wide-field Infrared Survey Explorer (WISE) data, and searched for the intrinsic optical quasar variability (in particular time lags) in OGLE counterparts to ‘red’ mid-IR sources, without prior spectroscopic analysis. The details of our search and results will be presented in subsequent sections.

This paper is organized as follows. In Section 2, we introduce our data sample and present selection method, then describe the lensed quasar candidate (that we subsequently confirm) in Section 3. In Section 4, we discuss results of modelling and lens and source properties. The paper is concluded in Section 5. Throughout this paper we assume a flat Λ cold dark matter (ΛCDM) concordance cosmological model of the Universe with parameters |$\Omega _\Lambda = 0.7$|⁠, Ω_M = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹, h = 0.70.

2 THE OGLE SURVEY AND LENSED QUASARS DETECTION STRATEGY

2.1 The OGLE variability survey

The key data used in this work were collected during the fourth phase of the OGLE project (hereafter OGLE-IV), between the years 2010 and 2017, with the 1.3-m Warsaw telescope located at Las Campanas Observatory, Chile, operated by the Carnegie Institution for Science. The telescope is equipped with a 32 chip mosaic CCD camera, with a total field of view of 1.4deg² (Udalski et al. 2015). The entire area of the Magellanic Clouds region is covered by 475 OGLE-IV fields that amounts to about 670deg². Only two filters I and V band have been used by OGLE-IV. About 90 per cent observations were carried out in the I band with magnitude range from 13 to 21.5 mag. The number of collected epochs varies between 100 and 700, and observations of a given field are repeated every 2–4 d. The photometry was obtained using the difference image analysis (DIA) method, as implemented by Woźniak (2000).

2.2 Lensed quasar candidates

Since 2010, The Magellanic Quasars Survey (MQS) has increased the number of quasars known behind the Magellanic Clouds by almost an order of magnitude using the OGLE-III data (Kozłowski et al. 2013). It was then a natural consequence to perform a similar extensive search of quasar candidates in the increased OGLE-IV area of the Magellanic System (Kozłowski et al. in preparation).

The primary sample of quasar candidates was obtained by locating all objects in the WISE survey (Wright et al. 2010; Cutri & et al. 2014) fulfilling any of the mid-IR colours criteria for quasars proposed by Stern et al. (2012) and Assef et al. (2013), then cross-matching them with objects in the OGLE data base, and finally performing a variability analysis of the OGLE objects, as a result isolating a sample of quasar candidates (∼15 000 objects).

We used this sample of quasar candidates in further identification of lensed quasar candidates. First, we investigated their environments within a radius of r < 6 arcsec (approximately W1 WISE FWHM) to identify neighbouring objects in the OGLE data base. Of the 15 000 sources, 63 objects turned out to be ‘our early’ lens candidates based on the similar rms. in the light curves. Next, for each such a lens candidate (now two or more neighbouring objects), we visually compared both their colour and variability properties. The highest probability candidates were required to have similar variability patterns in both images (and so plausibly the time delay) as reported by the Javelin¹ software (Zu, Kochanek & Peterson 2011; Zu et al. 2013, 2016). We also used a dedicated DIA (to prevent ‘ghost’ variability in the neighbour caused by the true variable object – we created two images, where we stacked a number of images when the quasar is bright and when it is faint, and then subtracted them to check if the residua due to flux change appear on both locations). Three most promising candidates were observed spectroscopically with the New Technology Telescope (NTT), of which only one turned out to be a genuine lensing system (the other two are unlensed quasar+star pairs, presented in Fig. A1). The false positive rate might be lower by, first, a visual inspection of both the light curves and difference images (ghosts), and then by the spectral energy distribution (SED) fitting. This method allowed for a straightforward identification of the lensed quasar candidate system without extensive spectroscopic observations.

3 THE LENS SYSTEM

Fulfilling the requirements from Section 2, the candidate considered in this paper had to have ‘red’ WISE colours (W1 − W2 > 0.8 mag from Stern et al. 2012 or a combination of W1 − W2 and W2 data from Assef et al. 2013) but also two optical counterparts well separated on the sky, both showing similar photometric variability, albeit allowed to be shifted in time. In Fig. 1, we present the postage-stamp cutouts for the candidate images A and B, at (RA, Dec.) =(02:18:12.34, −73:35:35.8) and (02:18:12.47, −73:35:37.2), respectively.

Figure 1.

The OGLE lensed quasar as seen in OGLE (top panels) and GROND (bottom panels) is shown. The two lensed quasar images are marked in the top left-hand panel, where A (B) is the brighter (fainter) image. They are separated by 1.5 arcsec and are readily resolvable in the best-seeing OGLE frames of 0.7 arcsec. The top right-hand panel shows the colour composition of I- and V-band frames. The bottom row presents g″r″i″z″JHK GROND observations (from the left- to right-hand panel). The lens galaxy becomes visible as a blend to image B in the infrared data. All panels cover approximately 8 × 8 arcsec². North is up and east is to the left-hand side.

3.1 Spectroscopic confirmation

We obtained a low resolution spectrum of our candidate with the ESO Faint Object Spectrograph and Camera (EFOSC2) on the 3.58-m Ritchey-Chrétien NTT on 2016 September 27. The spectrum of 300 s was taken with a 1.2 arcsec slit, aligned to contain both objects, and using the 13 grism, which gives a wavelength range of 3685 to 9315 Å. The data were binned by factors of 2 along the slit and in the spectral direction. This setup yields a dispersion of 5.54 Å per pixel and a full width at half-maximum (FWHM) of 21.24 Å. The seeing ranged between 0.5–0.6 arcsec. The pair of images in the candidate system was completely resolved. The 2D spectrum was reduced by fitting two Gaussians and extracting 1D spectrum. The confirmation spectrum is shown in Fig. 2. Several common AGN broad emission lines (Lyα, Si iv, C iii], and C iv) at the redshift of z_S = 2.16 are present. We were unable to recognize any features from a potential lensing galaxy, simply because it was expected to be too faint in this wavelength range (see SED in Fig. 3).

Figure 2.

NTT spectra of the two images (A in red, B in blue) of the lensed quasar. For comparison, in black, we also show a composite AGN spectrum from Vanden Berk et al. (2001). The source redshift is z_S = 2.16. Four common quasar emission lines are marked with the vertical dotted lines and labelled at the top.

Figure 3.

SED of the OGLE lens system (top panel) and the χ² distribution for the lens redshift (bottom panel). Top panel: the black, red, and blue points come from the OGLE, GROND, Spitzer, and WISE surveys (Table 1), where black is for the two images of the quasar being unresolved, while the red and blue points are for the resolved images A and B, respectively. Image A (red) is assumed to be a sum of the AGN and its host (at z_S = 2.16), while image B (blue) due to proximity of the lens galaxy is a sum of the AGN, its host, and the lensing galaxy. In both sub-panels, the solid black line is the best-fitting model, being the sum of the AGN and the host, and a lens galaxy (dashed line): spiral or elliptical. The magenta point is the constraint on the lens magnitude in I band from the PSF fitting. Bottom panel: Δχ² as a function of the lens redshift for two types of lensing galaxies and two types of host galaxies.

3.2 Estimation of the lens redshift

The estimate of the lens redshift comes from the SED fitting. In order to build an SED, we combined a number of data sets from the literature and observations, totalling 20 bands from ultraviolet (UV) to mid-IR. The near-UV (NUV) data came from Galex (Bianchi et al. 2011), while both V- and I-band magnitudes were obtained by OGLE-IV (Udalski et al. 2015). We performed dedicated observations of the candidate with GROND (PI: Jochen Greiner) in g″, r″, i″, z″, J, H, K filters (Greiner et al. 2008), operated as an imaging instrument at the MPI/ESO 2.2m telescope in La Silla, Chile. The mid-IR data were obtained both from the Spitzer IRAC and MIPS detectors (C1–C4 and C5 channels; the SAGE project; Meixner et al. 2006; Gordon et al. 2011) and WISE (W1–W4 bands; Wright et al. 2010; Cutri & et al. 2014) space-based telescopes, both covering the 3–24 micron range (Table 1). These broad-band magnitudes were converted to luminosities, as prescribed by Kozłowski (2015), using the standard ΛCDM cosmological model.

Table 1.

Combined observational properties of the system across the electromagnetic spectrum (SED).

Filterab	λ	Image	Magnitudec	∝ νFν
NUV (AB)	0.23	AB	21.66 ± 0.30	0.54 ± 0.13
u (AB)	0.35	AB	19.57 ± 0.11	2.41 ± 0.23
g’ (AB)⋆	0.48	A	19.77 ± 0.02	1.49 ± 0.03
V (Vega)	0.55	A	19.94 ± 0.02	1.17 ± 0.01
r’ (AB)	0.62	A	20.22 ± 0.02	0.76 ± 0.01
i’ (AB)	0.76	A	20.24 ± 0.02	0.62 ± 0.01
I (Vega)	0.79	A	19.52 ± 0.02	0.75 ± 0.01
z’ (AB)	0.91	A	20.11 ± 0.04	0.57 ± 0.01
J (Vega)	1.24	A	18.70 ± 0.11	0.68 ± 0.04
H (Vega)	1.66	A	18.68 ± 0.26	0.33 ± 0.10
K (Vega)	2.16	A	17.65 ± 0.27	0.42 ± 0.11
g’ (AB)	0.48	B	20.28 ± 0.02	0.93 ± 0.02
V (Vega)	0.55	B	20.35 ± 0.02	0.75 ± 0.04
r’ (AB)	0.62	B	20.58 ± 0.02	0.53 ± 0.01
i’ (AB)	0.76	B	20.45 ± 0.03	0.50 ± 0.01
I (Vega)	0.79	B	19.93 ± 0.02	0.48 ± 0.01
z’ (AB)	0.91	B	20.25 ± 0.03	0.50 ± 0.02
J (Vega)	1.24	B	19.00 ± 0.12	0.50 ± 0.05
H (Vega)	1.66	B	18.35 ± 0.15	0.44 ± 0.06
K (Vega)	2.16	B	17.37 ± 0.20	0.55 ± 0.08
W1 (Vega)	3.35	AB	15.42 ± 0.04	0.98 ± 0.03
C1 (Vega)	3.56	AB	15.50 ± 0.09	0.79 ± 0.06
C2 (Vega)	4.51	AB	14.70 ± 0.08	0.83 ± 0.06
W2 (Vega)	4.60	AB	14.69 ± 0.04	0.77 ± 0.03
C3 (Vega)	5.69	AB	13.88 ± 0.12	0.89 ± 0.61
C4 (Vega)	7.98	AB	13.24 ± 0.12	0.64 ± 0.07
W3 (Vega)	11.56	AB	12.39 ± 0.18	0.49 ± 0.08
W4 (Vega)	22.09	AB	<9.41	<1.02
C5 (Vega)	23.68	AB	9.72 ± 0.08	0.62 ± 0.16

^aC/W stands for Spitzer-IRAC(MIPS)/WISE bands.

^bThe photometry is provided in either AB or Vega systems.

^cThe magnitudes are extinction corrected.

*SDSS magnitude system as described in Fukugita et al. (1996).

To fit the SED (Fig. 3 and Table 1), we used the low-resolution templates of AGN and galaxies from Assef et al. (2010). We excluded the Galex NUV point from the fitting, because of its unreasonable low value as compared to models. The origin of this discrepancy lies in the fact that the UV flux in distant AGNs is absorbed by neutral hydrogen in the intergalactic medium, a feature that is not taken into account in both the templates used and our simplistic SED modelling. Both the OGLE and GROND data have sufficiently high resolution to measure fluxes for both lensed images, while the Spitzer and WISE data show only a single object due to blending or insufficient resolution. Whenever possible, we fitted the image A (as a sum of the AGN and the host galaxy: either spiral or elliptical) separately from the image B (a sum of the AGN, the host and the lensing galaxy that blends with the image B). We only roughly know the I-band magnitude of the lens galaxy (I ≈ 21.5 mag), that we used as an additional constraint in the SED fitting. We also know the source redshift, so we fixed it in the SED fitting to z_S = 2.16. We obtained the lens redshift of |$z_{\rm L}=0.94_{-0.28}^{+0.24}$| (90 per cent CL) for the spiral (insensitive to the AGN host type) and |$z_{\rm L}=0.84_{-0.21}^{+0.18}$| (90 per cent CL) for the elliptical galaxy (also insensitive to the host type). Best-fitting galaxy and AGN models as well as the χ² distributions are presented in Fig. 3. We have included the Galactic extinction model as proposed by Cardelli, Clayton & Mathis (1989) with R_V = 3.1, however, the results turned out to be non-physical and pointing to negative extinctions. Fixing the extinction to positive values provides a slightly lower redshift for the lens galaxy at a price of higher χ². From the SED fit, we can also estimate the lens galaxy stellar mass: (2.5 ± 0.5) × 10¹¹ M_⊙ for the elliptical and (3.3 ± 0.7) × 10¹¹ M_⊙ for the spiral, however, the latter is more than an order of magnitude too high for this type of galaxies (Schawinski et al. 2014).

3.3 Astrometry from OGLE

Basic properties of the lensing system are given in Table 2 and the finding chart is shown in Figs 1 and 4. The lensed system is clearly resolved as two point sources separated by 1.5 arcsec in the OGLE I-band imaging data. The image pixels were modelled to obtain the astrometry of the system. First, we determined the point spread function (PSF) by fitting a nearby star with a parametric Moffat profile. We used four different nearby stars (within ∼ 1 arcmin) and compared the residuals and the Moffat shape parameters. All the stars are brighter than the quasar images. The brightest star did not give a reasonably good fit consistent with shape parameters for the other stars and it was excluded from the final fit. The Moffat β parameter changes between fits but is ∼2.5. For the final lensed quasar fit, we accounted for the uncertainty in the PSF, hence we allowed the PSF parameters to vary within the measured values (of the other 3 stars) and then marginalize over them. The system was then modelled as a combination of two quasar images (two PSFs) and a galaxy profile, that was set to be a de Vaucouleurs profile convolved with the PSF profile. The free parameters of our model were the 2D positions of the quasar images and lensing galaxy, and the galaxy shape parameters (flattening and effective radius). The model parameters were determined, along with their uncertainties, via sampling with the Markov chain Monte Carlo (MCMC) ensemble sampler – emcee (Foreman-Mackey et al. 2013).

Figure 4.

The results of image modelling. The OGLE I-band image is shown (top left-hand panel). It was modelled as a sum of two PSFs and a galaxy. The two lensed quasar images obtained from this model were removed in the top right-hand panel, leaving only the galaxy (marked with the red crosshair in all panels). Residuals after subtracting the best-fitting model including two PSFs is shown only in the bottom left-hand panel. It is clear such a model was insufficient to correctly describe the data. The best model included both the two PSFs and the galaxy, then the residuals became smooth and showed no obvious additional light (bottom right-hand panel). Astrometry and basic parameters of this system are provided in Table 2. All panels cover 8 × 8 arcsec², north is up, east is to the left-hand side.

Table 2.

Basic properties of the OGLE lens system.

Object	RA	Dec.	ΔRA	Δ Dec.	V	I	Redshift
			(arcsec)	(arcsec)	(mag)	(mag)
Image A	02:18:12.34	−73:35:35.8	0	0	20.09	19.61	2.16
Image B	02:18:12.47	−73:35:37.2	0.54 ± 0.01	−1.44 ± 0.01	20.50	20.02	2.16
Galaxy	02:18:12.42	−73:35:36.8	0.32 ± 0.03	−1.01 ± 0.03	–	≈21.5	≈0.9

In Fig. 4, we show the original data (top left-hand panel). They were modelled as a sum of two PSFs and a galaxy. The top right-hand panel presents the lensing galaxy with I ≈ 21.5 mag, that becomes visible after subtracting the two model PSFs (and not subtracting the galaxy) from the image. Modelling of this image with two PSFs results only in bad residuals after subtracting the model from the data (bottom left-hand panel). The data are sufficiently modelled only if both the images and the galaxy are included. The residuals after subtracting such a model from the data are clean (bottom right-hand panel). The lens light profile was modelled with de Vaucouleurs profile and we obtained the shape parameters: the effective radius of 0.55 ± 0.08 arcsec, axis ratio 0.51 ± 0.07 and position angle 101 ± 4 (in degrees measured east of north).

3.4 The OGLE photometry

In Table 3, we present the I-band light curves (419 epochs) for the two lensed images of the quasar, along with a basic observational log: airmass, seeing, and background, enabling trimming the data from bad weather points. Note, that the light curves are only best-calibrated to the Vega system, but not standard, as they lack the (V − I) colour term, because there is much less data points in V band than in I band. The image subtraction method typically reports underestimated error-bars that we corrected using the following equation |$\Delta I_{\rm new} = \sqrt{\left(\eta \Delta I_{\rm DIA}\right)^2+\epsilon ^2}$|⁠, where η = 1.558, and ε = 0.004 (Skowron et al. 2016).

Table 3.

The OGLE light curves for the two images of the lensed quasar. (For the guidance purposes, we provided below the first five rows of the full light curves, that are available in its entirety from the electronic version of the journal).

HJD-2450000	Image A	Image A	Image B	Image B	Airmass	Seeing	Background
(days)	I (mag)	σI (mag)	I (mag)	σI (mag)		( arcsec)	(counts)
5377.90992	19.801	0.156	19.708	0.143	1.528	1.55	696
5386.92763	19.898	0.080	19.870	0.080	1.450	1.12	253
5388.92893	19.846	0.087	19.942	0.097	1.441	1.24	245
5389.93796	20.010	0.125	19.974	0.122	1.428	1.61	336
5392.92310	19.975	0.154	19.875	0.140	1.435	2.01	320

From the OGLE photometry of template images (now in the standard Johnson–Cousins Vega system), we find that the colours of the lensed images are (V − I)_A = 0.48 mag and (V − I)_B = 0.78 mag, so the image B is redder than A, what is also readily visible in the top right-hand panel of Fig. 1. This may point to two scenarios: (1) image B is a composite of the lensing galaxy and the lensed quasar, or (2) image B is redder due to extinction in the lensed galaxy. If scenario (1) is the case, then we can easily estimate the I-band magnitudes of the two sources (assuming the V-band lens galaxy flux is negligible). The lensing galaxy brightness is I = 21.23 mag and the lensed image has the brightness of I_B = 20.02 mag (Table 2). The lensing galaxy I-band brightness is confirmed by the measurement obtained from the PSF fitting. On the other hand, from the SED modelling, we find the ‘pure’ flux ratio between images A (61 per cent of light of the total AB light) and B (39 per cent) is ∼1.6 with the formal uncertainty of 0.1, therefore the magnitude difference between images is ∼0.5 mag, however, the lensing galaxy contributes significantly to the flux of image B due to its proximity. In the near-IR, the flux ratio decreases to below 1.0, where the image B appears to be brighter than A, but in fact the best SED model of the lensing galaxy peaks around 2 μm and dominates the light of image B (bottom panel of Fig. 1 and top panel of Fig. 3).

3.5 The time delay

Time delays in gravitational lenses may be the key to understanding the lens system, including its mass distribution, but also are probes of cosmological parameters such as the Hubble constant (e.g. Blandford & Narayan 1986; Gorenstein, Shapiro & Falco 1988; Kochanek 2002).

Time delay(s) are estimated from light curves of the two (or more) images of the quasar. The necessary condition is that all light curves cover the same variability pattern intrinsic to the quasar. This usually means a need for about several years-long, well-sampled light curves having a good signal-to-noise (i.e. sufficiently bright images), because the typical AGN variability amplitude is small, a few tenths of a magnitude, on a time-scale of months to years (e.g. MacLeod et al. 2010; Kozłowski 2016). If the separation between the image(s) and the galaxy is small, or simply the galaxy is large (as in the case of QSO 2237+0305; e.g. Huchra et al. 1985; Udalski et al. 2006), then one or more images may undergo secondary microlensing by stars of the lensing galaxy. This happens rather frequently and is the main reason of troubles with correctly estimating time delays (e.g. Goicoechea et al. 1998; Tie & Kochanek 2018) and the image flux ratios (e.g. Wyithe & Turner 2002), while on the other hand enabling the measurement of the intrinsic properties of the quasar accretion disc (its temperature profile and size, e.g. Kochanek 2004; Dai et al. 2010; Morgan et al. 2012; Mosquera et al. 2013).

A straightforward method to measure time delay(s) is the cross-correlation (CC) of the two (or more) light curves, however because the data points are rarely regularly sampled, one needs to interpolate (or sometimes extrapolate) one of the light curves to match the points of the first light curve. The exact choice of inter/extrapolating/modelling method is tricky – one can simply use the spline/straight line or similar method of interpolation or one can try to model the quasar variability as the damped random walk (DRW; e.g. Kelly, Bechtold & Siemiginowska 2009; Kozłowski et al. 2010; MacLeod et al. 2010) or even more sophisticated continuous-time autoregressive moving average (CARMA; e.g. Kelly et al. 2014) models and then simultaneously search for the time delays.

In this paper, to increase the signal-to-noise ratio in the light curves, we combine the data into 10-d contiguous bins, and then whenever necessary interpolate such light curves using DRW model (Fig. 5). From a structure function analysis, we know that the asymptotic variability is SF_∞ = 0.30 mag and that the structure function power-law slope is γ ≈ 0.5, as required by DRW model. The DRW time-scale was fixed to 1 yr rest-frame (Kozłowski 2016). We cross-correlated such a DRW-interpolated light curve to obtain a weakly constrained time delay of t_AB ≈ −83  d and allowing for a wide range of negative values (Fig. 6), where the image A leads the image B what is in line with expectations from the geometry as image A is farther from the galaxy and it should always lead.

Figure 5.

Top panels: 10-d binned light curves of the two images A (red) and B (blue) of the lensed quasar. The upper panel shows the data as observed, while the bottom one is for image B shifted by −102 d (see Section 3.5). Bottom panels: We present the DRW models for both images, where the data are interpolated and extrapolated to within 10 d of each original data point. Similarly to the upper panels, we show the best model light curves as observed (top panel) and with the image B shifted by −102 d (bottom panel).

Figure 6.

The observed-frame time delay density probability obtained using Javelin on the cleaned (black) and 10-d binned (red) OGLE light curves is shown. The cleaned light curves have the 90 per cent time delay limits spanning −151 < t_AB < −75  d with the median of t_AB = −108  d, while the binned light curves have −132 < t_AB < −76  d and t_AB = −102  d, respectively. The blue line is the resulting time delay distribution from our python code, that returned the median time delay of t_AB = −90  d, where the 90 per cent confidence limits are −119 < t_AB < −57  d. In green (right y-axis), we show the CC of the 10-d binned data interpolated with DRW model, that broadly peaks at t_AB = −83  d.

Kochanek et al. (2006b) provide a method of estimating the time delay between the images of the quasar, where the intrinsic variations of the source can be approximated as a Legendre series of order N_src and are present in all images shifted by the time delay, while the secondary variations due to microlensing are approximated by another Legendre series of order N_μ. The orders of both series are typically set by using the F-test. In Table 4, we present the results for three low-order variability polynomials N_src = 10, 20, and 30, and the lowest order microlensing polynomial with N_μ = 1. The time delays are typically within the 1σ range of −100 < t_AB < −70  d.

Table 4.

Time delays obtained using the Kochanek et al. (2006b) method (see Section 3.5 for explanations).

Nsrc	Nμ	Parameter	Value
10	1	best	tAB = −89  d
10	1	1σ	−106 < tAB < −72  d
10	1	2σ	−124 < tAB < −56  d
10	1	3σ	−142 < tAB < −41  d
10	1	χ2/dof	222.6/219
20	1	best	tAB = −84  d
20	1	1σ	−97 < tAB < −72  d
20	1	2σ	−111 < tAB < −60  d
20	1	3σ	−126 < tAB < −48  d
20	1	χ2/dof	185.6/209
30	1	best	tAB = −84  d
30	1	1σ	−99 < tAB < −71  d
30	1	2σ	−116 < tAB < −49  d
30	1	3σ	−132 < tAB < −32  d
30	1	χ2/dof	167.1/199

The Javelin software is designed to search for time delays between two or more time series, where the time series is intrinsically modelled as the DRW model (Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010). However, it does not simultaneously model any additional variability in each light curves due to microlensing. We used Javelin to search for time delays for both full (but cleaned of bad points) and 10-d binned light curves. The time delay density probability distributions are presented in Fig. 6. The median time delays are t_AB = −108  d and t_AB = −102  d, respectively, while the 90 per cent confidence limits are −151 < t_AB < −75  d and −132 < t_AB < −76  d. As a cross-check, we also used our own python DRW-modelling time-delay searching code. The results are comparable, although the time delay is a little shorter, with the median t_AB = −90  d and the 90 per cent confidence limits of −119 < t_AB < −57  d.

The CC, the Kochanek method, and the python code seem to prefer slightly shorter, but not disjoint, time delays to the ones derived from the Javelin code. For further discussions, throughout this paper we will adopt the following median time delay of t_AB = −102  d and the 90 per cent confidence level of −132 < t_AB < −76  d, which are the results from the Javelin modelling of the 10-d binned data.

4 DISCUSSION

4.1 Modelling the lensing system

The available data and their uncertainties limit our ability to fully model the lensing system. As an approximation, we modelled the lens by adopting the singular isothermal ellipsoid (SIE) model for the lens mass density profile (Kormann, Schneider & Bartelmann 1994), and using the gravlens software (Keeton 2011). We obtained the Einstein radius of R_E = 0.78 ± 0.01 arcsec, axis ratio 0.91 ± 0.05, and position angle 90 ± 20 (in degrees measured east of north) from exploring the parameter space with MCMC sampling. The basic SIE model leads to the discrepancy in fitted values of axis ratio between mass and light profile but one has to remember that the external shear was not taken into account and/or the SIE model may be too simplistic. For the measured time delays, the gravlens modelling results point to the lens redshift (used here as a free parameter and so sampled for) and the lens mass. The lens redshift is found to be 0.83–1.12 (90 per cent CL) and the lens mass within the Einstein radius 4–6 × 10¹¹ M_⊙ (see Fig. 7). The estimated lens redshift from the SED modelling is in excellent agreement with our prediction from the lens model; however, the uncertainties are large.

Figure 7.

The time delay (left-hand panel y-axis) and the approximate lens mass (right-hand panel y-axis) as a function of redshift. The inclined line shows the prediction based on the simple SIE lens model. The horizontal blue line indicates the time delay of |t_AB| = 102  d, while the light blue region represents 90 per cent confidence levels. The two vertical lines mark the best-fitted SED redshifts for the spiral (left-hand panel; z_L = 0.94) and elliptical galaxy (right-hand panel; z_L = 0.84), when the source redshift is fixed to the spectroscopically measured value (z_S = 2.16) in the SED fitting. The green areas mark the 90 per cent confidence levels for these redshifts, while in red, we mark the overlapping 90 per cent confidence levels for both the time delay and redshifts. The best estimate for the lens galaxy mass is ∼5 × 10¹¹ M_⊙ within the Einstein radius.

The best-estimated total lens galaxy mass implies a stellar mass of 2 × 10¹¹ M_⊙ assuming a typical stellar-to-total mass ratio within the Einstein radius of ∼0.5 (e.g. Auger et al. 2010). This number can be compared against the galaxy mass distributions for late-type (spiral and irregular) and early-type (elliptical and lenticular) galaxies from SDSS. Our derived mass is by far more likely found in elliptical galaxies than in spirals (Schawinski et al. 2014). Another confirmation comes from the SED modelling, where we estimate the stellar mass for the lens being elliptical or spiral galaxy with (2.5 ± 0.5) × 10¹¹ M_⊙ and (3.3 ± 0.7) × 10¹¹ M_⊙, respectively. Both values broadly agree with the expectation from the lens model.

4.2 Detectability of lenses in OGLE

According to predictions presented by Oguri & Marshall (2010), there should be 1 gravitational lens for I < 19 mag, 3 lenses for I < 20 mag, 8 lenses for I < 21 mag, and 20 lenses (both doubles and quads) for I < 21.5 mag (the detection limit in OGLE) within the area of about 670deg² around the Magellanic Clouds. We adopted the typical quasar colour of (i − I) ≈ 0.4 mag (fig. 9 in Kozłowski et al. 2012) to convert the Oguri & Marshall (2010) i-band magnitudes into I band. Our candidate selection was partly based, however, on a requirement of detecting quasar variability. Kozłowski (2016) showed that the rms of the optical variability of quasars reaches about 0.25 mag, provided the light curves are sufficiently long. The rms of the photometric noise in OGLE-IV is 0.06 mag for I = 19 mag, 0.13 mag for I = 20 mag, and 0.30 mag for objects with I = 21 mag (Udalski et al. 2015). So the expected number of found lensing systems will be actually lower due to the fact that only a fraction of quasars is sufficiently variable around I ≈ 20 mag to exceed the photometric noise in OGLE. Ignoring other factors, such as the high stellar crowding, we should, in principle, detect fewer than five gravitational lenses. This number still could be lower in our case, because a fraction of the analysed area has been observed only for three years at the time of this search, so a fraction of quasars show lower rms variability. This implies that our detection efficiency is at least 20 per cent (>1/5).

5 CONCLUSIONS

We presented the results from our search for lensed quasars in the OGLE Variability Survey's data in the extended area of the Magellanic Clouds. Our first candidate, that we confirmed spectroscopically, is a lensed quasar at a redshift of z_S = 2.16, with two quasar images separated by 1.5 arcsec on the sky. After subtracting the two lensed images from a deep OGLE frame, we detected the lensing galaxy with I ≈ 21.5 mag. From the SED modelling spanning 0.2–24 microns and having the fixed source redshift, we found that the lensing galaxy could be both a spiral or an elliptical at the redshift of z_L ≈ 0.9 ± 0.2, although an elliptical galaxy is preferred. This redshift is also expected from the simplest SIE lens model with the measured time delay of −132 < t_AB < −76 d (90 per cent CL).

If the lens was a spiral this would raise a potential question on the amount of the expected extinction to be present in the SED of the quasar images (but undetected), because the lines of sight pass approximately 8 kpc (for image A) and 3 kpc (for B) from the galaxy lens centre. For a very massive spiral with the stellar mass of 3 × 10¹¹ M_⊙, this would likely mean passing through the spiral arms containing dust (unless it is edge on and has a favourable position angle, the latter being the case here as the PA is roughly perpendicular to the line joining the images). However, the lensing mass is nearly circular, which is inconsistent with an edge-on disc. Furthermore, our derived stellar mass is typically seen in elliptical galaxies and is generally too high, by roughly an order of magnitude, for spirals (Schawinski et al. 2014).

To obtain the time delay of −132 < t_AB < −76 d (90 per cent CL), we used several available methods that included the DRW light curve modelling with Javelin, but also with our own python code, the Kochanek et al. (2006b) method using the Legendre polynomials, and also CC. All these methods provided consistent results, although the Javelin software preferred slightly longer time delays.

We proposed another detection route of gravitationally lensed systems, first through the ‘red’ mid-IR WISE colours, typical for quasars, and then including the identification of point source pairs in optical data showing similar variability patterns. Our search method does not require any prior spectroscopic observations.

ACKNOWLEDGEMENTS

We are grateful to the anonymous referee for providing us with comments that improved this manuscript. ZKR acknowledges support from European Research Council Consolidator Grant 647208. SK acknowledges the financial support of the Polish National Science Centre through the OPUS grant number 2014/15/B/ST9/00093. TA and YA acknowledge support by proyecto FONDECYT 11130630 and by the Ministry for the Economy, Development, and Tourism's Programa Inicativa Científica Milenio through grant IC 12009, awarded to The Millennium Institute of Astrophysics (MAS). OGLE project has received funding from the Polish National Science Centre, grant MAESTRO 2014/14/A/ST9/00121 to AU. Part of the funding for GROND (both hardware as well as personnel) was generously granted from the Leibniz-Prize to Prof. G. Hasinger (DFG grant HA 1850/28-1).

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

Footnotes

1

https://bitbucket.org/nye17/javelin

SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

Table3_online.txt

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

APPENDIX A: OTHER OBSERVED CANDIDATES

Here, we present basic information on the two spectroscopically observed lensed quasar candidates, which turned out to be unlensed quasar+star pairs. In Fig. A1, we present the finding charts, the NTT spectra and the OGLE light curves in I band. Both candidates were observed on 2016 December 5 with the EFOSC2. Both quasars show common AGN broad emission lines at the redshift of z = 1.465 and 1.49. The spectra of nearby interlopers are also presented for comparison.

Figure A1.

Two spectroscopically observed candidates that turned out to be regular quasars with a nearby interloper. In the left-hand panel, we present 1 by 1 arcmin finding charts centred at these sources, in the middle panels, we show the NTT spectra with the main quasar emission lines marked as vertical lines, and in the right-hand panels, we present the OGLE-IV light curves. The best estimate Javelin time delays are −2 d (top) and 15 d (bottom).