Free Access
Issue
A&A
Volume 630, October 2019
Article Number A108
Number of page(s) 11
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201935427
Published online 30 September 2019

© ESO 2019

1. Introduction

The formation of super-massive black holes (SMBHs) at the centres of galaxies is still an unclear process. According to the hierarchical structure formation scenario, SMBHs can be created as a result of a major merger of two galaxies, each with its own nuclear black hole (Begelman et al. 1980; Volonteri 2010). Such systems can have an important effect on the build-up of the stellar halo through mechanical and radiative feedback when both black holes are undergoing an active phase. Also, the merger of such black holes may result in extreme gravitational wave events in the early Universe, which are the primary targets of the Laser Interferometer Space Antenna (LISA, e.g. Enoki et al. 2004). However, the lack of observed active galactic nuclei (AGN) pairs suggests that there must be a fast spiralling of the two black holes when they reach the final merging-stage at pc-scales (Mayer et al. 2007), and detecting such systems with 1 to 100 pc separation is extremely difficult, with only one pc-scale dual AGN system being confirmed to date (Rodriguez et al. 2006). However, the low detection rate of active binary SMBHs seems to be in agreement with the theoretical expectation of dual AGN if only one of the two SMBHs accretes and emits radiation during the merger. Then, in order to become a double AGN, the system must undergo at least two other major mergers (Volonteri et al. 2003). Under these assumptions, numerical simulations based on the optical and X-ray emission from AGN show that the fraction of dual AGN increases from 0.1 percent at z = 0 to only a few percent at z = 2, (Volonteri et al. 2016; Rosas-Guevara et al. 2019).

Observationally, the most common approach to identify such pairs of active SMBHs is to detect emission lines with an offset in velocity of a several hundred km s−1. This velocity offset can be seen as a double peak in the lines that originate in the narrow line regions of the two AGN, if they are spatially unresolved (e.g. [O III] lines, Liu et al. 2018). However, it is known that the double peak in the emission lines in AGN can be also due to a wide range of phenomena, like outflows, inflows and unresolved rotation of the gas surrounding the SMBH. Recently, thanks to integral field unit spectrographs, it has been revealed with high detail that the complexity of the emission line profile can be attributed to these phenomena in most of the cases (e.g. Roche et al. 2016; Davies et al. 2017). Therefore, using the doubly peaked feature alone does not guarantee that the target is a dual AGN and a multi-wavelength approach is necessary to confirm the binary system. Complementary observations can be perfomed at X-rays, because the two SMBH should exhibit non-thermal X-ray emission and, therefore, are easy to recognize at these wavelengths (Lena et al. 2018). However, the limited angular resolution of X-ray instruments does not allow the identification of the closest pairs of AGN. If the two AGN are radio emitters, the high angular resolution of radio interferometers can spatially resolve the system. Therefore, radio interferometric observations provide one of the most direct methods to identify dual AGN (Burke-Spolaor et al. 2018).

In this context, gravitational lensing eases the confirmation of such close binary SMBH systems. The magnifying effect of gravitational lensing can spatially resolve the two AGN, especially if they are located in the region at highest magnification, namely close to the caustics. However, the gravitational lensing effect is a rare phenomenon, as the probability of observing a multiply-imaged quasar is of the order 10−3 (e.g. Turner et al. 1984). Therefore, detecting a gravitationally lensed dual AGN source is expected to be extremely unlikely.

In this paper, we investigate the gravitational lensing system MG B2016+112, whose peculiar properties have been puzzling since its discovery (e.g. Garrett et al. 1994). In particular, we compare two VLBI observations at 1.7 GHz separated by 14.359 years with the aim of better understanding the nature of the background radio source. We detect a significant positional shift in the lensed images between the two epochs, which can be interpreted as either a proper motion along the jets or an orbital motion of two radio-loud AGN in the source plane. In Sect. 2, we introduce the radio properties of MG B2016+112. A description of the VLBI observations and data reduction is provided in Sect. 3. We present the lens modelling and source reconstruction in Sect. 4, while the discussion and conclusions are presented in Sects. 5 and 6, respectively.

Throughout this paper, we assume H0 = 67.8 km s−1 Mpc−1, ΩM = 0.31, and ΩΛ = 0.69 (Planck Collaboration VI 2018). The spectral index α is defined as Sν ∝ να, where Sν is the flux density as a function of frequency ν.

2. The target MG B2016+112

The gravitational lensing system MG B2016+112 was discovered during the MIT-Green Bank survey (MG survey, Lawrence et al. 1984; Bennett et al. 1986). It consists of three images (A, B and C) of a background source at redshift z = 3.2773, which is gravitationally lensed by an elliptical galaxy at redshift z = 1.01 and its satellite galaxy (Lawrence et al. 1993; Yamada et al. 2001; Soucail et al. 2001; Koopmans et al. 2002).

From the first VLBI observations of MG B2016+112 at 1.7 GHz, it was evident that images A and B are more compact, while image C is resolved into four sub-components connected by a faint extended emission (Lawrence et al. 1984; Garrett et al. 1994, 1996; Koopmans et al. 2002). Later, high sensitivity and high angular-resolution observations with global VLBI at 1.7 GHz up to 8 GHz revealed that images A and B are resolved into 5 sub-components, where some sub-components have a flat spectral energy distribution, while others show a steep radio spectrum (More et al. 2009). Also, region C is resolved into multiple sub-components with both flat and steep radio spectra, where the two closest images, C12 and C22, show both compact and extended emission (More et al. 2009). Thanks to the high angular resolution of these VLBI observations, it was possible to measure that there is a significant asymmetry in the angular separation of the sub-components of the merging images in region C. The lensed images C11–C12 and C21–C22 should show a mirror inverted morphology and, therefore, should have the same angular separation and a similar flux density. Such an astrometric anomaly can be considered as an indication for a mass density perturbation, which in this case was attributed to the presence of a spectroscopically confirmed satellite galaxy (G1) that is south of region C (Koopmans & Treu 2002; Chen et al. 2007; More et al. 2009).

The lensing galaxy is an elliptical galaxy (called D) with a stellar velocity dispersion of 328 ± 32 km s−1 (Koopmans & Treu 2002). Galaxy D is not active, as it does not display any emission at radio or X-ray wavelengths. This lensing galaxy lies in a massive galaxy cluster, which was detected for the first time at X-ray wavelengths (Chartas et al. 2001; Toft et al. 2003).

Several gravitational lens mass models have been proposed to reproduce the observations of MG B2016+112 (e.g. Narasimha et al. 1987; Nair & Garrett 1997). Using the image positions given by early European VLBI Network (EVN) observations at 5 GHz, Koopmans et al. (2002) proposed a model where all of the lensed images belong to the same background source, which is assumed to be a core-jet AGN. In this model, the caustics pass through the AGN in a way that the core is doubly-imaged (region A and B) and part of the counter-jet (region C) is quadruply imaged. This model was revised and confirmed by More et al. (2009) using follow-up global VLBI observations at 1.7, 5 and 8.46 GHz. The multiple sub-components detected in the lensed images provided more constraints to the mass model. Moreover, More et al. (2009) included the faint satellite galaxy in the model, which accounts for most of the astrometric anomaly observed in region C.

3. Multi-epoch VLBI imaging

In this section, we present the new and archival VLBI observations used to investigate the proper motion of the lensed radio components.

3.1. Observations

