© ESO 2020

1. Introduction

B0128+437 is a galaxy-scale 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 model-independent 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 free-form lens model approach PixeLens (Saha & Williams 2004) predict. Compared to previous lens models, we took into account the current best-fit redshifts of the lens and the source.

This comparison is analogous to the one carried out in Wagner et al. (2018) on the galaxy-cluster scale lens CL0024 and shows that the model-independent 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 model-independent 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 model-independent 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 non-smooth irregularly shaped mass density distributions on a pixelised grid. Therefore, such methods usually require computationally intensive pixelisations to determine these small-scale 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 small-scale details of the lens. Thus, a more efficient way to obtain small-scale 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 model-independent, local lens properties and that were used to constrain our lens models in Sect. 3. In Sect. 4, we briefly outline the model-independent algorithm, which is further detailed in Wagner & Tessore (2018) before we show the model-independent, local lens characterisation of B0128. Analogously, Sects. 5 and 6 describe the lens modelling based on the parametric code Lensmodel, and the free-form 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 model-independent and the model-based approaches and only revealed for the comparison at the very end, after all modelling is complete. In Sect. 7, we apply the model-free 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 model-independent and the model-based lens reconstructions.

2. Related work on B0128

In the discovery paper, Phillips et al. (2000), high-resolution MERLIN observations at 5 GHz from B0128 were obtained in the course of the Cosmic Lens All-Sky 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 quad-lens formation were detected (see Fig. 1 (right)). Long-term observations showed no hint of a time-variability of the quasar.

Fig. 1. Left: HST I-band 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 (ST-ECF/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 I-band 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 root-mean-square (rms) noise of 277 μJy beam−1, beam full-width at half maximum (FWHM): 58.8 × 47.7 mas2 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 model-predicted 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 follow-up 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 SIE-lens 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, follow-up 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 scatter-broadened and hence, it is difficult to decompose it into subcomponents. They also performed VLBI-observations 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 sub-milliarcsecond 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 off-source rms noise in the map (90 μJy beam−1), beam FWHM: 1.2 × 0.7 mas2 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 late-type 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 sub-mas precision for the image positions is required. To investigate the most probable source of the discrepancies that enter at sub-mas 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 scatter-broadening is more likely to cause the flux anomaly than additional small-scale dark-matter halos.

Summarising the previous results, observations indicate that

– B0128 is a potentially lenticular or late-type 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 sub-components, the fourth, image B in Fig. 1 (right), is most likely scatter-broadened,

– 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 sub-components to sub-mas-precision 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).

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 three-image 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 mirror-inverted 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₁ in the north and D₃ in the south. Matching the subcomponents by their intensity led to this labelling. Yet, the multi-band 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 scatter-broadening. 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 model-independent 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 model-independent approach and the lens model as a consistency check because we did not employ them as constraints.

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 scatter-broadening due to the dust in the lens, as already noted by Biggs et al. (2004).

4. Model-independent reconstruction

4.1. The method

Effectively describing the gravitational lens as a two-dimensional 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 two-dimensional mass density at the position x ∈ ℝ² 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) = (γ₁(x),γ₂(x)) additionally distorts the images. Both κ(x) and γ(x) can be determined as second-order derivatives from the two-dimensional 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 model-based 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,

(1)

between all multiple images i, j and reduced shears,

(2)

at the positions of the multiple images. These local lens properties in Eqs. (1) and (2) are invariant under the mass-sheet 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_ijs and g_is, the magnification ratios,

(3)

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

(4)

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 C-implementation¹ 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 so-called “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 χ²-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 non-aligned 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 under-constrained 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 three-image 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 three-image 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 four-image 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 three-image 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

(5)

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 semi-major and semi-minor axis of the fitted Gaussian as further reference points in each subcomponent. The calculations how to obtain these positions from the values of the semi-major 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 image-scale local lens properties. The matching for the image- and subcomponent-scale is sketched in Fig. 3.

Fig. 3. Visualisation of image-scale 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 subcomponent-scale 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 image-scale 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 cross-checked. 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 image-scale 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. Image-scale 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 g1, g2, 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.

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. Image-scale 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 scatter-broadening of image B, which is not accounted for in our local lens reconstruction.

Swapping the labelling of subcomponent B₁ with B₃, 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 four-image 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₁ and B₃ and the increase of the measurement uncertainties in Appendix C.

4.2.3. Subcomponent-scale matching of subcomponents 1 and 3

After the image-scale local lens properties were obtained, we turned to the subcomponent-scale 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 semi-major and semi-minor 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 image-scale local lens properties, all f- and g-values 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 four-image 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 galaxy-cluster-scale 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 scatter-broadening 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 multi-band observations, it is hard to disentangle the impact of the individual effects.

Comparing the image-scale 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 four-image configurations, if no sub-mas precision of the positions of the subcomponents is assumed.

5. Model-based 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

(6)

where b is the normalisation constant, s is the core radius, which is set to a small non-zero 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

(7)

where ϵ is the ellipticity and θ_ϵ is its position angle². 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 image-scale, 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₁C₁D₁ and config A₁B₁C₁D₁ 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 model-predicted ones in the lens plane. If the lens centre was held fixed at its observed position, config A₁C₁D₁ could still reproduce the images perfectly (third column of the first row). If config A₁B₁C₁D₁ 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.

The resulting Lensmodel local lens properties for config A₁B₁C₁D₁ and config A₁C₁D₁, 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₁B₁C₁D₁ 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 milli-arcsecond 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 HST-sized pixels in the lens plane, there would be no indication that simple models do not fit.

