Issue 
A&A
Volume 635, March 2020



Article Number  A86  
Number of page(s)  16  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201936628  
Published online  11 March 2020 
Modelindependent and modelbased local lensing properties of B0128+437 from resolved quasar images
^{1}
Universität Heidelberg, Zentrum für Astronomie, Astron. RechenInstitut, Mönchhofstr. 12–14, 69120 Heidelberg, Germany
email: j.wagner@uniheidelberg.de
^{2}
School of Physics and Astronomy, University of Minnesota, 116 Church Street, Minneapolis, MN 55455, USA
email: llrw@umn.edu
Received:
3
September
2019
Accepted:
13
January
2020
The galaxyscale gravitational lens B0128+437 generates a quadrupoleimage configuration of a background quasar that shows milliarcsecondscale subcomponents in the multiple images observed with VLBI. As this multipleimage configuration including the subcomponents has eluded a simple parametric lensmodel characterisation so far, we determined local lens properties at the positions of the multiple images with our modelindependent approach. Using PixeLens, we also succeeded in setting up a global freeform mass density reconstruction, including all subcomponents as constraints. We compared the modelindependent local lens properties with those obtained by PixeLens and those obtained by the parametric modelling algorithm Lensmodel. A comparison of all three approaches and a modelfree analysis based on the relative polar angles of the multiple images corroborate the hypothesis that elliptically symmetric models are too simplistic to characterise the asymmetric mass density distribution of this lenticular or latetype galaxy. Determining the local lens properties independently of a model, the sparsity and the strong alignment of the subcomponents yield broad 1σ confidence intervals ranging from 8% to over 1000% of the local lens property values. The lens model approaches yield comparably broad confidence intervals. Within these intervals, there is a high degree of agreement between the modelindependent local lens properties of our approach based on the subcomponent positions and the local lens properties obtained by PixeLens. In addition, the modelindependent approach efficiently determines local lens properties on the scale of the quasar subcomponents, which are computationally intensive to obtain by freeform modelbased approaches. Relying on the quadrupole moment of each subcomponent, these smallscale local lens properties show tighter 1σ confidence bounds by at least one order of magnitude on the average with a range of 9% to 535% of the of the local lens property values. As only 40% of the smallscale subcomponent local lens properties overlap within the confidence bounds, mass density gradients on milliarcsecond scales cannot be excluded. Hence, aiming at a global reconstruction of the deflecting mass density distribution, increasingly detailed observations require flexible freeform models that allow for density fluctuations on milliarcsecond scale to replace parametric ones, especially for such lenses as B0128, which have an asymmetric mass density distribution that may include localised inhomogeneities.
Key words: dark matter / gravitational lensing: strong / methods: analytical / galaxies: individual: B0128+437 / galaxies: luminosity function / mass function / quasars: general
© ESO 2020
1. Introduction
B0128+437 is a galaxyscale gravitational lens that has been discussed with some controversy in the literature. It was discovered by Phillips et al. (2000) as a gravitational lensing configuration of four multiple images of a quasar. Three of these images show three subcomponents each in the radio band (Norbury 2002; Biggs et al. 2004). So far, no lens model with a simple, symmetric smooth mass distribution has been found that can explain the positions of all four quasar images and their subcomponents. Using the data from observations described in Biggs et al. (2004), we further investigated this problem with the modelindependent approach as developed in Wagner & Tessore (2018) to determine the ratios of scaled mass densities (convergences) and the reduced shear components at the angular positions of the multiple images. Subsequently, we compared these values to the ratios of scaled mass densities and reduced shear components that the parametric lens model approach Lensmodel (Keeton 2001, 2004) and the freeform lens model approach PixeLens (Saha & Williams 2004) predict. Compared to previous lens models, we took into account the current bestfit redshifts of the lens and the source.
This comparison is analogous to the one carried out in Wagner et al. (2018) on the galaxycluster scale lens CL0024 and shows that the modelindependent approach can be applied to gravitational lensing configurations of any size in the same way. For both lens scales, the ratios of convergences and the reduced shear components of the modelindependent approach show a high degree of agreement to the same local lens properties obtained by lens modelling approaches. In this work, we investigate, in addition, whether the modelindependent approach is able to further constrain the local lens properties on the level of the subcomponents in the multiple images of the quasar behind B0128. Such lens reconstructions as PixeLens describe nonsmooth irregularly shaped mass density distributions on a pixelised grid. Therefore, such methods usually require computationally intensive pixelisations to determine these smallscale properties. For lens reconstructions based on analytical models like Lensmodel, solving the lens equation as an optimisation problem may also require a highly resolved sampling grid to determine smallscale details of the lens. Thus, a more efficient way to obtain smallscale local lens properties is highly desired for the increasing amount of data obtained by upcoming surveys.
The paper is organised as follows: Section 2 introduces information about B0128 that has become available so far. It is mainly based on the works of Phillips et al. (2000), Norbury (2002), Biggs et al. (2004), and Lagattuta et al. (2010). Subsequently, we list all observational data that we employed to calculate the modelindependent, local lens properties and that were used to constrain our lens models in Sect. 3. In Sect. 4, we briefly outline the modelindependent algorithm, which is further detailed in Wagner & Tessore (2018) before we show the modelindependent, local lens characterisation of B0128. Analogously, Sects. 5 and 6 describe the lens modelling based on the parametric code Lensmodel, and the freeform code PixeLens, respectively, before applying it to B0128. To avoid any potential confirmation bias due to exchanging the values to be compared at an early stage, the evaluation is blinded in the same way as the evaluation performed in Wagner et al. (2018). This means that the values of the local lens properties are determined independently for the modelindependent and the modelbased approaches and only revealed for the comparison at the very end, after all modelling is complete. In Sect. 7, we apply the modelfree comparison as established in Woldesenbet & Williams (2012, 2015) and Gomer & Williams (2018) to B0128 to investigate its degree of asymmetry with respect to lenses with double mirror symmetry based on the relative polar angles between the multiple image positions. In Sect. 8, we compare all results obtained in Sects. 4–7. Finally, we conclude by assembling a consistent picture of the lensing configuration in B0128 and summarising the methodological results that can be deduced from the comparison of the modelindependent and the modelbased lens reconstructions.
2. Related work on B0128
In the discovery paper, Phillips et al. (2000), highresolution MERLIN observations at 5 GHz from B0128 were obtained in the course of the Cosmic Lens AllSky Survey (CLASS). Four unresolved multiple images, meaning those without visible subcomponents down to the scale of 30 mas, with a maximum image separation of 0.54 arcsec, arranged in a classic quadlens formation were detected (see Fig. 1 (right)). Longterm observations showed no hint of a timevariability of the quasar.
Fig. 1. Left: HST Iband observation (HST Proposal 9133, based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (STECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA)) with B0128 (red circle) and the shear galaxy (yellow circle); the image detail shows the cleaned Iband image of B0128 as obtained by the CASTLe Survey, https://www.cfa.harvard.edu/castles/. Right: MERLIN 5 GHz observation by Phillips et al. (2000) showing the four multiple images; map peak: 0.0174 Jy beam^{−1}, contours: 0.000277 Jy beam^{−1} × (−3 3 6 12 24 48) amounting to multiples of 48 times the rootmeansquare (rms) noise of 277 μJy beam^{−1}, beam fullwidth at half maximum (FWHM): 58.8 × 47.7 mas^{2} at 50.9°. 
B0128 was modelled by a singular isothermal ellipse (SIE), using the angular positions and flux ratios of the four (unresolved) multiple images to constrain the lens model. Redshifts of z_{l} = 0.5 for the lens and z_{s} = 1.5 for the source were assumed in a flat universe with Ω_{m0} = 1 as today’s matter density parameter. High deviations between the observed and the modelpredicted image positions were found. Including external shear alleviated the discrepancies, which indicated that the deflecting mass distribution might be less smooth or symmetric than previously assumed.
To investigate this phenomenon further, optical followup observations were performed as detailed in Norbury (2002) (see Fig. 1 (left) for an HST observation showing B0128 and its environment, including a galaxy that could be responsible for the external shear that was introduced in the SIElens model in Phillips et al. 2000). They revealed that all signals from B0128 are highly reddened, such that the source is most probably at a high redshift or the lens contains a large amount of dust. In addition, followup VLBA 5 GHz data were obtained, showing that three of the four images (images A, C, and D in Fig. 1 (left)) could be further decomposed into three subcomponents.
Subsequently, Biggs et al. (2004) summarised the results found in Norbury (2002) to conclude that image B in Fig. 1 (right) is most likely scatterbroadened and hence, it is difficult to decompose it into subcomponents. They also performed VLBIobservations at 2.3, 5, and 8.4 GHz to obtain further details about the subcomponents and found that the observations could rather show a core and jet instead of a compact symmetric object, as originally assumed in Phillips et al. (2000). The detailed structure of all four images with the labelling of subcomponents according to Norbury (2002) and Biggs et al. (2004) is shown in Fig. 2. As with the lens models by Norbury (2002), the lens models by Biggs et al. (2004), based on the relative image positions as constraints, were not able to fully describe the lensing configuration to submilliarcsecond precision, even when different algorithms and modelling principles were employed.
Fig. 2. VLBA 8.4 GHz details of all four multiple images from Biggs et al. (2004), map peaks for images A, C, and D: around 1.4 mJy beam^{−1}, for image B: 0.39 mJy beam^{−1}, contours: 0.000225 Jy beam^{−1} × (−1 1 2 4 8 16 32 64 128 256 512), amounting to multiples of 3σ where σ is the offsource rms noise in the map (90 μJy beam^{−1}), beam FWHM: 1.2 × 0.7 mas^{2} at −21.6°. 
McKean et al. (2004) finally determined the redshift of the source to be z_{s} = 3.1240 ± 0.0042 in a KECK observation and found another emission line, which, depending on its origin, allows the lens to be at redshifts z_{l} = 1.145, 0.645, or 0.218. The assumption that the lens redshift is either z_{l} = 0.645 or 1.145 was corroborated by further KECK observations as detailed in Lagattuta et al. (2010). This work also corroborated all previous findings and set up the hypothesis that B0128 could be a lenticular or a latetype galaxy.
Later, Sluse et al. (2012) modelled B0128 as an SIE with external shear as Phillips et al. (2000), but used Ω_{m0} = 0.3 and Ω_{Λ} = 0.7 as today’s cosmological parameters, in accordance with recent observations, requiring the cosmological model to introduce Ω_{Λ} ≠ 0. They assumed the lens at a redshift of z_{l} = 0.6, which is less likely than z_{l} = 1.145, according to Lagattuta et al. (2010). Sluse et al. (2012) conclude that the model fits the observation until submas precision for the image positions is required. To investigate the most probable source of the discrepancies that enter at submas precision, Xu et al. (2015) determined the probability that the flux anomalies between the images A and B are caused by yet unobserved dark matter substructures in the lens and concluded that the oversimplified or improper lens model in addition to the scatterbroadening is more likely to cause the flux anomaly than additional smallscale darkmatter halos.
Summarising the previous results, observations indicate that
– B0128 is a potentially lenticular or latetype galaxy most likely located at redshift z_{l} = 1.145 that generates four multiple images of a quasar at redshift z_{s} = 3.124,
– three of these multiple images, A, C, and D in Fig. 1 (right), can be resolved into three clearly identifiable subcomponents, the fourth, image B in Fig. 1 (right), is most likely scatterbroadened,
– no simple, highly symmetric lens model, even including the shear galaxy located 7.8 arcsec away, which is shown in Fig. 1 (left), can fully explain the positions of all subcomponents to submasprecision simultaneously.
3. Observed data from the multiple images
For our analysis, we used the data, as already stated in Phillips et al. (2000), Biggs et al. (2004), and Lagattuta et al. (2010). Tables 1 and 2 summarise the image positions and flux densities for the unresolved and resolved observations. Since it is very difficult to determine the relative angular positions for the hardly visible subcomponents in image B, we assumed an uncertainty in these offsets of 1 mas and also investigated the impact, if we increased the uncertainty to 3 mas. A similar approach was pursued in Biggs et al. (2004).
MERLIN 5 GHz measurements as in Table 2 of Phillips et al. (2000) (Cols. 2–4), VLA 8.4 GHz measurements as in Table 1 of Phillips et al. (2000) (Cols. 5–7), and as in Table 2 of Lagattuta et al. (2010) using the Kpfilter of the Near Infrared Camera 2 on the Keck II telescope along with the LGS AO system (Cols. 8–10), both showing unresolved multiple images.
VLBI observations as in Table 1 from Biggs et al. (2004).
In order to align the coordinates of the KECK observations in the optical from Lagattuta et al. (2010) with the MERLIN observations from Phillips et al. (2000), we had to identify a reference point in both observations. To do so, we chose the angular position of image C, since the resolved components lie closest together and, as the counter image opposite to the threeimage configuration B, A, D, it is farther from the critical curve as those images and is therefore subject to the least amount of magnification.
The spectroscopic observation of McKean et al. (2004) did not find any hints of a second source and the similarity of the subcomponents in the images A, C, and D is high, so that we treated the images A, B, C, and D as multiple images from the same source at z_{s} = 3.124. We assumed that the images B, A, and D form a cusp configuration and that image C is of the same parity as image A, see Wagner & Tessore (2018) for further details about fold and cusp configurations and their parity.
Mapping the subcomponents of the images onto each other should then work analogously to the mapping of reference points in CL0024, as done in Wagner et al. (2018). The subcomponents in the images in B0128 are much more aligned, which causes higher uncertainties. In addition, three subcomponents are the minimum number of reference points required in each image to apply the method outlined in Wagner & Tessore (2018) and Wagner et al. (2018). Considering images A and B as a fold configuration and images A and D as another fold configuration of two images mirrorinverted at the critical curve between them, we noticed that the matching of subcomponents between A and D as shown in Fig. 2 did not yield a fold configuration on the scale of the subcomponents. To be a fold configuration, the order of subcomponents in image D should be reversed with D_{1} in the north and D_{3} in the south. Matching the subcomponents by their intensity led to this labelling. Yet, the multiband observations carried out in Norbury (2002) and Biggs et al. (2004) strongly hint at dust in the lens, so that image B might not be the only image that is subject to attentuation due to scatterbroadening. As a consequence, the matching of the subcomponents of image A and D could be different than proposed in Biggs et al. (2004). In Sect. 4, we systematically investigated the impact of different matchings on the local modelindependent properties and the lens model.
Table 3 shows the flux density ratios of images B, C, and D with respect to image A for all available bands, ℱ_{i}, i = B, C, D. Biggs et al. (2004) assumed less than 5% measurement uncertainties on their flux ratios. Lagattuta et al. (2010) obtained Δℱ_{B} = 3%, Δℱ_{C} = 1%, and Δℱ_{D} = 3%. These flux density ratios were compared to the magnification ratios determined by the modelindependent approach and the lens model as a consistency check because we did not employ them as constraints.
Observed flux density ratios ℱ_{i} with respect to image A for images i = B, C, D in different bands.
Natural weighting yields a high signal to noise ratio with a low angular resolution, such that it focuses on smooth, extended subcomponent brightness structures. Comparing the first two rows in Table 3, we observed that the flux ratios are higher for this weighting scheme than for the uniform weighting which has a lower signal to noise ratio but a higher angular resolution. Hence, images B, C, and D must be better visible, that is, have a higher (flux) density, when the natural weighting scheme is applied. Furthermore, the fact that ℱ_{B} is much smaller for the uniform than for the natural weighting scheme, strongly hints at scatterbroadening due to the dust in the lens, as already noted by Biggs et al. (2004).
4. Modelindependent reconstruction
4.1. The method
Effectively describing the gravitational lens as a twodimensional mass distribution in a single lens plane at redshift z_{l}, the standard gravitational lensing formalism (see, for instance, Petters et al. 2001 or Schneider et al. 1992 for an introduction) treats this projected deflecting mass distribution in terms of the convergence κ(x). This is the twodimensional mass density at the position x ∈ ℝ^{2} in the lens plane scaled by the critical density Σ_{cr} which is the sufficient mass density to generate multiple images, given a lens at z_{l} and a source at z_{s}. While κ(x) only enlarges or diminishes the source, the shear γ(x) = (γ_{1}(x),γ_{2}(x)) additionally distorts the images. Both κ(x) and γ(x) can be determined as secondorder derivatives from the twodimensional gravitational deflection potential ψ(x) as detailed in Wagner & Tessore (2018).
Usually, global lens reconstructions, as detailed in Sects. 5 and 6 are set up, using the observables from the multiple images as constraints to reconstruct the deflecting mass density distribution as a whole. As observables, relative image positions, the quadrupole moment of the images around their centre of light, the flux ratios, and the time delays, if available, are employed. These lens models are subject to a lot of degeneracies, see Wagner (2018a) for a mathematical derivation and Wagner (2019a) for the physical explanation of all degeneracies arising in the general lensing formalism. To avoid the modelbased degeneracies, Tessore (2017) derived the most general information which can be obtained from multiple images of a background source without assuming a specific gravitational lens model. He found that transforming the multiple images onto each other yields ratios of convergences,
between all multiple images i, j and reduced shears,
at the positions of the multiple images. These local lens properties in Eqs. (1) and (2) are invariant under the masssheet transformation. Thus, as further elaborated in Wagner (2019b), they represent the information about the lens that all lens models should have in common at leading order. Derived from the f_{ij}s and g_{i}s, the magnification ratios,
for image pairs i, j can be calculated because the magnification μ_{i} of an image i is given as the inverse of the determinant of the distortion matrix
As minimum requirements, the brightness profiles of the multiple images must contain clearly identifiable subcomponents, for instance, at least three linearly independent reference points like star forming regions, that can be matched across all multiple images. Furthermore, at least three multiple images with these identifiable subcomponents are needed which are not aligned like in an axisymmetric deflection potential.
Our Cimplementation^{1} employs the centroid of all reference points as the position of a multiple image, also called the “anchor point”. Apart from this choice, any point in the convex hull of the reference points can be used as anchor point, since the approach assumes that the convergence and the shear are approximately constant over the extension of a multiple image. The position of the anchor point of one image, the socalled “reference image” remains fixed. Then, a system of equations is set up by linearly mapping the reference image onto all other multiple images. Solving the system of equations by a χ^{2}parameter estimation as detailed in Wagner et al. (2018) yields all f_{ij} and g_{i} for all multiple images, and the most likely anchor point positions in all remaining images. Wagner & Tessore (2018) and Wagner et al. (2018) further detailed the algorithmic implementation and the procedure to obtain confidence bounds on the local lens properties by sampling their covariance matrix close to the most probable values of the f_{ij}, g_{i}, and the anchor points.
4.2. Results for B0128
Next, we applied the method as outlined in Sect. 4.1 to the case of B0128. We first noted that the four nonaligned multiple images with their three subcomponents fulfil the requirements stated in Sect. 4.1. The subcomponents are almost linearly aligned, yet, not in a mathematically rigorous way to cause the optimisation problem which determines the local lens properties to be underconstrained and thus degenerate. But, from the systematic analyses performed in Wagner et al. (2018), we expected broad confidence bounds for the local lens properties determined by the subcomponents due to their alignment and the small area that they span.
We chose image A as the reference image and set up transformations between all remaining images to image A in order to determine the local lens properties. To investigate the impact of choosing image A as reference image, we also determined the local lens properties in the threeimage configuration of images A, C, and D mapping the three identifiable subcomponents onto each other and using image C or image D as reference image. In all three cases, we obtained the same local lens properties, yet with different confidence bounds. The analysis shows that using image C as reference image yields smaller confidence bounds for the threeimage configuration ACD. Yet, it yields comparable to slightly worse confidence bounds when matching the individual Gaussian fits of the subcomponents in image A, C, and D onto each other and we obtained larger confidence bounds for the fourimage configuration ABCD. Thus, image C is only slightly less suitable as reference image than image A. Contrary to that, using image D as reference image yields smaller confidence bounds for the threeimage configuration, yet with a very low effective number of samples in the importance sampling process.
To simplify the notation, we omit A in the subscripts of the local lens properties and simply write
Instead of using the positions of the three subcomponents in each multiple image as reference points, we reconstructed local lens properties on the scale of the subcomponents themselves. For this, we chose image A again as reference image and used the centre of light of the Gaussian fitted to the subcomponent i as first reference point. Then, we used the end points of the semimajor and semiminor axis of the fitted Gaussian as further reference points in each subcomponent. The calculations how to obtain these positions from the values of the semimajor axis, the axis ratio, and the position angle can be found in Appendix A. In this way, we matched the subcomponents 1 and 3 across images A, C, and D to infer the individual local lens properties at their positions. Subsequently, we compared these local lens properties obtained for the subcomponents with the ones obtained on image scale. If reduced shear and convergence are constant over the area of a multiple image, the local lens properties of the subcomponents should coincide with each other and with the imagescale local lens properties. The matching for the image and subcomponentscale is sketched in Fig. 3.
Fig. 3. Visualisation of imagescale matching: using the positions of the three subcomponents (red, blue, green dots) in each image, the multiple images A, B, C, and D are mapped onto each other (red arrows, indicating that the dots of each colour are mapped onto each other). This mapping determines the local lens properties on the image scale (indicated by the larger red areas for visualisation purposes, as the convex hull spanned by the reference points in each image is too small to be drawn). Visualisation of subcomponentscale matching (highlighted in blue): using the Gaussians fitted to a subcomponent (here: for the first subcomponent denoted by 1, highlighted in dark grey), the subcomponents can be matched onto each other (as indicated by mapping the dots of each colour onto each other at the subcomponent scale) to determine the local lens properties within the area of the Gaussians. The first reference point in each Gaussian is its centre of light (yellow dots). The two reference points that are further required are the end points of the axes of the elliptical isocontour (brown, white dots). 
In the following, we start with determining the imagescale local lens properties by matching the subcomponents for images A, C, and D onto each other (see Sect. 4.2.1). To investigate the impact of a relabelling of the subcomponents, we systematically interchanged the 1 and 3 labels for images C and D. In this way, the choice of subcomponent labelling according to Biggs et al. (2004) was crosschecked. Then, in Sect. 4.2.2, we included image B and systematically interchanged its subcomponents 1 and 3 to determine the most likely local lens properties for all four images. As image B is subject to scatter broadening and its subcomponents are hard to detect, a comparison between the local lens properties obtained from matching images A, C, and D onto each other with the local lens properties obtained by mapping all images onto each other, reveals how reliably the subcomponents of image B were identified. In Sect. 4.2.3, we determine the local lens properties at subcomponent scale. We matched the Gaussians fitted to subcomponents 1 and 3 across the images A, C, and D to investigate their relationship with respect to each other and to the imagescale local lens properties. In addition, we discuss potential biases due to dust in the lens plane, as assumed in Norbury (2002), Biggs et al. (2004), and Lagattuta et al. (2010).
4.2.1. Imagescale matching of images A, C, and D
Employing the positions of the three subcomponents 1, 2, and 3 in images A, C, and D listed in Table 2 as reference points, we determined the local lens properties on image scale. The most likely lens properties, the mean values and the 1σ confidence bounds are listed in the second, third, and fourth column of Table 5, respectively. These results were also plotted as the set of first four points in Fig. 4.
Fig. 4. Graphical representation of the results presented in Tables 5 and 7. Each panel corresponds to one image of B0128, as labelled. Within each panel, we plot the values of g_{1}, g_{2}, f, and 𝒥 in 4 columns, separted by thin black vertical dotted lines. In each of these columns, there are 10 points, with confidence bounds, representing the four models from Table 5, and six models from Table 7 For image A, all values of f and 𝒥 are, by definition, 1. 
As mentioned in Sect. 3, the labelling of the subcomponents according to Table 2 is not compatible with image A and D being a fold configuration on the level of subcomponents. To systematically investigate which subcomponents across the images should be matched, we interchanged the labelling of the subcomponents 1 and 3 systematically as indicated in Table 4. The resulting lens properties can be found in Appendix B. We found that only the labelling according to Biggs et al. (2004) yields the correct relative parities between the images A, C, and D in the signs of the relative magnifications 𝒥_{i}, i = C, D. None of the configurations with other labelling are able to obtain the correct relative parities between images.
Configurations of differently matched subcomponents across images A, C, and D.
Thus, the labelling of subcomponents according to Biggs et al. (2004) is the only one yielding results which are consistent with leading order lensing theory. Considering the ratios of convergences f_{i}, we noted that the highly negative value for f_{D} could indicate that images A and D lie on opposite sides of the isocontour κ(x) = 1. The positive sign for f_{C} could indicate, that images A and C are located at the same side of κ(x) = 1. Since f_{C} and f_{D} are both subject to broad confidence bounds several times their absolute value, these statements are further investigated in Sects. 4.2.2 and 4.2.3. In Sect. 4.2.4, they are also checked for consistency in a comparison with our lens models as set up in Sects. 5 and 6.
4.2.2. Imagescale matching of all images
Building upon the configuration of three multiple images A, C, and D with the three subcomponent positions as shown in Table 2, we included the three subcomponent positions in image B and determined the local lens properties from it. The results can be read off Cols. 5–7 in Table 5 (also plotted in Fig. 4) and are mostly subject to confidence bounds that exceed their absolute values. These large confidence bounds can partly be caused by the scatterbroadening of image B, which is not accounted for in our local lens reconstruction.
Swapping the labelling of subcomponent B_{1} with B_{3}, we found that both results agree within their confidence bounds except for g_{A}. We favoured the labelling as proposed by Biggs et al. (2004) due to the slightly lower confidence bounds of the resulting local lens properties and because the most likely value for 𝒥_{B} has the correct parity, which is not the case when interchanging the labels.
In the fourimage configuration, images A, C, and D most likely lie on the same side of the isocontour κ(x) = 1 and image B seems to be located on the opposite side. Yet, analogously to the previous results, f_{B} has a large confidence bound, which also includes the possibility for it to lie on the same side as image A.
Increasing the uncertainty in the positions of the subcomponents of image B from 1 mas to 3 mas leads to an increase in the confidence bounds. The local lens properties still agree within their confidence bounds except for g_{A}. Assuming an uncertainty of 3 mas in the position of all subcomponents, the confidence bounds do not increase, yet, the effective number of samples in the importance sampling is reduced below 10. In this case, the local lens properties except for g_{A, 1} coincide with the ones using 1 mas as uncertainty in the positions of the B_{i}, i = 1, 2, 3. We show all results for the interchange of B_{1} and B_{3} and the increase of the measurement uncertainties in Appendix C.
4.2.3. Subcomponentscale matching of subcomponents 1 and 3
After the imagescale local lens properties were obtained, we turned to the subcomponentscale local lens properties. As the Gaussians that are fitted to the subcomponents are elliptically symmetric, several choices of reference points are available for each subcomponent, apart from the ones shown in Fig. 3. All these choices of reference point configurations should be equally well suited for determining the local lens properties on subcomponent level. Therefore, as a consistency check, we first investigated whether the different configurations of reference points taken from the end points of the semimajor and semiminor axes yield the same local lens properties. As detailed in Appendix A, this is the case for both subcomponents 1 and 3. Subsequently, we determined the local lens properties at the positions of subcomponent 1 in images A, C, and D as listed in Cols. 8–10 of Table 5 and the ones at the positions of subcomponent 3 as listed in Cols. 11–13 of Table 5. These values are also plotted in Fig. 4.
Contrary to the imagescale local lens properties, all f and gvalues have confidence bounds that are smaller than their absolute value except for g_{C, 1} and g_{D, 1} of subcomponent 3. We found that half of the lens properties at the position of subcomponent 1 agree with the ones at the position of subcomponent 3 within their confidence bounds.
Concerning their relative positions with respect to the isocontour κ(x) = 1, the most likely f_{i}, i = C, D in Table 5 imply that all images are supposed to lie on the same side to a much higher degree of confidence because the confidence bounds do not include negative values anymore.
4.2.4. Synopsis and comparison to previous results
From Sect. 4.2.1, we drew the conclusion that it is possible to determine local f_{i} and g_{i}, i = A, C, D, such that their values are constant over the area spanned by the three subcomponents in each image. The matching of the subcomponents according to Table 2 is the only one yielding the correct relative parities in the relative magnifications 𝒥_{C} and 𝒥_{D}. These findings are in agreement with previous results as summarised in Sect. 2.
Including the subcomponent positions of image B, we found that the local lens properties for images C and D agree within their confidence bounds and all relative confidence bounds do not significantly increase. These results indicate that it should be feasible to construct a smooth lens model that explains the fourimage configuration with its subcomponents. So far, only the specific approaches based on highly symmetric lens models detailed in Biggs et al. (2004) have been unsuccessful.
In addition, we noted that we obtained similarly large relative confidence bounds for the case of the galaxyclusterscale lens CL0024 (Wagner et al. 2018), when we reduced the number of reference points to four that are almost linearly aligned and span only a small area. Thus, apart from the scatterbroadening due to the relatively large amount of dust in the lens, the sparsity of the data and their alignment also contributes significantly to the large confidence bounds. The smaller confidence bounds determined for the local lens properties of the subcomponents, due to the orthogonally oriented vectors between the reference points, support this hypothesis. Without further, improved multiband observations, it is hard to disentangle the impact of the individual effects.
Comparing the imagescale local lens properties as obtained in Sect. 4.2.1 to the local lens properties at the positions of the subcomponents 1 and 3, as derived in Sect. 4.2.3, we found that the local lens properties at the positions of the two subcomponents agree with the ones determined over their entire convex hull with the exceptions of g_{A} at both positions and 𝒥_{C} and 𝒥_{D} at the position of subcomponent 1. Comparing the local lens properties at the positions of the two subcomponents with each other, we observed that half of them agree within their confidence bounds. This puts a very weak upper limit to the scale of potential higher order perturbations, like gradients, in the surface mass density at the position of the images A, C, and D.
Local lens properties obtained at the positions of the subcomponents 1 and 3 indicate that all images lie on the same side of the isocontour κ(x) = 1. Most probably, they lie outside κ(x) = 1, given that κ(x) = 1/2 for an SIE and that such a smooth, symmetric lens model can be fitted to the fourimage configurations, if no submas precision of the positions of the subcomponents is assumed.
5. Modelbased reconstruction: parametric
5.1. The method
The most common way to model galaxy lenses is to fit image positions using a simple parametric form for the galaxy mass distribution. We used publicly available software, Lensmodel, described in Keeton (2001, 2004). Lensmodel offers a range of lens models; we chose to work with analytic potentials, instead of mass distributions, because the former have analytic expressions for all relevant lensing quantities. The reconstructions presented in this section use a very different ansatz from the ones in Sect. 6, in that the latter recover a mass distribution, instead of a potential, and do not assume any fixed parameteric form. The alphapot lensing potential is a softened power law potential given by
where b is the normalisation constant, s is the core radius, which is set to a small nonzero value, and , with q being the axis ratio of the potential. The boxy power law potential, called boxypot, expressed in polar coordinates x = (r, θ) has the form
where ϵ is the ellipticity and θ_{ϵ} is its position angle^{2}. In the reconstructions, both potentials are augmented with external shear.
5.2. Results for B0128
Employing the Lensmodel lens models as outlined in Sect. 5.1, we reconstructed the parametric mass density distribution on imagescale, with the positions of the subcomponents of the multiple images as constraints. We ran several lens reconstructions for each of the two image configurations (abbreviated as config ABCD and config ACD; we used the naming convention as introduced in Table 5), and for the cases with and without fixing the galaxy lens centre. The reconstructions use different initial parameter guesses, and alphapot and boxypot potentials.
If only subcomponent 1 was used as input, config A_{1}C_{1}D_{1} and config A_{1}B_{1}C_{1}D_{1} could be fit perfectly with either alpha pot or boxypot, if the lens centre coordinates was left as free model parameters. This is summarised in the second column of the first two rows of Table 6, which shows typical rms deviations between the observed image positions and the modelpredicted ones in the lens plane. If the lens centre was held fixed at its observed position, config A_{1}C_{1}D_{1} could still reproduce the images perfectly (third column of the first row). If config A_{1}B_{1}C_{1}D_{1} was used, typical rms became 0.022″. When the other two subcomponents, 2 and 3, were included, for a total of nine and twelve subcomponents for each of config ACD and config ABCD (last two rows), neither alphapot nor boxypot provided good fits, for floating or fixed lens centre. These lens plane rm are larger than our assumed uncertainty of 0.1 mas for images A, C, and D.
Lensmodel results: typical lens plane image rms.
The resulting Lensmodel local lens properties for config A_{1}B_{1}C_{1}D_{1} and config A_{1}C_{1}D_{1}, as defined in Eqs. (1) and (2), were extracted at the positions of subcomponent 1 and listed in the Cols. 4–7 of Table 7, and summarised in Fig. 4. The tightest confidence bounds were found for config A_{1}B_{1}C_{1}D_{1} with fixed lens centre. This is probably because, in this case, the modelling has the least freedom to chose the lens mass distribution.
It is interesting to note that in B0128, ruling out simple lens models is possible because of excellent astromentric precision provided by the radio data (0.01 milliarcsecond uncertainty), and further, by the presence of subcomponents in the multiple images on scales smaller than 0.01″, that is, substantially smaller than HST pixel size. If the source in B0128 were extended, with lensed images covering one or more HSTsized pixels in the lens plane, there would be no indication that simple models do not fit.
6. Modelbased reconstruction: freeform
6.1. The method
Because simple parametric models, like elliptical mass distributions with external shear presented in Sect. 5, cannot reproduce all twelve subcomponents in B0128, here, we used a freeform method, called PixeLens (Saha & Williams 2004), to reconstruct the mass density distribution in B0128. PixeLens is publicly available, and has an easytouse GUI interface^{3}.
For any given basis set, the lens equation can be written as a set of linear equations in the unknowns, which are the weights of the basis functions, and the source positions. PixeLens breaks up the lens plane into equal size square mass pixels (its basis set), and imposes a few constraints. The positions of the images are specified with respect to the centre of the lens, which serves as the centre of the reconstruction. The parities of the images are also specified as input and are strictly enforced. The mass gradient should point not more than ±45° away from purely radial. Except for the central pixel, the mass of no other pixel can exceed twice the average mass of all its neighbours. These constraints act to regularise the mass distribution. Because many sets of pixel weights, that is, many mass distributions, can reproduce the images exactly, there is a plethora of solutions. In this paper we ran PixeLens for 250 models, and discarded the first fifty “burnin” models, as is sometimes done for MCMC runs. The final PixeLens lens reconstruction was then taken as an average over 200 individual solutions.
The confidence bounds for the local lens properties were calculated as the rms dispersion, around the mean, between twenty sets of ten individual models. We used ten, instead of one, because, averaging over a handful of individual mass density models suppresses oneoff astrophysically unrealistic features and enhances common features. We experimented with averaging ten, twenty, and forty individual models to obtain one final PixeLens model, with the corresponding number of twenty, ten, and five sets of PixeLens models to determine the rms. As expected, we found that, the rms between the sets decreases. In the limit, if averaging is done over a very large number of individual reconstructions, the difference between sets go to zero, and so does the rms. Our choice of averaging over ten individual models to obtain one final model is somewhat conservative, because the corresponding rms is on the high side.
6.2. Results for B0128
First, we carried out PixeLens mass density reconstructions using all twelve subcomponents of all four images in B0128. Relative images fluxes were not used. A region of radius 0.675 arcsec around the central reference point in Table 2 was divided into 41 by 41 pixels, such that each PixeLens pixel covers an area with an edge length of 33 milliarcseconds. The average projected lensing convergence map is shown in the left panel of Fig. 5. It is the average over 200 PixeLens solutions. The recovered mass distribution is not very symmetric, and would be hard to represent with a simple parametric model. This is not surprising, and is consistent with parametric models not being able to reproduce all twelve subcomponents. Put differently, to reproduce all twelve subcomponents, one requires significant deviations from a purely elliptical projected mass distribution.
Fig. 5. PixeLens reconstructed maps of lensing convergence. Left: using all twelve subcomponents of B0128. Right: using only the nine subcomponents of images A, C, and D. The subcomponents used as constraints in each case are indicated with magenta circles. The grey isoconvergence contours are spaced logarithmically. The blue contour is the isoconvergence contour at κ(x) = 1. 
Mindful of the fact that mass reconstructions depend critically on the quality of the image data, we also carried out a reconstruction that did not include any subcomponents of the mostlikely scatterbroadened image B. The reconstruction based on just the three subcomponents of each of images A, C, and D is shown in the right panel of Fig. 5. The local lens properties of config ABCD and config ACD as set up in Eqs. (1) and (2) at the pixels, sized 33 milliarcseconds, that contain subcomponent 1, are displayed in the second and third column of Table 7 (see also Fig. 4).
To investigate whether the differences in the two mass maps in Fig. 5 are significant, we calculated the statistical significance at each pixel location in the lens plane as
where the σ’s are the rms obtained using twenty sets of ten individual PixeLens reconstructions. In fact, none of the differences are statistically significant; the value of S is less than 1 for every pixel. Consistently, the f and gvalues for the two models (see the first two models in Table 7 and Fig. 4) all agree within their 1σ confidence bounds. Using alternative labelling of subcomponents, for example Conf. 1 in Table 4, yielded arrival time surfaces with very contorted contours, which is consistent with the results of Sect. 4.2.1.
7. Modelfree analysis
7.1. The method
Finally, we performed another type of analysis on the three sets of subcomponents in images A, B, C, and D of B0128. This is not a mass reconstruction, and so does not yield values of surface mass density or shear. The analysis was described in Woldesenbet & Williams (2012, 2015) and Gomer & Williams (2018). It is based solely on the relative image polar angles of quadrupolyimaged quasars (quads), as viewed from the centre of the galaxy lens.
In a quad, there are three relative angles between images. If one labels each of the four images by their arrival time, then the three relative angles, θ_{ij}, between the ith and jth arriving images, are θ_{12}, θ_{23}, and θ_{34}. In principle, other combinations i and j can be used, but this is the set that has proven to be useful in previous works. Woldesenbet & Williams (2012) analysed quads arising from double mirror symmetric lenses, that is, those whose mass and lensing potential distributions have two mutually orthogonal symmetry axes, like elliptical lenses, or circular lenses with external shear. These authors showed that regardless of the ellipticity or density profile slope of the lens, all quads from such lenses lie on a nearly invariant surface in the 3D space of the three relative image angles, called the Fundamental Surface of Quads (FSQ).
Instead of using the 3D representation, it is easier to use a 2D projection, and plot θ_{23} versus the distance of quads from the FSQ, Δθ_{23}, which is measured parallel to the θ_{23} axis. Note that θ_{23} is singled out because second and third arriving images are the ones that approach each other in the lens plane and vanish when the source moves further away from the lens centre and a quad becomes a double. Therefore, these two images distinguish a quad from a double.
7.2. Results for B0128
Figure 6 is the 2D projection introduced above, and shows the three quads formed by the three sets of subcomponents in B0128 as green squares, together with forty galaxyscale quads (black circles) presented in Woldesenbet & Williams (2012).
Fig. 6. Distribution of quads in the space of relative image angles. Here, a 2D projection of that 3D space is used. The Fundamental Surface of Quads is the horizontal line at Δθ_{23} = 0. Galaxyscale quads from Woldesenbet & Williams (2012) are shown as black circles with errorbars. The three quads of B0128 are represented by green squares, labeled with the corresponding subcomponent number. The orange points are quads from a synthetic lens, plotted here for reference, with projected density profile ∝r^{−0.8}, ellipticity ϵ = 0.25, and external shear γ = 0.25, which is misaligned with the ellipticity position angle by 80°. 
As a reference, we plotted a few thousand quads from a synthetic lens with projected mass density profile ∝ r^{−0.8}, ellipticity ϵ = 0.25, and external shear γ = 0.25, which is misaligned with the ellipticity position angle by 80°. This lens is not double mirror symmetric because the mass ellipticity axis and the external shear are not orthogonal to each other. Deviations from double mirror symmetry can be of several types: one can add a misaligned external shear to an elliptical lens, as we did here, or one can add nonelliptical mass density perturbations to an otherwise elliptical lens. Examples of such deviations are the boxy and disky perturbations observed in the isophotes of some galaxies. Presence of satellite galaxies, or offsets between the dark matter and baryonic components of the lens mass distribution also result in deviations from ellipticity. Strictly speaking, randomly distributed small scale ΛCDM substructure also break the double mirror symmetry. However, the corresponding perturbations to the projected mass are small and localised, and usually not large enough to affect positions of images.
In general, the more a given lens deviates from being purely double mirror symmetric, the more its quads deviate from the FSQ. By comparing the location of the B0128 quads to the quads of the synthetic lens (orange points), we estimated that, if B0128 is fitted with a simple lens model, its ellipticity or shear is approximately 0.25. This conclusion is based on a number of other mock galaxies, with various ellipticities, density profile slopes and external shears, which we do not show here. This approximate value of ellipticity or shear is consistent with the findings of Biggs et al. (2004), whose two models had γ = 0.26 and 0.22, as well as our own findings using Lensmodel (Sect. 5.2), where ellipticity and external shear had magnitudes similar to these.
Instead of requiring the galaxylens to have very high shear, an alternative interpretation to explain why the three quads of B0128 are located so far from FSQ in Fig. 6 is that the galaxy lens has nonelliptical density perturbations, such as the ones described above, and shown in Fig. 14 of Gomer & Williams (2018). Those authors explored the role of such perturbations on the relative image angles of galaxyscale observed quads. They concluded that observed deviations from FSQ by Δθ_{23} ∼ 5 ° −10° are possible if perturbations of the density profile away from a purely elliptical model are included in the mass model.
The possible presence of nonnegligible perturbations from ellipticity in the case of B0128, indicated by Fig. 6 is consistent with the mass distribution produced by freeform PixeLens, in Sect. 6, and shown in Fig. 5.
8. Conclusion
We presented four different types of analyses and reconstructions of the galaxyscale lenticular or latetype gravitational lens B0128 constrained by the quadrupoleimage configuration of a background quasar. Previous multiband observations revealed that each of the four quasar images shows three bright subcomponent features. Radio observations resolved these subcomponents, which are separated by less than ten milliarcseconds, giving constraints on the lensing mass distribution on very small scales. All approaches to find a global mass density reconstruction of B0128 based on lens models to reproduce the subcomponent configurations within the multiple images, mainly pursued by Biggs et al. (2004), have been unsuccessful. In contrast to that, the four multiple images observed at a lower resolution at which no subcomponents are resolved can be reproduced by a singular isothermal ellipse lens plus external shear.
Section 8.1 summarises the methodological progress that could be made to explain the multipleimage configuration at subcomponentscale by analysing B0128 with different lens characterisation approaches and comparing their results with each other. Subsequently, we conclude in Sect. 8.2 by summarising the consistent lens description that can be set up with all methods discussed and compared in Sect. 8.1. Lastly, we put our findings in the context of similar cases and give an outlook of further fields of application.
8.1. Methodological results
Using the modelindependent approach, as detailed in Wagner & Tessore (2018) and Wagner et al. (2018), we found lensmodelindependent leadingorder ratios of convergences and reduced shear values for all multiple images in B0128 based on the positions of the three subcomponents in the images (see Sect. 4.2.4, Table 5, and Fig. 4). Subsequently, we succeeded in setting up a global freeform PixeLens lens model using the positions of the subcomponents as constraints (see Sect. 6.2). Due to scatterbroadening and a strong alignment of the subcomponents, the local lens properties of all approaches are subject to broad confidence bounds. The fact that both, the local and the global lens reconstructions suffer from the high degree of alignment of the subcomponents shows that an improvement in the astrometric precision for future observations will only yield a limited improvement in the lens reconstruction precision for highly aligned multiple images or highly aligned subcomponents. Our comparisons highlight the constraining power of the spatial configuration of the multiple images on the reconstruction precision, as also found in Wagner et al. (2018).
Within the 1σ confidence bounds, the modelindependent ratios of convergences and reduced shear values agree with the values obtained by the PixeLens reconstruction in all but one. So there is a similarly high degree of agreement between the modelindependent local lens properties and the modelbased values as was found in Wagner et al. (2018) for the galaxyclusterscale lens CL0024. Similar to Wagner et al. (2018), we conclude that the overall width of confidence intervals of the local lens properties is larger than those of the PixeLens modelbased reconstruction. This result shows the influence of the additional global regularisation constraints employed by PixeLens compared to the modelindependent method that forgoes a regularisation by being a local approach. The tendency of tightening confidence intervals for an increasing amount of additional model assumptions and of regularisation constraints is supported by the findings of Williams & Liesenborgs (2019) on galaxycluster scale as well.
Determining the ratios of convergences and reduced shear values at the individual positions of the subcomponents, that is, on milliarcsecond scale, we found larger mean reduced shear values than on imagescale. In addition, the local lens properties between the subcomponents within one image only overlap in 50% of the cases within their 1σ confidence bounds. The same degree of agreement was found when comparing the modelindependent local lens properties at subcomponent 1 to the ones at the same positions obtained by the parametric Lensmodel reconstruction using only subcomponent 1 in images A, C, and D as multiple image constraints. The high amount of required milliarcsecondsized pixels in PixeLens prevented us from setting up a PixeLens model to determine the local lens properties at the individual subcomponents with their confidence bounds. Hence, PixeLens is more robust than the modelindependent approach in returning tighter confidence bounds due to additional regularisation constraints. Vice versa, the modelindependent approach has the advantage over PixeLens that it is highly efficient in returning local lens properties and their confidence bounds at any scale with a minimum amount of computational effort.
On the whole, we conclude that the suitability of different global lens reconstruction approaches decisively depends on the resolution and the quality of the observations: on the scale of unresolved multiple images (for instance, for the data summarised in Table 1), the mass density distribution in B0128 still has elliptical symmetry, so that parametric lens models like Lensmodel are able to reproduce the multipleimage configuration, if the lens centre is not fixed. The resolved subcomponent structures in the multiple images reveal asymmetries in the deflecting mass density distribution which require more sophisticated freeform modelling approaches like PixeLens to explain the multipleimage configuration. Constraining local lens properties at the milliarcsecond scale of the subcomponents to probe smallscale dark matter properties is computationally more efficient to pursue with the modelindependent approach. It only yields local lens properties, meaning that it does not pursue a global reconstruction, but the local lens properties give the maximum information at leading order, which all lens models agree upon. Statistical screening methods that probe the symmetries of the observables of the multipleimage configuration like the one put forward in Woldesenbet & Williams (2012) can serve as consistency checks or provide initialisations for the global lens reconstruction approaches to increase the modelling efficiency.
We also note that, without including complementary nonlensing observables, the effective lens description employed here does not allow us to decide whether the asymmetries arising in the PixeLens reconstruction are caused by a smooth but asymmetric mass density distribution or by inhomogeneities in the mass density, for instance due to additional mass clumps located anywhere along the line of sight. Further details how to break this masssheet degeneracy can be found in Wagner (2019a).
8.2. Astrophysical conclusions for B0128
So far, only a few multipleimage configurations have been observed at milliarcsecond level and have been found to show subcomponents on this scale. The quadconfiguration in B0128 is one of these rare cases. B0128 is also special in a second way because its deflecting mass density distribution is a highredshift lenticular or latetype galaxy and not an earlytype one. As many recent works have consistently shown, see, for instance, Hsueh et al. (2018), Gomer & Williams (2018), Nightingale et al. (2019), and Gilman et al. (2019), smooth symmetric parametric lens models may not be a sufficient means to describe observed highly resolved multipleimage configurations on galaxyscale much longer. The findings summarised in Sect. 8.1 consistently show that B0128 is such an example.
Based on the PixeLens model of Sect. 6.2, we confirm the hypothesis stated in Xu et al. (2015) that the lens models as set up in Norbury (2002) or Biggs et al. (2004) are too simplistic to resolve the asymmetric deflecting mass density distribution that causes the multipleimages including their subcomponents on milliarcsecond scale. Even without assuming further deflecting mass clumps along the line of sight, it is not surprising that mass density isocontours change their morphology for increasing distance from the galactic centre, given the type of the deflecting galaxy. It also makes the higher reduced shear values at the subcomponentlevel and the potentially arising variations between the subcomponents look plausible. Yet the sparse constraints from the multiple images alone do not suffice to make more detailed statements about the smoothness of the mass density distribution. This is prevented by the masssheet degeneracy as detailed in Wagner (2018a,b, 2019a). Therefore, as stated in Sect. 8.1, complementary observational data is required to corroborate the assumption of Xu et al. (2015) that it is unlikely that smallscale dark matter inhomogeneities on top of a smooth mass density distribution cause the asymmetries observed in the mass density of B0128. Observations with increased sensitivity and resolution may also shed light on the contribution of the baryonic part to the asymmetry of the deflecting mass distribution.
Comparing the magnification ratios obtained by the modelindependent approach (see 𝒥_{i}, i = A, B, C, D, in Table 5), the magnification ratios as determined by PixeLens (see 𝒥_{i} in Table 7, and Fig. 4), and the observed flux ratios (see Table 3), we found a high degree of agreement between the observed flux ratios and the PixeLens values within the broad confidence bounds of the PixeLens reconstruction. The values based on the modelindependent approach have tighter confidence bounds and only agree for the subcomponent 3 in image D with the observed ones. While the comparison between the flux ratios and the PixeLens magnification ratios is considered over areas that have the same order of magnitude, the modelindependent approach determines the magnification ratios over the triangle spanned by the three subcomponent positions, which is less than 1% of the area of a PixeLens pixel. For the subcomponents, the differences between observed flux ratios and modelindependent magnification ratio values shrink, which corroborates the hypothesis that the different sizes over which the quantities are determined causes discrepancies between the results. Yet further investigations about the way that the flux ratios are calculated are necessary to confirm this. In addition, the observed flux ratios can be influenced by microlensing, scatterbroadening and absorption, see Biggs et al. (2004) and Lagattuta et al. (2010) for further details.
On the whole, we conclude that, at the current observational accuracy and precision, we have arrived at a consistent reconstruction of the deflecting mass density distribution and modelindependent local lens properties of B0128 which are able to explain the observed multipleimage configuration including the subcomponent structure on milliarcsecond scale. The findings are also in accordance with observations and modelling results of previous works.
The highprecision astrometry and existence of subcomponents in the radio bands give an unprecedented view of the galaxyscale lens. Evidence for deviations from a simple symmetric lens is found at arcsecond and milliarcsecond scales. On milliarcsecond scale, there are high shear values and shear gradients, which imply gradients in mass density. On arcsecond scale, simple parametric mass distributions, like Lensmodel models cannot reproduce images within the astrometric precision if the lens centre is fixed (see Sect. 5.2). Furthermore, the external shear required (as determined by Lensmodel and the modelfree approach in Fig. 6) is 0.22−0.26. Such large shears are unlikely to arise from nearby galaxies. For comparison, Bolton et al. (2008) modelled 63 SLACS lenses and found that shears range from 0 to 0.27, but the median is only 0.05. For B0128, the large shear value in a parametric symmetric lens model can equally well be modelled as a freeform asymmetric mass density distribution. Thus, it is possible that large shear values in other cases also subsume in them other complexities of the mass distribution, such as those suggested by Gomer & Williams (2018), and illustrated in their Fig. 14. Taking into account that Wagner (2018b) showed that splitting the total mass density into a main lens and a perturber introduces an additional masssheetlike degeneracy which allows for the redistribution of the overall mass between the two parts, this explanation seems highly probable. We will investigate further cases for a corroboration.
A similar case like B0128 is B1933+503, Cohn et al. (2001), Suyu et al. (2012), which is a spiral galaxy with ten multiple images of a threecomponent source. Observations at the resolution of milliarcsecond scale are not yet available. A second, similar case is the galaxygroup of B1349+154, Rusin et al. (2001), in which an unprecedented triangle of galaxies forms a gravitational lens that generates six multiple images of a background quasar. VLBA 1.7 GHz observations hint at potentially resolvable substructures in the brightest multiple images, so that this unique lensing configuration could also be an informative target. Further candidates for future smallscale radio observations and analyses could be selected from the forty quadrupoleimage configurations in Woldesenbet & Williams (2012) that show high distances to the Fundamental Surface of Quads.
As we showed in this work and in its predecessor on galaxycluster scale, Wagner et al. (2018), a joint analysis of multiple image configurations with complementary lens reconstruction methods offers the unique opportunity to obtain the most encompassing, consistent mass density reconstruction possible: extracting local lens properties at subcomponent scale is the computationally most efficient way to detect smallscale mass density gradients. Subsequently fitting a symmetric, parametric lens model or a freeform model to the sparse set of multiple images extends the local lens properties to the entire lensing region and allows to estimate the degree of symmetry of the lens. For a consistency check, the modelfree method of Woldesenbet & Williams (2012, 2015) can also be applied to constrain the degree of symmetry of the deflecting mass density distribution. Symmetric mass density distributions may be caused by relaxed mass agglomerations or by a low resolution observation that is unable to detect the deviations. Asymmetries in the mass density distribution can hint at merging mass distributions, further masses along the line of sight, or at a nonnegligible contribution of the baryonic part of the lens. While the asymmetries may be inferred from multiple image configurations alone, further, complementary observations, for instance, velocity dispersions along the line of sight, are necessary to determine the detailed origin of these asymmetries.
Thus, combining the observational data from current and future sky surveys, for instance, from the optical over infrared to the radiobands, and jointly employing complementary lens reconstruction approaches, as done for the case of B0128, we expect to gain an intriguing insight into a large ensemble of galaxyscale and galaxycluster scale gravitational lenses to infer smallscale dark matter properties as well as cosmological parameters.
Available at https://github.com/ntessore/imagemap
Lensmodel actually uses ϕ_{b} = br^{α}[1 − ϵcos(2(θ − θ_{ϵ}))]^{αβ} to do the calculations, which is not exactly the same form as presented in the Keeton (2004) manual.
Acknowledgments
The authors would like to thank Jori Liesenborgs for kindly agreeing to be the mediator of the unblinding process. JW gratefully acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG) WA3547/13.
References
 Biggs, A. D., Browne, I. W. A., Jackson, N. J., et al. 2004, MNRAS, 350, 949 [NASA ADS] [CrossRef] [Google Scholar]
 Bolton, A. S., Burles, S., Koopmans, L. V. E., et al. 2008, ApJ, 682, 964 [NASA ADS] [CrossRef] [Google Scholar]
 Cohn, J. D., Kochanek, C. S., McLeod, B. A., & Keeton, C. R. 2001, ApJ, 554, 1216 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, D., Birrer, S., Treu, T., Nierenberg, A., & Benson, A. 2019, MNRAS, 487, 5721 [NASA ADS] [CrossRef] [Google Scholar]
 Gomer, M. R., & Williams, L. L. R. 2018, MNRAS, 475, 1987 [NASA ADS] [CrossRef] [Google Scholar]
 Hsueh, J.W., Despali, G., Vegetti, S., et al. 2018, MNRAS, 475, 2438 [NASA ADS] [CrossRef] [Google Scholar]
 Keeton, C. R. 2001, ArXiv eprints [arXiv:astroph/0102341] [Google Scholar]
 Keeton, C. 2004, gravlens 1.06 Software for Gravitational Lensing, 9th edn., http://www.physics.rutgers.edu/~keeton/gravlens/manual.pdf [Google Scholar]
 Lagattuta, D. J., Auger, M. W., & Fassnacht, C. D. 2010, ApJ, 716, L185 [NASA ADS] [CrossRef] [Google Scholar]
 McKean, J. P., Koopmans, L. V. E., Browne, I. W. A., et al. 2004, MNRAS, 350, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Nightingale, J. W., Massey, R. J., Harvey, D. R., et al. 2019, MNRAS, 489, 2049 [NASA ADS] [Google Scholar]
 Norbury, M. A. 2002, PhD Thesis, The University of Manchester, UK [Google Scholar]
 Petters, A. O., Levine, H., & Wambsganss, J. 2001, in Singularity Theory and Gravitational Lensing (Birkhäuser), Prog. Math. Phys., 21 [CrossRef] [Google Scholar]
 Phillips, P. M., Norbury, M. A., Koopmans, L. V. E., et al. 2000, MNRAS, 319, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Rusin, D., Kochanek, C. S., Norbury, M., et al. 2001, ApJ, 557, 594 [NASA ADS] [CrossRef] [Google Scholar]
 Saha, P., & Williams, L. L. R. 2004, AJ, 127, 2604 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P., Ehlers, J., & Falco, E. E. 1992, in Gravitational Lenses (New York: Springer), Astron. Astrophys. Lib. [Google Scholar]
 Sluse, D., Chantry, V., Magain, P., Courbin, F., & Meylan, G. 2012, A&A, 538, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suyu, S. H., Hensel, S. W., McKean, J. P., et al. 2012, ApJ, 750, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Tessore, N. 2017, A&A, 597, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wagner, J. 2018a, A&A, 620, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wagner, J. 2018b, A&A, 615, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wagner, J. 2019a, MNRAS, 487, 4492 [NASA ADS] [CrossRef] [Google Scholar]
 Wagner, J. 2019b, Universe, 5, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Wagner, J., & Tessore, N. 2018, A&A, 613, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wagner, J., Liesenborgs, J., & Tessore, N. 2018, A&A, 612, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Williams, L. L. R., & Liesenborgs, J. 2019, MNRAS, 482, 5666 [NASA ADS] [CrossRef] [Google Scholar]
 Woldesenbet, A. G., & Williams, L. L. R. 2012, MNRAS, 420, 2944 [NASA ADS] [CrossRef] [Google Scholar]
 Woldesenbet, A. G., & Williams, L. L. R. 2015, MNRAS, 454, 862 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, D., Sluse, D., Gao, L., et al. 2015, MNRAS, 447, 3189 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Derivation of the reference point positions on the subcomponentscale