MG B2016+112 was observed with the ten 25 m antennas of the Very Long Baseline Array (VLBA) at a central frequency of 1.65 GHz (L-band) on 2016 July 2 (Epoch 2; ID: BS251, PI: Spingola). The experiment was carried out in phase-reference mode for a total observing time of 12 h. Scans on the phase-reference calibrator J2018+0831 of ∼2 min were alternated by observations of ∼3.5 min on the target (see Table 1 for details). The fringe finder and bandpass calibrator for this experiment was 3C454.3. The data were recorded at 2 Gbit s−1 and correlated at the VLBA correlator in Socorro to obtain two intermediate frequencies (IFs) of 128 MHz each, divided in 256 channels, at both circular RR and LL polarizations.

Table 1.

Summary of the VLBA observations at 1.7 GHz for MG B2016+112 at Epoch 1 and Epoch 2.

The archival global VLBI observations were taken on 2002 February 25 (Epoch 1; ID: GP0030, PI: Porcas), making the two epochs separated by Δt = 14.359 years. The setup of the observation is summarized in Table 1. The fringe finder for these observations was B2029+121, which was also the phase-reference calibrator. The observations were correlated to obtain 4 IFs of 8 MHz bandwidth each, which were divided in 16 spectral channels. The observing antennas for this experiment were Effelsberg, Jodrell Bank, Medicina, Onsala, Torun, Robledo, Goldstone, all the antennas of the VLBA and the phased Very Large Array (VLA). Further details of these observations are reported by More et al. (2009). For our analysis, we calibrate the VLBA-only subset of these observations, to match the uv-coverage between the two epochs.

3.2. Data reduction

We perform the calibration for both epochs following the standard VLBA procedures using the VLBAUTILS tasks built in the NRAO Astronomical Image Processing System (AIPS; Greisen 2003). The amplitude calibration is based on the a priori knowledge of the system temperature and gain curve of each antenna. The initial calibration steps include corrections for instrumental offsets, Earth rotation, atmospheric opacity, ionospheric dispersive delay and parallactic angle for the rotation of the antenna feed. Then, we correct for the fringe phases as a function of time and frequency (fringe fitting). Finally, we apply the bandpass calibration to correct for the response of the receiver as a function of frequency. All of these corrections are performed on the calibrators and then the solutions are interpolated to the target MG B2016+112.

The imaging and self-calibration of both observations have been performed within the Common Astronomy Software Applications package (CASA; McMullin et al. 2007). We apply phase-only self-calibration by starting with a solution interval as long as the scan length, which is iteratively decreased to 60 s. We do not use the first and last channels of the IFs. We use a Briggs weighting scheme during imaging (Robust = 0), which is a compromise between natural and uniform visibility weighting. The final restoring beam for the Epoch 2 observations is 11 mas × 5 mas at a position angle of 10°, while for the Epoch 1 (VLBA-only) observations is 11 mas × 4.5 mas at a position angle of 8°. Even if the difference in angular resolution between the two observations is small, it can lead to a possible incorrect identification of the Gaussian centroids of the sub-components of the lensed images at the two epochs. In order to suppress the angular resolution effects in the estimate of the lensed image positions, we use the same weighting scheme and restoring Gaussian clean beam for imaging the target at the two epochs. This step allows us to recover the same angular scales and, therefore, identify the same sub-components in the lensed images. Nevertheless, the non-identical uv-coverage of the two observations may also lead to differences in the imaging. To minimize possible effects due to the different uv-coverage, we use only the VLBA antennas for imaging Epoch 1.

The off-source rms noise level is ∼70 μJy beam−1 for the first epoch and ∼40 μJy beam−1 for the second epoch. This difference in sensitivity is to be expected, given the larger bandwidth of the new VLBA observations. The final self-calibrated VLBA images are shown in Fig. 1. The observations from both epochs clearly resolve image A into two sub-components (A1 and A2+A3), with an indication of two other sub-components (A4 and A5) that were detected in the global VLBI observations at 5 GHz by More et al. (2009). The sub-components of image B are more distorted and blended together with respect to image A. Moreover, image C is resolved into four distinct sub-components at both epochs, and shows a faint extended emission in the tangential direction that connects images C12 and C22 in the observations taken during Epoch 2 (see Fig. 1).

thumbnail Fig. 1.

Multi-epoch VLBA observations at 1.7 GHz of the gravitational lens MG B2016+112. The central image shows the entire system as observed at 1.7 GHz with VLBI. The black contours are the observations taken on 2002 February 25 (Epoch 1), the greyscale map and the red contours are the new observations taken on 2016 July 2 (Epoch 2). The greyscale map is in units of mJy beam−1, as indicated by the colour bar in each image. On VLBI-scales, image A is resolved into four components (A1, A2+A3, A4 and A5 following the nomenclature of More et al. 2009), image B is resolved into two components (B1 and B2+B3+B5) with an indication for a possible third component (B4), while image C is resolved into four components (C11, C12, C22 and C21). The shift is more visible in region C, which is at higher magnification (μ ∼ 270 to 350). Contours are at (0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.5, 0.75, 1) × the peak flux density of each individual image, which is ∼22 mJy beam−1 for Epoch 1 and 23 mJy beam−1 for Epoch 2. The restoring beam is shown in the bottom left corner and is 11 mas × 5 mas at a position angle of 10°. All of the images are obtained using a Briggs weighting scheme, with Robust = 0.

Open with DEXTER

Finally, the total flux density of the system and the flux density of each sub-component are in agreement within the errors between the two epochs, indicating no significant flux density variability in this period at 1.7 GHz.

3.3. Measurement of the lensed image positions

In order to measure the position of the lensed images, we fit the brightness distribution with two-dimensional Gaussian components using the task IMFIT within CASA. Images A and B are fitted with two Gaussian components, while the four sub-components of image C could be fitted by a single Gaussian component each. Note that image A is resolved into at least three sub-components (see Fig. 1). However, only two Gaussian components are clearly found when performing the fit (images A1 and A2+A3). As discussed in the previous section, the sub-components of image B are difficult to disentangle and clearly associate with the sub-components of image A. We use the Gaussian centroid as the position of the lensed images. The uncertainty on the position is estimated in the standard way, and depends on the major and minor axes of the elliptical Gaussian and the signal-to-noise ratio, under the assumption that the component is unresolved. This is a major assumption, considering the significant blending of some sub-components at this frequency, which affects especially the lensed image B. Therefore, the positional uncertainties estimated using this method should be considered as lower limits on the real ones. However, the two observations have almost the same uv-coverage, as we selected the same antennas (VLBA only) and the observing time is comparable. As a result, the morphology of the lensed images is found to be completely consistent between the two epochs (see Fig. 1). For this reason, any deviation from a single 2D Gaussian function in the lensed images would be alike in the two epochs.

3.4. Alignment of the two Epochs

Self-calibration was fundamental for recovering both the compact and the extended structure of image C, which is at low surface brightness. However, this resulted in the precise absolute coordinate information being lost. Moreover, as the two observations have different phase-reference calibrators, this also resulted in a different absolute position of the lensed images at the two epochs. For these reasons, we need to align the two images to a common reference frame. This was picked as the frame defined by the self-calibrated image of Epoch 2. Positions measured in Epoch 1 are then transformed to the Epoch 2 reference frame by means of a linear transformation that is computed based on the flux density-weighted positions of A1 and B1 in the two epochs. Two sources are enough for this purpose, as VLBI observations provide a distortion-free reference frame, with the two images already having the same scale and rotation1.