6. Model-based reconstruction: free-form

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 free-form method, called PixeLens (Saha & Williams 2004), to reconstruct the mass density distribution in B0128. PixeLens is publicly available, and has an easy-to-use GUI interface³.

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 “burn-in” 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 one-off 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 milli-arcseconds. 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 iso-convergence contours are spaced logarithmically. The blue contour is the iso-convergence 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 most-likely scatter-broadened 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 milli-arcseconds, 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

(8)

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 g-values 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. Model-free analysis

7.1. The method

Finally, we performed another type of analysis on the three sets of sub-components 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 quadrupoly-imaged 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 θ₁₂, θ₂₃, and θ₃₄. 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 θ₂₃ versus the distance of quads from the FSQ, Δθ₂₃, which is measured parallel to the θ₂₃ axis. Note that θ₂₃ 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 sub-components in B0128 as green squares, together with forty galaxy-scale 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. Galaxy-scale 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 non-elliptical 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 galaxy-lens 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 non-elliptical 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 galaxy-scale observed quads. They concluded that observed deviations from FSQ by Δθ₂₃ ∼ 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 non-negligible perturbations from ellipticity in the case of B0128, indicated by Fig. 6 is consistent with the mass distribution produced by free-form PixeLens, in Sect. 6, and shown in Fig. 5.

8. Conclusion

We presented four different types of analyses and reconstructions of the galaxy-scale lenticular or late-type gravitational lens B0128 constrained by the quadrupole-image configuration of a background quasar. Previous multi-band 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 milli-arcseconds, 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 multiple-image configuration at subcomponent-scale 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 model-independent approach, as detailed in Wagner & Tessore (2018) and Wagner et al. (2018), we found lens-model-independent leading-order 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 free-form PixeLens lens model using the positions of the subcomponents as constraints (see Sect. 6.2). Due to scatter-broadening 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 model-independent 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 model-independent local lens properties and the model-based values as was found in Wagner et al. (2018) for the galaxy-cluster-scale 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 model-based reconstruction. This result shows the influence of the additional global regularisation constraints employed by PixeLens compared to the model-independent 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 galaxy-cluster scale as well.

Determining the ratios of convergences and reduced shear values at the individual positions of the subcomponents, that is, on milli-arcsecond scale, we found larger mean reduced shear values than on image-scale. 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 model-independent 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 milli-arcsecond-sized 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 model-independent approach in returning tighter confidence bounds due to additional regularisation constraints. Vice versa, the model-independent 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 multiple-image 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 free-form modelling approaches like PixeLens to explain the multiple-image configuration. Constraining local lens properties at the milli-arcsecond scale of the subcomponents to probe small-scale dark matter properties is computationally more efficient to pursue with the model-independent 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 multiple-image 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 non-lensing 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 mass-sheet degeneracy can be found in Wagner (2019a).

8.2. Astrophysical conclusions for B0128

So far, only a few multiple-image configurations have been observed at milli-arcsecond level and have been found to show subcomponents on this scale. The quad-configuration in B0128 is one of these rare cases. B0128 is also special in a second way because its deflecting mass density distribution is a high-redshift lenticular or late-type galaxy and not an early-type 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 multiple-image configurations on galaxy-scale 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 multiple-images including their subcomponents on milli-arcsecond 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 subcomponent-level 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 mass-sheet 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 small-scale 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 model-independent 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 model-independent 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 model-independent 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 model-independent 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, scatter-broadening 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 model-independent local lens properties of B0128 which are able to explain the observed multiple-image configuration including the subcomponent structure on milli-arcsecond scale. The findings are also in accordance with observations and modelling results of previous works.

The high-precision astrometry and existence of subcomponents in the radio bands give an unprecedented view of the galaxy-scale lens. Evidence for deviations from a simple symmetric lens is found at arcsecond and milli-arcsecond scales. On milli-arcsecond 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 model-free 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 free-form 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 mass-sheet-like 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 three-component source. Observations at the resolution of milli-arcsecond scale are not yet available. A second, similar case is the galaxy-group 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 small-scale radio observations and analyses could be selected from the forty quadrupole-image 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 galaxy-cluster 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 small-scale mass density gradients. Subsequently fitting a symmetric, parametric lens model or a free-form 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 model-free 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 non-negligible 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 radio-bands, 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 galaxy-scale and galaxy-cluster scale gravitational lenses to infer small-scale dark matter properties as well as cosmological parameters.

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/1-3.

References

Appendix A: Derivation of the reference point positions on the subcomponent-scale

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 semi-major axis a, the axis ratio between the semi-minor and the semi-major 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 (α₀, δ₀), we obtained the four end points at the axes of the Gaussian by the following trigonometric relations:

(A.1)

(A.2)

(A.3)

(A.4)

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 semi-major and semi-minor 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 α₀, δ₀, 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 semi-major and semi-minor axes according to Appendix A.1 and calculating the local lens properties for all configurations listed in Table A.1, we corrobated this assumption.

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

Appendix C: Local lens properties for the four-image configuration ABCD

All Tables

All Figures

Fig. 1. Left: HST I-band 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 (ST-ECF/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 I-band 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 root-mean-square (rms) noise of 277 μJy beam−1, beam full-width at half maximum (FWHM): 58.8 × 47.7 mas2 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 off-source rms noise in the map (90 μJy beam−1), beam FWHM: 1.2 × 0.7 mas2 at −21.6°.
In the text

Fig. 3. Visualisation of image-scale 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 subcomponent-scale 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 g1, g2, 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 iso-convergence contours are spaced logarithmically. The blue contour is the iso-convergence 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. Galaxy-scale 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 semi-major and semi-minor 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