A.1. Derivation of coordinate positions
Figure A.1 depicts the quantities measured in the elliptical Gaussians fitted to the subcomponents: the length of the semimajor axis a, the axis ratio between the semiminor and the semimajor axis r, and the position angle θ measured with respect to the north axis. Denoting the centre of light of the Gaussian at the relative coordinates (Δα, Δδ) from a global reference point (see Table 2) as (α_{0}, δ_{0}), we obtained the four end points at the axes of the Gaussian by the following trigonometric relations:
Fig. A.1. Determining the reference points from an elliptical Gaussian fitted to the subcomponents: the centre of light is given by the relative coordinates Δα and Δδ in Table 2, here denoted by (α_{0}, δ_{0}). For images A and C with the same parity, we use the positions of the semimajor and semiminor axes as denoted by (α_{1}, δ_{1}) and (α_{2}, δ_{2}). Since image D has opposite parity, we employ (α_{1}, δ_{1}) and (α_{3}, δ_{3}) as reference point positions in addition to (α_{0}, δ_{0}). 
To estimate the uncertainty of the (α_{i}, δ_{i}), i = 1, …, 4, we assumed that the uncertainties in α_{0}, δ_{0}, r, and a were uncorrelated. Due to the fitting procedure, this is not the case. Yet it yields the order of magnitude to which the reference points can be determined. Given the uncertainties in Table 2, we found uncertainties on the order of 0.1 to 0.2 mas for the subcomponents in Table 2.
A.2. Impact of the uncertainties in the positions on the local lens properties
When mapping one elliptically Gaussian subcomponent to the one of another multiple image, the mapping should be independent of the reference points used. Only the relative parity between the multiple images must be obeyed. Thus, for the case of B0128, the configurations of reference points as listed in Table A.1 should all yield the same local lens properties. Determining the reference points from the semimajor and semiminor axes according to Appendix A.1 and calculating the local lens properties for all configurations listed in Table A.1, we corrobated this assumption.
Configurations of reference points in the subcomponents i = 1, 3 across images A, C, and D that should all yield the same local lens properties.
Without any loss of generality, we used configuration 0 of Table A.1 for all analyses described in Sect. 4.2.3.
Appendix B: Local lens properties for systematically interchanged subcomponent labels
Local lens properties as obtained using the three subcomponents in images A, C, and D with the matching across the images as listed in Table 4.
Appendix C: Local lens properties for the fourimage configuration ABCD
Local lens properties, their most likely value, mean, and the 1σ uncertainty bound (in Cols. 1, 2, and 3 of each configuration, respectively), as obtained using the three subcomponents in images A, B, C, and D with the matching across the images as listed in Table 4.
All Tables
MERLIN 5 GHz measurements as in Table 2 of Phillips et al. (2000) (Cols. 2–4), VLA 8.4 GHz measurements as in Table 1 of Phillips et al. (2000) (Cols. 5–7), and as in Table 2 of Lagattuta et al. (2010) using the Kpfilter of the Near Infrared Camera 2 on the Keck II telescope along with the LGS AO system (Cols. 8–10), both showing unresolved multiple images.
Observed flux density ratios ℱ_{i} with respect to image A for images i = B, C, D in different bands.
Configurations of reference points in the subcomponents i = 1, 3 across images A, C, and D that should all yield the same local lens properties.
Local lens properties as obtained using the three subcomponents in images A, C, and D with the matching across the images as listed in Table 4.
Local lens properties, their most likely value, mean, and the 1σ uncertainty bound (in Cols. 1, 2, and 3 of each configuration, respectively), as obtained using the three subcomponents in images A, B, C, and D with the matching across the images as listed in Table 4.
All Figures
Fig. 1. Left: HST Iband observation (HST Proposal 9133, based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (STECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA)) with B0128 (red circle) and the shear galaxy (yellow circle); the image detail shows the cleaned Iband image of B0128 as obtained by the CASTLe Survey, https://www.cfa.harvard.edu/castles/. Right: MERLIN 5 GHz observation by Phillips et al. (2000) showing the four multiple images; map peak: 0.0174 Jy beam^{−1}, contours: 0.000277 Jy beam^{−1} × (−3 3 6 12 24 48) amounting to multiples of 48 times the rootmeansquare (rms) noise of 277 μJy beam^{−1}, beam fullwidth at half maximum (FWHM): 58.8 × 47.7 mas^{2} at 50.9°. 