The position of the lensed images in the common reference frame and their uncertainties estimated with this method are listed in Table 2.

Table 2.

Position of the lensed images (Col. 1) of MG B2016+112 at Epoch 1 (Cols. 2 and 3) and Epoch 2 (Cols. 4 and 5) relative to the lensing galaxy, which is at (0, 0).

3.5. Positional offsets

We measure positional offsets in the range 0.4−5.7 mas in right ascension and 0.6−3 mas in declination between Epoch 1 and 2 for the various sub-components. These offsets are much larger than the astrometric uncertainties, which are between 8 and 30 μas for the group of sub-components associated with image C, and of the order of hundreds of μas for those making up images A and B (see Table 2). The lensed images with the largest positional offsets are C11, C12, C22 and C21, which are also at the highest magnification region. The positional shift of this group of images is clearly visible in the image plane, as shown in Fig. 1, and it is along the direction of highest magnification, which is a direct evidence for motion. Moreover, the direction of the shift is consistent with the symmetry expected by the gravitational lensing, that is, images C11–C12 and C21–C22 have moved in opposite directions along the highest magnification direction.

4. Lens modelling

We model MG B2016+112 using the software GRAVLENS and Monte Carlo realizations to estimate the errors on the mass model parameters and the source positions. To date, a change in the position of gravitationally lensed images over time has been only detected in two doubly-imaged radio sources. A tentative detection of positional change in the lensed images has been observed in JVAS B1030+074, where there is no clear correspondence among the source sub-components because the source lies in a region at low magnification (Zhang et al. 2007). While high frequency VLBI observations of PKS 1830–211 revealed a change in the lensed images position that has been attributed to a motion in the source (Jin et al. 2003). Therefore, MG B2016+112 represents one of the rare cases where proper motion has been clearly detected in a lensing system. Theoretically, the observed change in the image positions between Epoch 1 and 2 could be due to either a change in the lensing galaxy position (Birkinshaw 1989; Kochanek et al. 1996; Wucknitz & Sperhake 2004) or a movement of one (or more) radio components in the source (Jin et al. 2003; Biggs 2005). In this section, we explore both scenarios.

4.1. The lensing galaxies have moved

A possible explanation for the change in the image positions is that the lensing galaxy and/or the satellite galaxy have moved. For example, given that the lensing galaxy is in a cluster, with a velocity dispersion of σv ≃ 771 km s−1 (Soucail et al. 2001), we would expect a positional change of just ∼1.5 μas2 in the 14 year period between Epoch 1 and 2. Although small, such a change in the position of the caustics could result in a significant change in the position of the lensed images, particularly for image C.

In order to test how much the lensing galaxies could have moved from Epoch 1 to reproduce the image positions observed at Epoch 2, we proceed in the following way. By using the image positions at Epoch 1 and the lens mass model “scenario C” of More et al. (2009), we keep all the lens mass model parameters fixed and we use the lensing galaxies positions produced as follows. We generate random positions for the lensing galaxies within a radius of 0.058 mas, which is the distance travelled by the galaxy if it moved at the speed of light for the temporal baseline Δt between the two epochs, and is, therefore, a largely conservative assumption. Since we do not have any information on the direction of motion that the two lensing galaxies (the main galaxy and its satellite) could have, the circles where we generate the positions for the lensing galaxies are uniformly filled. Then, we forward ray-trace the source components to the image plane for all the simulated positions of the lensing galaxies and determine if the predicted positions of the lensed images match the observations at Epoch 2.

We find that the model-predicted positions cannot fit simultaneously images A and B, and image C; either the model reproduces the doubly-imaged source or the quadruply-imaged source, with offsets between the observed and the model-predicted positions of the order of 30σ on average; none of the simulated positions for the lensing galaxies can reproduce the position of the lensed images at Epoch 2. Therefore, we reject this scenario as a possible explanation for the positional shift observed in the lensed images between the two epochs.

4.2. The source has moved

The second scenario involves the source components in the source-plane moving with respect to the lensing galaxy position over the two epochs. In order to reconstruct the position of the source components, we again assume the mass density distribution proposed by More et al. (2009) as a starting model. In particular, we adopt the model where images A1–B1 and A2–B2 are doubly imaged, while the pairs C11–C21 and C12–C22 are quadruply imaged (scenario C). We choose this model because the morphology of region C and the separation between the sub-components is typical of a pair of merging images in a four-image system (e.g. Biggs et al. 2004, Hsueh et al. 2016). The main lensing galaxy (D) is parameterized as an elliptical power law mass density distribution [ρ(r)∝rγ]. The mass density distribution has 6 variables: the mass scale b; lens position (x, y); ellipticity e and its position angle, ϑ, and power-law slope γ. We keep the satellite galaxy (G1) fixed, assuming the same mass model of More et al. (2009). This model consists of a singular isothermal sphere (SIS) with mass strength b = 0.145 arcsec, at a position (xG1, yG1) = (0.840, −2.150) arcsec relative to the lensing galaxy D. We take into account the perturbation to the mass model due to the cluster of galaxies in the external shear term, which adds two other variables to the mass model, namely the shear strength Γ and its position angle Γϑ.

We simultaneously use the position of the lensed images listed in Table 2 to constrain the mass density distribution. In this way, we are implicitly assuming that the same mass density distribution reproduces the lensed images at the two epochs. This approach provides double the constraints to the lens model than using the two epochs separately, as each epoch provides an independent source distribution for the same lensing potential. Due to possible substructures within the lensing galaxy or along the line-of-sight, we do not use the relative flux-densities as additional constraints (Metcalf 2002; Mao et al. 2004; McKean et al. 2007; Rumbaugh et al. 2015; Hsueh et al. 2018; Despali et al. 2018). Since the counter-images of components A4 and A5 are not resolved in image B, we do not use these components as constraints.

The best model parameters are presented in Table 3, a schematic representation of the lens mass model is shown in Fig. 2 and the probability density distribution for each parameter is shown in Fig. 3. Our mass model did not deviate from the model proposed by More et al. (2009), which already converged to a global minimum of the χ2 function. In their model, More et al. (2009) fixed the ellipticity and position angle of the main deflector (D) to the values estimated from the surface brightness profile at near-infrared and optical wavelengths. Our model confirms that e and ϑ are consistent with the parameters derived from the stellar emission within 1σ. Therefore, there is a good alignment between the stellar and the dark matter components within the Einstein radius, which is generally not observed for lens systems with a strong external shear (Γ = 0.10 ± 0.02; Spingola et al. 2018; Shajib et al. 2019). We find the power-law slope to be consistent with an isothermal profile (γ = 2.0 ± 0.1), which is consistent with the results obtained by Treu & Koopmans (2002), who combined gravitational lensing and stellar dynamics (γTK2002 = 2.0 ± 0.1).

Table 3.

Best recovered lens model parameters for the mass density distribution of the main deflector (D).

thumbnail Fig. 2.

Convergence map for the lens mass model of MG B2016+112. The field-of-view is 6 arcsec × 6 arcsec. The white contours are the 1.7 GHz observations at Epoch 2, while the black contours are the critical and caustic curves. The filled magenta circle indicates the location of the source components. The white labels indicate the groups of lensed images (A, B and C), and the black labels identify the two lensing galaxies (D and G1), also shown by the convergence peaks.

Open with DEXTER

thumbnail Fig. 3.

Diagonal histograms: probability density functions (PDFs) for the lens model parameters of the main lensing galaxy (D), which were obtained by marginalizing over all the other model parameters, with two dashed vertical lines indicating the 1σ limits. The other panels show the 2-dimensional contours of the PDF for each pair of model parameters, where the contours indicate the 1σ region. The meaning of the parameters, their maximum-likelihood model value and their uncertainties are presented in Table 3.

Open with DEXTER

Some of the model-predicted positions of the lensed images were found to differ from the observed positions by 2 to 10σ, which was also noted by More et al. (2009). These positional residuals are not as critical as, for example, in the cases of CLASS B0128+437 and MG J0751+2716 (Biggs et al. 2004; Spingola et al. 2018). Therefore, the astrometric anomaly in MG B2016+112 is not completely solved by the inclusion of G1, but could be due to an extra mass component that is currently not part of the model (e.g. see Spingola et al. 2018 for discussion). More et al. (2009) also tested a model with three lensing galaxies, but found that this did not improve the model-predicted positions of the lensed images. Therefore, a more complex model for the mass density distribution is needed to fully explain the image positions of MG B2016+112.

Our best model predicts the position and flux density of the counter-images of region C, at the position of region A and B, finding a flux density between 2 and 10 μJy for the image pair C11–C21, and less than 1 μJy for the image pair C12–C22. These flux densities are lower by at least a factor of two when compared to the rms noise level of our imaging data. Therefore, the non-detection of these counter images in regions A and B is consistent with our best model, but future deeper observations that detect these counter images are needed to test the validity of the mass model.

In Fig. 4, we show the reconstructed source plane, given our best mass model, and we list the position of the source components in Table 4. Source 1 corresponds to images A2–B2, source 2 corresponds to images A1–B1, source 3 corresponds to images C11–C21 and source 4 corresponds to images C12–C22. The uncertainty on each source position is estimated via Monte Carlo realizations using the following procedure. We simulate 1000 lensed images by randomly extracting them from a Gaussian distribution with a standard deviation that corresponds to the observed uncertainty and the expectation value of the observed position. We then keep our best lens mass model fixed and use these simulated lensed images to perform backward ray-tracing to obtain 1000 realizations for each source. Finally, the uncertainty on each component is given by the standard deviation of the 1000 source positions. In this way, the magnification (μ) is taken into account in the estimate of the uncertainty, as opposed to just the uncertainty in the observed image position. As a result, regions with a higher magnification will have a lower positional uncertainty. Indeed, the source components associated with the highly magnified region C (sources 3 and 4; μ ∼ 270 and ∼350, respectively) have a positional uncertainty of the order of 8 to 20 μas, while the astrometric uncertainty on the source components corresponding to images A and B (sources 1 and 2; μ ∼ 2 for both the sources) is between 1 and 4 mas (see Table 4). The positional uncertainty is larger in declination than in right ascension, which reflects the shape of the synthesized beam (see Fig. 1). We highlight that source 1 has the largest positional uncertainty, not only because of the low μ, but also because it is associated with the lensed image components A2–B2, which are the most difficult images to de-blend in the image plane, especially for image B (see Fig. 1).

thumbnail Fig. 4.

Source-plane model for the VLBA observations of MG B2016+112. Epoch 1 (blue) and 2 (red) observations aligned on the lensing galaxy position, which is at (0, 0). The solid blue line represents the caustics, which divides the source-plane into the double- and quadruple-image regions. The position of the sources is indicated by the filled circles. Source 1 corresponds to images A2–B2, source 2 corresponds to images A1–B1, source 3 corresponds to images C11–C21, and source 4 corresponds to images C12–C22. The labels “core” and “jet” are based on the radio spectral energy distribution of each image pair, as reported by More et al. (2009). The grey scale bar at the bottom left corner represents 100 pc at redshift z = 3.2773. The error bars take into account the uncertainties in the lens model. In the case of sources 1 and 2 (related to the lensed images A and B), the error is dominated by our ability to de-convolve the different source components, while for sources 3 and 4 (related to the lensed region C), the errors are very small because the sub-components are clearly spatially resolved due to their extremely large magnification (μ ∼ 270 to 350).

Open with DEXTER

Table 4.

Properties of the source-plane components of MG B2016+112.

Our source reconstruction finds that components 3 and 4 have moved in the same direction (south) by ∼40 ± 25 μas, while component 2 has shifted by ∼2 ± 1 mas in the north-east direction. As shown in Fig. 4, the positional shift of source component 2 is statistically significant only in the right ascension direction. The positional shift in source components 3 and 4 is significant in both directions. We would like to emphasize that the lack of a significant change in flux density between the two epochs is consistent with our source reconstruction. This is particularly important for the source components at high magnification, as they are the most sensitive to any small change in the lensing configuration, due to the steep rise of the magnification when going close to the caustics. For example, if the source components associated with C11–C21 and C12–C22 moved towards (away from) the caustics, their flux density would have significantly increased (decreased) at the second epoch. Therefore, their measured flux densities are an indirect evidence that the motion of such components must have been parallel to the caustics. Finally, we note that any error on the lens model will not absorb the proper motion, but will modify the relative positions in the source-plane.

5. Discussion

We have found evidence for proper motion in the lensed images of MG B2016+112 by analyzing two VLBI observations at 1.7 GHz that are separated by 14.359 years (see Fig. 1). In Sect. 4, we ruled out the possibility that the proper motion is due to a shift in the lensing galaxy position, and we attribute it to a motion in the source. The source-plane reconstruction (see Fig. 4) shows that the de-lensed radio-loud object is quite complex, with two sets of two components moving in different directions. In this section, we investigate two possible interpretations for explaining the source morphology. First, we will explore the scenario where all of the sub-components belong to the same AGN. In this case, the motion is attributed to knots moving along the jets. Second, we will examine the possibility of having two separate radio-loud AGN in the source plane that are interacting with each other.

5.1. Single AGN scenario

Most jetted AGN show only one jet (Padovani et al. 2017). This is due to relativistic boosting, which enhances the radiation in the forward direction due to an approaching jet, and reduces the emission in the backward direction due to a receding jet (Scheuer & Readhead 1979). However, according to the unified AGN model, if the jet and counter-jet are seen under a large viewing angle, it is possible to detect both of them, as for example in FR I type radio galaxies (Urry & Padovani 1995). In these cases, it is expected that the counter-jet moves at sub-luminal velocities in an opposite direction with respect to the approaching jet, as seen for example in Centaurus A (Jones et al. 1996). MG B2016+112 shows two optically thin components, which can be potentially associated with a jet and a counter-jet (sources 2 and 3, respectively), as one is moving at super-luminal velocity (v2 = 2.9c ± 3.9c), while the other has a sub-luminal velocity (v3 = 0.06c ± 0.04c) at an angular diameter distance of DA = 1576.2 Mpc.

To test this scenario, we use the apparent motion of these two optically thin source components (sources 2 and 3) to estimate the theoretical ratio between the flux density of the possible jet and counter jet. This value can then be compared with the intrinsic flux density ratio between the source components associated with the jet and counter-jet, namely the flux densities corrected for the magnification. Since source component A1–B1 moves at a higher velocity than component C11–C21, we assume that source 2 consists of a knot in the approaching jet, while source 3 could be a knot in the receding counter-jet.

The ratio between the jet and counter jet flux densities can be written as

(1)

where α is the spectral index, β is the velocity expressed in units of c and θ is the viewing angle (Scheuer & Readhead 1979). By assuming an intrinsically symmetric ejection of the two optically thin components, the factor βcos(θ) can be expressed in terms of the proper motion of the jet (μj) and counter-jet (μcj),