In the text 
Fig. 2. VLBA 8.4 GHz details of all four multiple images from Biggs et al. (2004), map peaks for images A, C, and D: around 1.4 mJy beam^{−1}, for image B: 0.39 mJy beam^{−1}, contours: 0.000225 Jy beam^{−1} × (−1 1 2 4 8 16 32 64 128 256 512), amounting to multiples of 3σ where σ is the offsource rms noise in the map (90 μJy beam^{−1}), beam FWHM: 1.2 × 0.7 mas^{2} at −21.6°. 

In the text 
Fig. 3. Visualisation of imagescale matching: using the positions of the three subcomponents (red, blue, green dots) in each image, the multiple images A, B, C, and D are mapped onto each other (red arrows, indicating that the dots of each colour are mapped onto each other). This mapping determines the local lens properties on the image scale (indicated by the larger red areas for visualisation purposes, as the convex hull spanned by the reference points in each image is too small to be drawn). Visualisation of subcomponentscale matching (highlighted in blue): using the Gaussians fitted to a subcomponent (here: for the first subcomponent denoted by 1, highlighted in dark grey), the subcomponents can be matched onto each other (as indicated by mapping the dots of each colour onto each other at the subcomponent scale) to determine the local lens properties within the area of the Gaussians. The first reference point in each Gaussian is its centre of light (yellow dots). The two reference points that are further required are the end points of the axes of the elliptical isocontour (brown, white dots). 