(2)

(Fender et al. 1999).

From Eq. (2), assuming β = 1 and α = 0.83, we find a maximum viewing angle of θmax ≃ 17°, and a theoretical flux density ratio of R ≃ 37 500. However, the observed flux density ratio, when corrected for the lensing magnification, between the jet (A1–B1) and the counter-jet (C11–C21) is ∼270, which is two orders of magnitude less than the predicted ratio. However, this criterion is based on the strong assumption of a symmetric ejection of the knots in the jet and counter-jet. Moreover, given the large light travel time between jet and counter-jet, and the likely not simultaneous changes in the two radio ouflows, our R value should be taken only as an indication that the single AGN scenario may not be completely compatible with the observations, rather than a conclusive statement. Also, the projected direction of the motion indicates that the two flat-spectrum components (i.e. the cores; sources 1 and 4) are moving perpendicularly to each other (see Fig. 4), even though the positional uncertainty is large for source component 1. This motion would imply an exotic jetted-AGN, or a possible reverse shock in the emission at pc-scales, as observed for example in the powerful radio jets of M 87 and 3C345 (Unwin & Wehrle 1992).

5.2. Dual AGN scenario

The multi-wavelength properties of MG B2016+112, when taken together, are also consistent with a possible dual AGN (DAGN) interpretation (defined as a pair of AGN separated by less than 10 kpc, while binary AGN consists of a pair of SMBH that are separated by less than < 100 pc, Burke-Spolaor et al. 2014). DAGN show specific morphological and spectral features, such as multiple flat-spectrum radio-cores and misaligned/disturbed kpc-scale jets with a S- or X- shaped morphology (e.g. Deane et al. 2014; Burke-Spolaor et al. 2018); jet-dominated radio emission (Frey et al. 2012; An et al. 2018); double peaked optical spectral emission lines separated by a few hundred km s−1 (e.g. Hβ, [OIII], Comerford et al. 2009; Liu et al. 2018); and multiple X-ray point source components (Koss et al. 2012).

5.2.1. Evidence in favour of the DAGN scenario from previous studies

Yamada et al. (2001) showed that the narrow-line spectra of images B and C have different properties. These spectra, obtained with the Canada-France-Hawaii Telescope (CFHT), show typical emission lines from active galaxies (e.g. C III, C IV, He II, N V). Yamada et al. (2001) found that they could not fit a single photo-ionization model that could explain simultaneously the line ratios with He II and those with N V for both images B and C. This led to the interpretation that the excitation between these different parts of the background source is also different, concluding that MG B2016+112 is likely a partially dust-obscured low-luminosity narrow-line AGN.

These differences in the optical spectra of images B and C could be due to two separate narrow-line regions, one associated with an un-obscured AGN (images A and B, which show a quasar morphology at optical wavelengths) and the other associated with a dust-obscured AGN (image C, which has an extended optical morphology). If so, the same emission lines should show a velocity offset. However, the low spectral-resolution of the CFHT observations does not provide an accurate enough measurement of the relative velocities of the narrow lines in images B and C, and the line properties of image A are still unknown. Alternatively, as region C is close to the caustics, the difference in the emission line flux-ratios could be due to a large differential magnification across a complex narrow line region. Therefore, even though the different line flux-ratios could be interpreted as evidence for a DAGN, further observations to measure the relative line velocities and their positions relative to the lensing caustics are needed.

Observations with Chandra showed that all three lensed images are X-ray sources (Hattori et al. 1997; Chartas et al. 2001). When correcting for the distortion due to gravitational lensing, the source corresponding to images A and B has a 2–10 keV luminosity of 3 × 1043–1.4 × 1044 erg s−1, but the authors do not investigate the intrinsic X-ray properties of the source related to image C. The images are quite faint in X-rays, with only 6, 5 and 12 photon counts for images A, B and C, respectively (Chartas et al. 2001).

The detection of multiple X-ray components associated with images A and B, and image C (Chartas et al. 2001), which also may have different intrinsic luminosities, could be explained by the presence of two possible distinct accretion disks within MG B2016+112. Nevertheless, Chartas et al. (2001) explain these differences in the X-ray properties as images A and B being associated with the AGN, and the emission from region C being related to inverse Compton emission associated with the radio jets. However, from the current data, it is not clear which interpretation is correct for the X-ray emission from MG B2016+112.

Based on the radio spectral energy distribution (More et al. 2009), there is evidence of two flat-spectrum components and two steep-spectrum components. Classically, the flat-spectrum component is considered the core (i.e. the emission at the base of the jet, closest to the black hole), while the steep-spectrum component consists of the jet(s) of the AGN. Therefore, there are two possible cores and two possible jets in the source plane of MG B2016+112. These two candidate core-jet AGN are intrinsically faint, with flux densities of the individual sub-components between 0.01 and 10 mJy. These properties at radio wavelengths can be taken as evidence in favour of the DAGN scenario.

5.2.2. Evidence from proper motion

The measurement of proper motion, and the direction of this proper motion in the source plane for the two different parts of the background source can also be taken as evidence for the DAGN interpretation. The source-plane consists of four components; sources 1 and 4 are the two flat-spectrum components (α ∼ 0.2 between 1.7 and 5 GHz), while sources 2 and 3 have a steeper spectral index (α ∼ 0.8 between 1.7 and 5 GHz; More et al. 2009). Given their relative projected distance in our reconstructed source-plane (see Fig. 4), they seem to form two separate core-jet structures. Therefore, we associate sources 1 and 4 with candidate radio cores, while sources 2 and 3 are identified as candidate jet components, as discussed briefly in the previous section. The separation between the two core-jet AGN is about 175 pc in projection, which is a strong indication that the two objects should be gravitationally bound, potentially forming a DAGN. The relative position of the optically thin components seems to indicate a misalignment between the radio jets. The presence of two flat-spectrum components and multiple misaligned jets is generally a criterion used for identifying DAGN at radio and X-ray wavelengths (Owen et al. 1985; Lal & Rao 2007), making this morphology consistent with the DAGN scenario. Clearly, more precise positional measurements of the source components 1 and 2 are needed to confirm the differences between the jet alignment.

The most unusual feature of the core component associated with images C12–C22 is the proper motion (source 4; see Fig. 4). Generally, the core is stationary in single AGN galaxies (e.g. Marscher 2009). Therefore, a movement of the radio core, as may be seen here, would imply a shift of the entire AGN system. Moreover, the two optically thin components (especially source 3) are moving in a similar direction as their associated core components. This could be due to the two AGN dragging their jets while they move. Sources 3 and 4 are found to be moving with an apparent sub-luminal velocity in the southern direction, with v3 = 18900 ± 15000 km s−1 and v4 = 20100 ± 13000 km s−1, for the candidate jet and core, respectively. Source 1 does not show significant motion within the uncertainties (see Table 4 and Fig. 4), while source 2 is moving with an apparent super-luminal velocity of v2 = 2.9c ± 3.9c. This velocity indicates the presence of Doppler boosting, and hence, requires the jet to be oriented at a small viewing angle. Therefore, the motion of this component might be a combination of the proper motion of the entire AGN, and the motion of the optically thin outflow with respect to the main core. Even if the velocities of the two AGN are higher than those typical of galaxies in clusters (< 1000 km s−1), our findings are still consistent with a scenario where the two AGN candidates are at a later stage of the merging and, therefore, accelerating. This is even more plausible given their small angular separation.

We find that both candidates AGN show low intrinsic flux densities, but they have different radio properties (More et al. 2009). The flux density of the possible AGN comprising sources 1 and 2 is dominated by the emission from the jet, whereas the candidate AGN composed of sources 3 and 4 is core-dominated. This kind of difference between the two radio-loud SMBH can be attributed to the different orientations of the two interacting AGN, which may be further evidence in favour of the DAGN scenario.

5.3. Implications of a dual AGN scenario

To date, two main avenues have been proposed for the formation of SMBHs in galaxies: the accretion of gas from a directly collapsed star (Begelman 2002) or the merging of multiple black hole seeds (Volonteri et al. 2003). The presence of a DAGN in MG B2016+112 would be in favour of the merger-driven formation scenario, as DAGN represent an intermediate evolutionary stage of such a process. According to the hierarchical formation scenario, this is expected to be observed at high redshift (Volonteri et al. 2016). The timescales on which multiple SMBHs can coalesce are not known, but it is expected to be short given the low detection rate of such systems. Therefore, observations of small separation DAGN are needed in order to probe the final stages of the merging process.

As the lensing probability of radio-loud AGN is ∼10−3 and the expected fraction of AGN with a dual SMBH system is ∼10−2 at high redshift, the combined probability of detecting a gravitationally lensed DAGN system is of order 10−5. Therefore, the detection of a gravitationally lensed DAGN here would imply that the overall probability is higher than expected, meaning that DAGN are more common in the early Universe than first thought. As the MG survey detected ∼6000 sources at signal-to-noise ratio above 5 (Bennett et al. 1986) and found six strong gravitational lensing systems, this gives approximately a strong lensing probility of 10−3, which is a typical strong lensing rate (e.g. Spingola et al. 2019). If MG B2016+112 is a genuine DAGN, then we can roughly assume that 1 every 6 strongly lensed sources is a DAGN, resulting in a probability of detecting a gravitationally lensed DAGN of 0.16 percent. Therefore, the probability of finding a gravitationally lensed DAGN system in the MG survey is ∼2 × 10−4, which is an order of magnitude higher than the expectations from the current simulations.

To some extent, this is expected given the increased merger rate (expected at high redshift), but to catch one in the act of forming would also mean that the in-spiral time needs to be slower than what is currently predicted (e.g. Rafikov 2016). As this conclusion is currently based on the statistics of this single possible detection, further studies of other lensed radio sources at high angular resolution with VLBI are needed to determine if there are other cases with proper motion in the radio jets, and whether such motion is also consistent with a DAGN system.

6. Conclusions and future work

In this paper, we have presented a clear detection of proper motion from a gravitationally lensed radio source at high redshift. From analyzing two VLBI datasets separated by 14.359 years that were taken with the same array and at the same frequency, we detected shifts up to 6 mas in the position of the lensed images. We test the possibility that the cause of these shifts is due to a motion of the lensing galaxies, which we find is unlikely. Therefore, we conclude that the observed positional shift seen in the lensed images is due to proper motion in the source plane, which we could reconstruct with a precision down to about 1 μas yr−1. Such an outstanding precision is more than a factor of 10 better than that of the best optical proper motions currently available (Massari et al. 2018; Libralato et al. 2018), and demonstrate the power of combining radio interferometry with gravitational lensing techniques.

To explain the observed proper motion, we investigated two possible scenarios. Assuming that the source consists of a single AGN, a possible explanation for the proper motion is given by the movement of knots moving along the radio jets. However, the de-magnified flux densities of the components are apparently not consistent with knots moving along an approaching and a receding jet. The second and more exotic scenario consists of two radio-loud AGN separated by ∼175 pc in projection, both with a core-jet morphology, which form a DAGN system. In this scenario, which is mainly driven by the motion of the flat-spectrum radio components, the two core-jet AGN are moving relative to each other and the jet components are misaligned. If genuine, identifying a DAGN at redshift 3 would have important implications for our understanding of galaxy formation at high redshift, as it would be evidence in favour of the merging scenario for the formation of SMBHs.

The relative position of the candidate DAGN in MG B2016+112 depends on the lens mass model. Therefore, any error in the lens model translates to an incorrect estimate of the proper motion. Our lens mass model indicates that there is the presence of an astrometric anomaly, even when the companion satellite galaxy is explicitly taken into account. This implies that the mass density distribution is more complex than the model presented here, which includes the main lensing galaxy and its closest satellite galaxy. For example, a Bayesian grid-based analysis (e.g. Dye & Warren 2005; Koopmans 2005; Vegetti & Koopmans 2009) can help in understanding whether this parametric model is overly simplified. Moreover, by performing pixelated potential corrections to the smooth potential associated with the main lensing galaxy, it will be possible to quantify the mass of the satellite galaxy and better constrain the substructure mass density profile, which in our model is fixed (Vegetti & Vogelsberger 2014). Also, modelling simultaneously the multi-wavelength extended emission from MG B2016+112 will give many more observational constraints to the mass model (Suyu et al. 2012).

Together with a more sophisticated lens mass modelling, we also need additional observations to further constrain the source scenarios for MG B2016+112. Further radio observations at higher angular resolution and dynamic range will improve the precision of the image positions, as they would resolve all the subcomponents in the lensed images, and, therefore, better determine the relative motions of sources 1 and 2. Moreover, the detection of the counter-images of region C is a direct test to the validity of the lens mass model presented here. Also, optical spectroscopic observations can potentially reveal velocity offsets between the optical emission lines associated with the AGN activity (e.g. Lyα) in the different lensed images, which would add to the case for the DAGN scenario. Spectroscopic observations at high spectral resolution would also clearly discern between the presence of two narrow line regions or multiple photo-ionization levels within a single complex narrow line region associated with MG B2016+112.


1

We also searched for NVSS sources within the field-of-view that could be used for aligning the two images. However, the two closest NVSS sources are at > 6 arcmin from the phase center and the distortion due to bandwith and time smearing is too strong to make them reliable astrometric references.

2

As the proper motion μ is the distance traveled by the object divided by the time difference between the two epochs, the positional change (in arcsec) between the two epochs corresponds to v * Δt/(4.74 * D); where v = 771 km s−1, Δt = 14.359 years and D = 1.576 × 106 pc.

3

This is the average spectral index between 1.7 and 5 GHz for the optically thin components (More et al. 2009).

Acknowledgments