In the text 
Fig. 4. Graphical representation of the results presented in Tables 5 and 7. Each panel corresponds to one image of B0128, as labelled. Within each panel, we plot the values of g_{1}, g_{2}, f, and 𝒥 in 4 columns, separted by thin black vertical dotted lines. In each of these columns, there are 10 points, with confidence bounds, representing the four models from Table 5, and six models from Table 7 For image A, all values of f and 𝒥 are, by definition, 1. 

In the text 
Fig. 5. PixeLens reconstructed maps of lensing convergence. Left: using all twelve subcomponents of B0128. Right: using only the nine subcomponents of images A, C, and D. The subcomponents used as constraints in each case are indicated with magenta circles. The grey isoconvergence contours are spaced logarithmically. The blue contour is the isoconvergence contour at κ(x) = 1. 

In the text 
Fig. 6. Distribution of quads in the space of relative image angles. Here, a 2D projection of that 3D space is used. The Fundamental Surface of Quads is the horizontal line at Δθ_{23} = 0. Galaxyscale quads from Woldesenbet & Williams (2012) are shown as black circles with errorbars. The three quads of B0128 are represented by green squares, labeled with the corresponding subcomponent number. The orange points are quads from a synthetic lens, plotted here for reference, with projected density profile ∝r^{−0.8}, ellipticity ϵ = 0.25, and external shear γ = 0.25, which is misaligned with the ellipticity position angle by 80°. 

In the text 
Fig. A.1. Determining the reference points from an elliptical Gaussian fitted to the subcomponents: the centre of light is given by the relative coordinates Δα and Δδ in Table 2, here denoted by (α_{0}, δ_{0}). For images A and C with the same parity, we use the positions of the semimajor and semiminor axes as denoted by (α_{1}, δ_{1}) and (α_{2}, δ_{2}). Since image D has opposite parity, we employ (α_{1}, δ_{1}) and (α_{3}, δ_{3}) as reference point positions in addition to (α_{0}, δ_{0}). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.