We thank the anonymous referee for the useful suggestions on this manuscript. CS would like to thank Z. Paragi, R. Schulz and R. Morganti for the useful discussions about this work. CS thanks JIVE for their help with the data reduction of the Epoch 1 observations. CS and JPM acknowledge support from a NWO-CAS grant (project number 629.001.023). DM acknowledges financial support from a Vici grant from NWO. LVEK is supported through an NWO-VICI grant (project number 639.043.308). The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code: GP0030. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

  1. An, T., Mohan, P., & Frey, S. 2018, Radio Sci., 53, 1211 [NASA ADS] [CrossRef] [Google Scholar]
  2. Begelman, M. C. 2002, ApJ, 568, L97 [NASA ADS] [CrossRef] [Google Scholar]
  3. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bennett, C. L., Lawrence, C. R., Burke, B. F., Hewitt, J. N., & Mahoney, J. 1986, ApJS, 61, 1 [NASA ADS] [CrossRef] [Google Scholar]
  5. Biggs, A. 2005, in Future Directions in High Resolution Astronomy, eds. J. Romney, & M. Reid, ASP Conf. Ser., 340, 433 [NASA ADS] [Google Scholar]
  6. Biggs, A. D., Browne, I. W. A., Jackson, N. J., et al. 2004, MNRAS, 350, 949 [NASA ADS] [CrossRef] [Google Scholar]
  7. Birkinshaw, M. 1989, in Gravitational Lenses, eds. J. M. Moran, J. N. Hewitt, & K. Y. Lo (Berlin: Springer Verlag), Lect. Notes Phys., 330, 59 [NASA ADS] [CrossRef] [Google Scholar]
  8. Burke-Spolaor, S., Brazier, A., Chatterjee, S., et al. 2014, ArXiv e-prints [arXiv:1402.0548] [Google Scholar]
  9. Burke-Spolaor, S., Blecha, L., Bogdanović, T., et al. 2018, in Science with a Next Generation Very Large Array, ed. E. Murphy, ASP Conf. Ser., 517, 677 [NASA ADS] [Google Scholar]
  10. Chartas, G., Bautz, M., Garmire, G., Jones, C., & Schneider, D. P. 2001, ApJ, 550, L163 [NASA ADS] [CrossRef] [Google Scholar]
  11. Chen, J., Rozo, E., Dalal, N., & Taylor, J. E. 2007, ApJ, 659, 52 [NASA ADS] [CrossRef] [Google Scholar]
  12. Comerford, J. M., Griffith, R. L., Gerke, B. F., et al. 2009, ApJ, 702, L82 [NASA ADS] [CrossRef] [Google Scholar]
  13. Davies, R. L., Groves, B., Kewley, L. J., et al. 2017, MNRAS, 470, 4974 [NASA ADS] [CrossRef] [Google Scholar]
  14. Deane, R. P., Paragi, Z., Jarvis, M. J., et al. 2014, Nature, 511, 57 [NASA ADS] [CrossRef] [Google Scholar]
  15. Despali, G., Vegetti, S., White, S. D. M., Giocoli, C., & van den Bosch, F. C. 2018, MNRAS, 475, 5424 [NASA ADS] [CrossRef] [Google Scholar]
  16. Dye, S., & Warren, S. J. 2005, ApJ, 623, 31 [NASA ADS] [CrossRef] [Google Scholar]
  17. Enoki, M., Inoue, K. T., Nagashima, M., & Sugiyama, N. 2004, ApJ, 615, 19 [NASA ADS] [CrossRef] [Google Scholar]
  18. Fender, R. P., Garrington, S. T., McKay, D. J., et al. 1999, MNRAS, 304, 865 [NASA ADS] [CrossRef] [Google Scholar]
  19. Frey, S., Paragi, Z., An, T., & Gabányi, K. É. 2012, MNRAS, 425, 1185 [NASA ADS] [CrossRef] [Google Scholar]
  20. Garrett, M. A., Muxlow, T. W. B., Patnaik, A. R., & Walsh, D. 1994, MNRAS, 269, 902 [NASA ADS] [CrossRef] [Google Scholar]
  21. Garrett, M. A., Porcas, R. W., Nair, S., & Patnaik, A. R. 1996, MNRAS, 279, L7 [NASA ADS] [CrossRef] [Google Scholar]
  22. Greisen, E. W. 2003, in Information Handling in Astronomy – Historical Vistas, ed. A. Heck, Astrophys. Space Sci. Lib., 285, 109 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hattori, M., Ikebe, Y., Asaoka, I., et al. 1997, Nature, 388, 146 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hsueh, J.-W., Fassnacht, C. D., Vegetti, S., et al. 2016, MNRAS, 463, L51 [NASA ADS] [CrossRef] [Google Scholar]
  25. Hsueh, J.-W., Despali, G., Vegetti, S., et al. 2018, MNRAS, 475, 2438 [NASA ADS] [CrossRef] [Google Scholar]
  26. Jin, C., Garrett, M. A., Nair, S., et al. 2003, MNRAS, 340, 1309 [NASA ADS] [CrossRef] [Google Scholar]
  27. Jones, D. L., Tingay, S. J., Murphy, D. W., et al. 1996, ApJ, 466, L63 [NASA ADS] [CrossRef] [Google Scholar]
  28. Kochanek, C. S., Kolatt, T. S., & Bartelmann, M. 1996, ApJ, 473, 610 [NASA ADS] [CrossRef] [Google Scholar]
  29. Koopmans, L. V. E. 2005, MNRAS, 363, 1136 [NASA ADS] [CrossRef] [Google Scholar]
  30. Koopmans, L. V. E., & Treu, T. 2002, ApJ, 568, L5 [NASA ADS] [CrossRef] [Google Scholar]
  31. Koopmans, L. V. E., Garrett, M. A., Blandford, R. D., et al. 2002, MNRAS, 334, 39 [NASA ADS] [CrossRef] [Google Scholar]
  32. Koss, M., Mushotzky, R., Treister, E., et al. 2012, ApJ, 746, L22 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lal, D. V., & Rao, A. P. 2007, MNRAS, 374, 1085 [NASA ADS] [CrossRef] [Google Scholar]
  34. Lawrence, C. R., Schneider, D. P., Schmidt, M., et al. 1984, Science, 223, 46 [NASA ADS] [CrossRef] [Google Scholar]
  35. Lawrence, C. R., Neugebauer, G., & Matthews, K. 1993, AJ, 105, 17 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lena, D., Panizo-Espinar, G., Jonker, P. G., Torres, M. A. P., & Heida, M. 2018, MNRAS, 478, 1326 [NASA ADS] [CrossRef] [Google Scholar]
  37. Libralato, M., Bellini, A., van der Marel, R. P., et al. 2018, ApJ, 861, 99 [NASA ADS] [CrossRef] [Google Scholar]
  38. Liu, X., Guo, H., Shen, Y., Greene, J. E., & Strauss, M. A. 2018, ApJ, 862, 29 [NASA ADS] [CrossRef] [Google Scholar]
  39. Mao, S., Jing, Y., Ostriker, J. P., & Weller, J. 2004, ApJ, 604, L5 [NASA ADS] [CrossRef] [Google Scholar]
  40. Marscher, A. P. 2009, ArXiv e-prints [arXiv:0909.2576] [Google Scholar]
  41. Massari, D., Breddels, M. A., Helmi, A., et al. 2018, Nat. Astron., 2, 156 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mayer, L., Kazantzidis, S., Madau, P., et al. 2007, Science, 316, 1874 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  43. McKean, J. P., Koopmans, L. V. E., Flack, C. E., et al. 2007, MNRAS, 378, 109 [NASA ADS] [CrossRef] [Google Scholar]
  44. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), ASP Conf. Ser., 376, 127 [Google Scholar]
  45. Metcalf, R. B. 2002, ApJ, 580, 696 [NASA ADS] [CrossRef] [Google Scholar]
  46. More, A., McKean, J. P., More, S., et al. 2009, MNRAS, 394, 174 [NASA ADS] [CrossRef] [Google Scholar]
  47. Nair, S., & Garrett, M. A. 1997, MNRAS, 284, 58 [NASA ADS] [CrossRef] [Google Scholar]
  48. Narasimha, D., Subramanian, K., & Chitre, S. M. 1987, ApJ, 315, 434 [NASA ADS] [CrossRef] [Google Scholar]
  49. Owen, F. N., O’Dea, C. P., Inoue, M., & Eilek, J. A. 1985, ApJ, 294, L85 [NASA ADS] [CrossRef] [Google Scholar]
  50. Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&ARv, 25, 2 [NASA ADS] [CrossRef] [Google Scholar]
  51. Planck Collaboration VI. 2018, A&A, submitted [arXiv:1807.06209] [Google Scholar]
  52. Rafikov, R. R. 2016, ApJ, 827, 111 [NASA ADS] [CrossRef] [Google Scholar]
  53. Roche, N., Humphrey, A., Lagos, P., et al. 2016, MNRAS, 459, 4259 [NASA ADS] [CrossRef] [Google Scholar]
  54. Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49 [NASA ADS] [CrossRef] [Google Scholar]
  55. Rosas-Guevara, Y. M., Bower, R. G., McAlpine, S., Bonoli, S., & Tissera, P. B. 2019, MNRAS, 483, 2712 [NASA ADS] [CrossRef] [Google Scholar]
  56. Rumbaugh, N., Fassnacht, C. D., McKean, J. P., et al. 2015, MNRAS, 450, 1042 [NASA ADS] [CrossRef] [Google Scholar]
  57. Scheuer, P. A. G., & Readhead, A. C. S. 1979, Nature, 277, 182 [NASA ADS] [CrossRef] [Google Scholar]
  58. Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649 [NASA ADS] [CrossRef] [Google Scholar]
  59. Soucail, G., Kneib, J.-P., Jaunsen, A. O., et al. 2001, ApJ, 367, 741 [Google Scholar]
  60. Spingola, C., McKean, J. P., Auger, M. W., et al. 2018, MNRAS, 478, 4816 [NASA ADS] [CrossRef] [Google Scholar]
  61. Spingola, C., McKean, J. P., Lee, M., Deller, A., & Moldon, J. 2019, MNRAS, 483, 2125 [NASA ADS] [CrossRef] [Google Scholar]
  62. Suyu, S. H., Hensel, S. W., McKean, J. P., et al. 2012, ApJ, 750, 10 [NASA ADS] [CrossRef] [Google Scholar]
  63. Toft, S., Soucail, G., & Hjorth, J. 2003, MNRAS, 344, 337 [NASA ADS] [CrossRef] [Google Scholar]
  64. Treu, T., & Koopmans, L. V. E. 2002, ApJ, 575, 87 [NASA ADS] [CrossRef] [Google Scholar]
  65. Turner, E. L., Ostriker, J. P., & Gott, J. R., III. 1984, ApJ, 284, 1 [NASA ADS] [CrossRef] [Google Scholar]
  66. Unwin, S. C., & Wehrle, A. E. 1992, ApJ, 398, 74 [NASA ADS] [CrossRef] [Google Scholar]
  67. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  68. Vegetti, S., & Koopmans, L. V. E. 2009, MNRAS, 392, 945 [NASA ADS] [CrossRef] [Google Scholar]
  69. Vegetti, S., & Vogelsberger, M. 2014, MNRAS, 442, 3598 [NASA ADS] [CrossRef] [Google Scholar]
  70. Volonteri, M. 2010, A&ARv, 18, 279 [NASA ADS] [CrossRef] [Google Scholar]
  71. Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559 [NASA ADS] [CrossRef] [Google Scholar]
  72. Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979 [NASA ADS] [CrossRef] [Google Scholar]
  73. Wucknitz, O., & Sperhake, U. 2004, Phys. Rev. D, 69, 063001 [NASA ADS] [CrossRef] [Google Scholar]
  74. Yamada, T., Yamazaki, S., Hattori, M., Soucail, G., & Kneib, J.-P. 2001, ApJ, 367, 51 [Google Scholar]
  75. Zhang, M., Jackson, N., Porcas, R. W., & Browne, I. W. A. 2007, MNRAS, 377, 1623 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Summary of the VLBA observations at 1.7 GHz for MG B2016+112 at Epoch 1 and Epoch 2.

Table 2.

Position of the lensed images (Col. 1) of MG B2016+112 at Epoch 1 (Cols. 2 and 3) and Epoch 2 (Cols. 4 and 5) relative to the lensing galaxy, which is at (0, 0).

Table 3.

Best recovered lens model parameters for the mass density distribution of the main deflector (D).

Table 4.

Properties of the source-plane components of MG B2016+112.

All Figures

thumbnail Fig. 1.

Multi-epoch VLBA observations at 1.7 GHz of the gravitational lens MG B2016+112. The central image shows the entire system as observed at 1.7 GHz with VLBI. The black contours are the observations taken on 2002 February 25 (Epoch 1), the greyscale map and the red contours are the new observations taken on 2016 July 2 (Epoch 2). The greyscale map is in units of mJy beam−1, as indicated by the colour bar in each image. On VLBI-scales, image A is resolved into four components (A1, A2+A3, A4 and A5 following the nomenclature of More et al. 2009), image B is resolved into two components (B1 and B2+B3+B5) with an indication for a possible third component (B4), while image C is resolved into four components (C11, C12, C22 and C21). The shift is more visible in region C, which is at higher magnification (μ ∼ 270 to 350). Contours are at (0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.5, 0.75, 1) × the peak flux density of each individual image, which is ∼22 mJy beam−1 for Epoch 1 and 23 mJy beam−1 for Epoch 2. The restoring beam is shown in the bottom left corner and is 11 mas × 5 mas at a position angle of 10°. All of the images are obtained using a Briggs weighting scheme, with Robust = 0.

Open with DEXTER
In the text
thumbnail Fig. 2.

Convergence map for the lens mass model of MG B2016+112. The field-of-view is 6 arcsec × 6 arcsec. The white contours are the 1.7 GHz observations at Epoch 2, while the black contours are the critical and caustic curves. The filled magenta circle indicates the location of the source components. The white labels indicate the groups of lensed images (A, B and C), and the black labels identify the two lensing galaxies (D and G1), also shown by the convergence peaks.

Open with DEXTER
In the text
thumbnail Fig. 3.

Diagonal histograms: probability density functions (PDFs) for the lens model parameters of the main lensing galaxy (D), which were obtained by marginalizing over all the other model parameters, with two dashed vertical lines indicating the 1σ limits. The other panels show the 2-dimensional contours of the PDF for each pair of model parameters, where the contours indicate the 1σ region. The meaning of the parameters, their maximum-likelihood model value and their uncertainties are presented in Table 3.

Open with DEXTER
In the text
thumbnail Fig. 4.

Source-plane model for the VLBA observations of MG B2016+112. Epoch 1 (blue) and 2 (red) observations aligned on the lensing galaxy position, which is at (0, 0). The solid blue line represents the caustics, which divides the source-plane into the double- and quadruple-image regions. The position of the sources is indicated by the filled circles. Source 1 corresponds to images A2–B2, source 2 corresponds to images A1–B1, source 3 corresponds to images C11–C21, and source 4 corresponds to images C12–C22. The labels “core” and “jet” are based on the radio spectral energy distribution of each image pair, as reported by More et al. (2009). The grey scale bar at the bottom left corner represents 100 pc at redshift z = 3.2773. The error bars take into account the uncertainties in the lens model. In the case of sources 1 and 2 (related to the lensed images A and B), the error is dominated by our ability to de-convolve the different source components, while for sources 3 and 4 (related to the lensed region C), the errors are very small because the sub-components are clearly spatially resolved due to their extremely large magnification (μ ∼ 270 to 350).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.