Modelling strong lenses from wide-field ground-based observations in KiDS and GAMA Free

ABSTRACT

1 INTRODUCTION

Strong gravitational lensing is an essential probe of galaxy structure, enabling mass measurements in the centre-most regions of the foreground lensing galaxy without assumptions about stellar populations. Numerous studies have shown that lensing galaxies are, in every other respect, indistinguishable from other galaxies in the observed mass range; therefore, their study offers insight into the larger global population of galaxies at similar mass and redshift (Auger et al. 2009). Complementary to kinematic and stellar population synthesis measurements, strong lensing allows the decoupling of internal mass components (dark and baryonic) and accurate central stellar population mass-to-light ratios (Auger et al. 2010b; Hopkins 2018).

A fundamental issue in astronomy is relating the predicted dark matter halo mass function to the observed galaxy mass function. The masses of dark matter haloes are not well constrained by the amount of visible matter in their constituent galaxies. Few single galaxies corresponding to the lowest and highest halo masses are found, due to feedback processes stopping star formation (Behroozi , Conroy & Wechsler 2010; Behroozi et al 2020). The galaxy stellar-to-halo mass relation represents a fundamental barometer for accretion and feedback processes in galaxy formation. Subhalo abundance matching assigns a galaxy with a specific stellar mass to a specific subhalo but does not consider (i) the enveloping host halo mass or (ii) whether the galaxy is a central or a satellite; thus it suffers from assembly bias (Zentner, Hearin & van den Bosch 2014; Chaves-Montero et al. 2016).

Assembly bias is a secondary halo property that is related to the clustering strength of haloes (Matthee et al. 2017; Zehavi et al. 2018), where the clustering of dark matter haloes depends on their mass and formation epoch. Investigations with cosmological simulations have revealed that dark matter halo concentration, formation time, and environment all play a role in the relation between the central galaxy’s stellar mass and the mass of the dark matter halo it occupies. Assembly bias appears to be mostly independent of the cosmological parameters assumed (Contreras et al. 2021). At high masses typical of lensing elliptical galaxies, inefficiency in the stellar occupancy of dark matter haloes is ascribed to the effects of active galactic nuclei (Somerville, Popping & Trager 2015). Weak-lensing studies (Mandelbaum et al. 2013; Velander et al. 2014; Lin et al. 2016) and velocity field studies (McCarthy et al. 2021; Posti & Fall 2021) appear to show the effects of assembly bias on the scales of galaxy populations (Cui et al. 2021); assembly bias is especially noticeable at ∼1 Mpc scales (Hearin et al. 2016): i.e. groups of galaxies. While well-explored with hydrodynamical simulations (e.g. Hearin, Watson & van den Bosch 2015; Hearin et al. 2016; Matthee et al. 2017; Artale et al. 2018; Zehavi et al. 2018, 2019; Contreras et al. 2019), observational studies have thus far been limited by the need to average over large numbers of similar-mass central elliptical galaxies to obtain a weak lensing or velocity signal.

With strong lensing, one has the opportunity to directly measure stellar and halo masses in elliptical galaxies. Relations between the environment and internal structure of elliptical galaxies have been explored using Sloan Lens ACS (SLACS; Bolton et al. 2006) by Treu et al. (2010). They find that the SLACS lenses are slightly biased towards overdense environments (12 of 70 are associated with known groups or clusters), which is consistent with the expectation for the most massive of elliptical galaxies. They find this result to be unbiased when compared to similar massive galaxies from SDSS, again showing lens galaxies to be representative of the overall elliptical galaxy population. They find the contribution of the external environment to have little effect on the local potential (except in extreme overdensities) and the internal structure of lens galaxies. SLACS and other lens studies have been conducted using detailed observations of individual lenses with HST-quality data. The application of lens modelling methods to larger wide-field surveys offers an alternative avenue with advantages for conducting experiments relating galaxy properties to environment and group properties. This motivates the need for exploring lens modelling methods in the context of large surveys.

In this paper, we explore what can be done with ground-based imaging and spectroscopy to model lens candidates after they have been identified in imaging surveys. We discuss strategies for ensuring quality control at each stage while extracting meaningful measurements from ground-based data. Using Galaxy and Mass Assembly (GAMA) (Driver et al. 2009, 2011; Liske et al. 2015) survey single-aperture spectroscopy, we explore the utility of automated redshift determination as a tool for identifying the background-source redshifts of strong lenses by applying this technique to lens candidates that were identified in the ground-based imaging of the Kilo-Degree Survey (KiDS) (de Jong et al. 2013, 2015, 2017; Kuijken et al. 2019) using machine learning techniques (Petrillo et al. 2019b; Li et al. 2020). With GAMA spectroscopic redshifts and other measurements in conjuction with KiDS cutout images, we construct lens models utilizing an automated lens modelling program called PyAutoLens (Nightingale et al. 2021b).

Our paper is organized as follows: Section 2 describes the KiDS and GAMA data used, as well as the parent samples used in our selection. Section 3 describes how background-source redshifts are identified in single-aperture spectra from GAMA Autoz catalogs to create a subsample for modelling. Section 4 describes the pyautolens software and the lens modelling methods we used to perform our analysis. Section 5 outlines the assessment of quality of the models and redshift determinations. Section 6 presents results for the four highest-quality models. Section 7 discusses some challenges that our prescription of second-redshift determination introduced, as well as recommendations for improving that method. Section 8 discusses galaxy environment and potential applications of a refined method to future studies. Section 9 lists our conclusions. Throughout this paper, we adopt a Planck Collaboration 2015 (2016) cosmology (H₀ = 67.7 km s⁻¹ Mpc⁻¹, Ω_m = 0.307).

2 DATA

2.1 GAMA spectroscopy and AUTOZ redshifts

GAMA (Driver et al. 2009, 2011; Liske et al. 2015) is a multiwavelength survey built around a deep and highly complete redshift survey of five fields with the Anglo-Australian Telescope. GAMA has three major advantages over SDSS in the identification of blended spectra: (i) the spectroscopic limiting depth is 2 mag deeper (m_r < 19.8 mag compared with SDSS main survey depth m_r < 17.7 (Eisenstein et al. 2001)¹), (ii) the completeness is close to 98 per cent (Liske et al. 2015), and (iii) the autoz redshift algorithm can identify spectral template matches with signal from two different redshifts (Baldry et al. 2014). These properties and the overlap of GAMA and KiDS fields make these two surveys exceptionally well suited to provide the data required for our study of lens modelling.

The autoz (Baldry et al. 2014) cross-correlation redshift software has been uniformly applied to the GAMA (Liske et al. 2015) spectroscopic data, resulting in a public data base that can be found in GAMA-DR3 aatspecautozall v27 (hereafter autoz catalog) and specall v27 tables (Liske et al. 2015, http://www.gama-survey.org/dr3/). The autoz algorithm outputs four flux-weighted cross-correlation peaks (denoted σ) of redshift matches to template spectra of emission line and passive galaxies (denoted ELG and PG, respectively) from SDSS-DR5. σ₁ corresponds to the highest cross-correlation or ‘best-fitting’ redshift match, σ₂ the second-best, etc. These matches have proven to be highly successful and are the base redshift measurement for GAMA objects. GAMA-DR3 also compiled SDSS-BOSS spectra for overlapping targets that are included in table specall, but these spectra did not utilize autoz for redshift determination.

The completeness of GAMA allows detailed environment measures, including population density and separation (Brough 2011; Alpaslan et al. 2014, 2015), and the GAMA team internal groupfinding catalogs include the total mass and placement of the galaxies in an identified group via a friends-of-friends algorithm (Robotham et al. 2011). In fact, GAMA was conceived to probe the effects of group environment on galaxy properties. In this context, we describe galaxies either as ‘group member’ galaxies or as those not in galaxy groups, which we designate as ‘isolated galaxies’. Stellar masses are taken from the GAMA-DR3 stellarmasseslambdar v20 catalog (Taylor et al. 2011).

2.2 Kilo-Degree Survey and machine learning strong lens samples

The KiDS (de Jong et al. 2013, 2015, 2017; Kuijken et al. 2019) is a VLT Survey Telescope (VST) program of medium-deep imaging in SDSS-ugri filters primarily to identify weak lensing. The deep imaging, high resolution (0.65 arcsec seeing in SDSS r-band), and wide sky-coverage (1350 deg²) also make this survey ideal for efforts to identify strong lens candidates from imaging. Image-based deep learning efforts have been the most promising of recent developments in automated lens-finding algorithms. Their efficiency and versatility make them ideal for astronomical classification problems involving large data sets, including the detection of strong gravitational lenses [e.g. in Subaru Hyper-Supreme Cam (Speagle et al. 2019), DESI DECam Legacy Survey (Huang et al. 2020), and Dark Energy Survey data (Jacobs et al. 2019)]. Petrillo et al. (2017) developed a machine learning method to identify strong lenses in KiDS using a convolutional neural network (CNN) with artificially constructed lens images as the training set. Training and target catalogs were intentionally constructed utilizing SDSS-luminous red galaxy (LRG; Petrillo et al. 2017) colour-magnitude selection cuts to return the largest of known strong lenses that result in the most readily identifiable lens features (Einstein radii close to and greater than 1 arcsec). The result is the Lenses in KiDS sample (LinKS; Petrillo et al. 2019a, b)² of some 1300 strong lensing candidates, 421 of which overlap with the equatorial regions of the GAMA survey. Knabel et al. (2020) compared data from LinKS objects with the autoz spectroscopic identifications in GAMA (Holwerda et al. 2015) as well as with KiDS-GalaxyZoo (Holwerda et al. 2019, Kelvin et al. in preparation) citizen science identifications in the overlapping equatorial fields. A disparity between the candidate samples in terms of stellar mass and redshift is attributed to selection effects. Of the subsample of 421 LinKS candidates in GAMA (hereafter referred to as ‘LinKS’ or ‘LinKS in GAMA’ sample), there was no overlap with the subsample of 47 GAMA spectroscopic candidates. Knabel et al. (2020) further subselected 47 LinKS candidates to represent the highest quality of the sample (hereafter referred to as ‘LinKS from Knabel-2020’ sample).

Li et al. (2020) followed a slightly modified approach from the lens-search prescription utilized by the LinKS team to search for strong lens candidates in KiDS-DR4. They included several more LRGs and applied their CNN to a sample of ‘bright galaxies’ (BG) that did not undergo LRG colour-magnitude cuts. Their search returned some LinKS candidates and resulted in 286 new candidates within the KiDS survey, 48 of which were identified in the GAMA equatorial regions. 39 of those have matches in the stellarmasseslambdar mass catalog, and there are no overlaps between this sample and that of GAMA spectroscopy or GalaxyZoo. This candidate sample, which we will refer to as ‘Li-BG’, shares essentially the same parameter space as LinKS candidates, even with the exclusion of the LRG selection (Knabel et al. 2020). This is not necessarily surprising considering Petrillo et al. (2019b) report no significant advantage to the inclusion of colour images in the CNNs, as the networks appear to focus more on morphological features and brightness than colour separation. Still, the extension beyond the typical red elliptical galaxy as candidate objects suggests the potential for more variability in candidate characteristics.

3 AUTOZ SECOND REDSHIFT SELECTION AND QUALITY CONTROL

We examine the Autoz cross-correlation values for each LinKS and Li-BG lens candidate. Each candidate has been matched to the closest GAMA object by right-ascension and declination within a positional tolerance of 2 arcsec. Not every object in the equatorial fields is featured in the autoz catalog; these candidates are removed from this study. Some of the objects feature duplicated entries, some of which have conflicting autoz outputs, which we retain for examination and selection.

3.1 Selection criteria

Since the candidates that remain have already been identified and vetted through machine learning methods, we adopt a more lenient selection criterion from the same σ₂ − R parameter space as that utilized in Holwerda et al. (2015). From a first look at the data, we select candidates with R ≥ 1.2. This is sufficient for a first selection and for characterizing the output of the autoz algorithm from already positively identified candidates. We show the selection in Fig. 1.

Figure 1.

Initial selection of candidates with autoz second-redshift determinations. The y-axis shows the σ₂ second cross-correlation peak, with higher values indicating a stronger match to the second galaxy template. The x-axis is the parameter R, given by equation (1). Black markers indicate 348 LinKS Autoz entries (300 unique candidates), 51 of which are overplotted with green markers to indicate autoz entries of the 47 unique LinKS candidates as selected in Knabel et al. (2020). Orange markers indicate the 53 autoz entries for the 32 unique Li-BG candidates. The dashed vertical line shows R = 1.2, and the dotted box encloses the area of parameter space used by Holwerda et al. (2015).

Red squares surround 59 LinKS and 8 Li-BG entries that satisfy R ≥ 1.2 and are followed with additional selection criteria. See Section 3.1.

Open in new tabDownload slide

From the 67 Autoz entries that pass the R selection, we remove those with stellar template matches (i.e. not a galaxy spectrum) and retain all those with galaxy-galaxy template match configurations regardless of the galaxy type. The distribution of foreground+background (lens+source) template type (PG + ELG, ELG + PG, ELG + ELG, PG + PG) is shown in the histogram of Fig. 2. Note that the majority of lens foreground objects match to PG templates, with the majority of background objects matching to ELG templates. This is expected. Massive elliptical galaxies tend to be the most readily observable strong lensing foreground objects, and bright emission lines from the background source are the most easily detected behind a PG continuum. This selection bias is further enhanced by the fact that the parent LinKS and Li-BG samples were both identified by CNNs trained with large elliptical galaxies. Configurations with foreground lens ELGs are possible, though much more difficult to detect. The reasons are: (i) ELGs are typically lower in mass, so the lensing is less pronounced, (ii) lensed background sources are often bright ELGs, so the blue light of each can blur together, and (iii) ELGs can include complex morphologies that can be mistaken for lens features. The majority of Li-BG candidates in the autoz catalog were matched to ELGs in the foreground, but further selection and assessment of the spectra showed this to be a false trend. As in Knabel et al. (2020), we select candidates with background source redshifts that are reasonably far away from the foreground lens redshifts, as well as those whose autoz foreground redshifts are greater than 0.05. autoz estimates the probability of success of the primary redshift match, and one candidate is removed due to its very low probability.

Figure 2.

autoz galaxy template combinations of candidate subsamples, listed as foreground + background. Black and orange refer to LinKS and Li-BG candidates, respectively, and are stacked to show total counts of each template combination. The green outline shows the subset of LinKS candidates that were examined in Knabel et al. (2020). As expected, the majority of lens foreground objects match to PG templates, while the majority of background objects match to ELG templates.

Open in new tabDownload slide

Described in more detail and discussed in the context of Knabel et al. (2020) in Appendix  E, our selection results in 42 lens candidates (39 from LinKS and 3 from Li-BG) with second redshift matches, which constitute the initial Autoz sample. In later work, it may be worth exploring machine learning algorithms to find better use of the parameter space for classification than our naive selection criteria.

3.2 AUTOZ selection quality control

An initial modelling and analysis strategy utilizing the autoz sample of 42 candidates and pyautolens revealed the need for further quality control. This assessment addresses the validity of our application of autoz output parameters to select second-redshift determinations for the use in strong lens modelling. We examined the 42 candidate spectra to understand the evidence for a second redshift match within each. We superimposed upon the candidate spectra important emission and absorption features redshifted by the lens and source redshifts determined by autoz. For ELG template matches, we inspected Hβ, [O ii] λ3727, and [O iii] λλ4959,5007 emission lines. For PG template matches, we looked at Calcium H and K, Hβ, Mg, and Na absorption lines. This illuminated some specific cases that pointed to dubious redshift matches. We identify four such cases:

3.2.1 When there are overlapping emission or absorption lines...

autoz template matches utilize strong absorption or emission features to identify passive and star-forming galaxies. We discovered a significant failure condition of our method that occurs when the second redshift match is identified by an emission or absorption line that is also attributed to the foreground lens galaxy redshift. This is particularly obvious in cases where Autoz determined ELG + ELG configurations (i.e. both redshifts are identified by emission lines), in which many of the higher-redshift matches were determined by an [O ii] λ3727 line that overlapped with the [O iii] λλ4959,5007 lines of the foreground lens galaxy. For configurations where both the lens and source are the same template type (i.e. PG+PG or ELG + ELG), overlapping line features usually indicate a poor second-redshift determination.

However, emission and absorption features are present in the spectra of both PGs and ELGs. For example, Treu et al. (2002) found that |$\sim 10{{\ \rm per\ cent}}$| of elliptical galaxies in the intermediate redshift range where our lens candidates lie have strong [O ii] λ3727 emission; in fact, the emission of [O ii] λ3727 is detectable in the spectra 14 of the PG matches in our sample, and absorption of Calcium H and K lines (as well as some Na and Mg) are identifiable in 15 of the ELG galaxy spectrum matches. In most cases, the overlap of lines draws suspicion of the background source spectrum match, not the foreground lens match, which is typically the primary template match (σ₁). In our sample, emission line overlaps occur exclusively between the background source [O ii] λ3727 and one of the [O iii] λλ4959,5007 couplet of the foreground lens PG. Absorption line overlaps are all between source H or K lines and lens Hβ, Na, or Mg lines.

3.2.2 When the lens is described as an emission line galaxy...

As shown in Fig. 2, several of the template matches selected by the strategy described in Section 3.1 are configurations with a foreground lens ELG (ELG + ELG or ELG + PG). We reiterate that there is no physical reason to mistrust these redshift matches on this fact alone. In fact, given that several of the ELG matches have strong absorption features at the same redshift, these may very well be elliptical galaxies with strong oxygen emission lines. However, upon examination of candidate spectra and the quality of initial models for these ELG + lens configurations, we found that several of these candidates included dubious template matches and were removed.

The overlapping ELG+ELG emission lines have been discussed. In addition, some of the spectra of ELG+ configurations included key emission line features that would have been observed in wavelength ranges of high noise. Many of the GAMA spectra have their most significant noise at the extremes of the wavelength range, e.g. at wavelengths shorter than ∼4500 Å. For lenses at redshifts lower than ∼0.2, the [O ii] |$\lambda 3727{\, \mathring{\rm A}}$| emission line lies in this region, which removes an important identifying emission line feature from consideration.

3.2.3 When source emission lines are redshifted beyond observed wavelength range...

The cases where a background source emission line overlaps with a foreground lens emission line often correlate with another case: some of the background source emission lines have redshifted to wavelengths beyond the observed wavelength range of the spectroscopic survey. For our purposes, this translates to upper limits of background source redshifts beyond which the emission line will not be present in the spectrum. For GAMA, which has an upper limit of 8850 Å, the emission lines we have discussed begin to disappear around z ∼ 0.77. SDSS-BOSS has an extended upper limit of wavelength to 10 400 Å (∼J-band), which corresponds to upper redshift limits of around z ∼ 1 where we begin to see missing emission lines. GAMA-DR3 specall table contains SDSS-BOSS spectra for several of the candidates, so we can look to these spectra for evidence of lines that have redshifted beyond the GAMA upper limit. The autoz catalog includes only GAMA spectra, so the outputs we use for selection would not benefit from the extended range of the SDSS-BOSS spectra.

3.2.4 When primary redshift is background source...

The primary redshift match corresponding to σ₁ is typically (but not always) the foreground lens galaxy. The background source redshift is the primary redshift match for 10 of the 42 autoz sample candidates. The results of the Autoz algorithm are largely flux-weighted, and an especially bright background source (e.g. a strong emission line) could be interpreted by the algorithm as the primary redshift solution instead of the lower-redshift lens object. However, 8 of these 10 are ELG + PG template configurations, where a PG template gives the primary redshift solution at higher redshift. The case of a bright continuum at higher redshift overshadowing the foreground ELG is unlikely.

All candidates in the autoz sample were modelled and examined in the manner described in Sections 4–6 before these four cases were identified. In Section 7, we discuss the results of modelling and assessment in the context of the cases described in this section and recommend alterations to the initial selection scheme. We note that a poor redshift match does not remove/negate the validity of the lens identification, nor does it question the accuracy of the uniform application of autoz to GAMA spectroscopic targets. These additional selection decisions were instituted following critical assessment of problems with our initial strategy that required time with human eyes on the spectra. With reasonable background source redshift determinations, modelling of the imaging data can yield meaningful physical measurements.

4 PYAUTOLENS

We use the open-source lens modelling software pyautolens (Nightingale et al. 2021b)³ to perform our analysis. The software is described in Nightingale, Dye & Massey (2018) and Nightingale et al. (2021b), building on the works of Warren & Dye (2003), Suyu et al. (2006), and Nightingale & Dye (2015). We refer readers to these works for a full description of our approach to lens modelling. Section 4.1 broadly describes the method as we apply it here, and a more technical description of the specifics of the implementation is given in Appendices  A and  B.

4.1 Lens modelling with pyautolens

pyautolens models the foreground lens galaxy’s light and mass as well as the background source galaxy’s light simultaneously. First, pyautolens assumes a profile for the foreground lens’s light (e.g. a Sérsic profile), producing a model image of the lens galaxy. A mass model then ray-traces a grid of image-pixels from the image-plane to the source-plane, with the source’s light evaluated on this deflected grid via another light profile. This creates an image of the lensed source, which is added to the lens galaxy image to create an overall model image of the strong lens. This image is convolved with the instrument point spread function (PSF) and compared to the data to evaluate the residuals and likelihood of that lens model. To fit a lens model to imaging data, pyautolens searches an N-dimensional parameter space so as to minimize the residuals (and therefore maximize the likelihood) between the model image and the observed image. The lens models fitted in this work consist of N = 7–14 parameters, corresponding to the parameters of the light and mass profiles that represent the lens and source galaxies. To sample parameter space, we use the nested sampling algorithm dynesty (Speagle et al. 2019), and we detail its specific implementation below.

The parameter spaces of a strong lens model are challenging to sample, and local maxima and unphyscial lens models are often inferred. Automating the model-fitting procedure is therefore difficult, and pyautolens approaches automation via a technique called non-linear search chaining. Here, a sequence of dynesty model-fits are performed that fit lens models of gradually increasing complexity, whereby the results of the initial searches are used to inform the search of more complex parameter spaces in the later searches. Through experimentation, we have designed a pipeline composed of a chain of three dynesty searches that we use as a template for fitting each lens. Our three-step pipeline consists of three sequential model-fits:

• Search 1 - Lens Light: models only the foreground lens elliptical light profile.

• Search 2 - Lens Mass and Source Light: focuses on the background source light profile and lensing deflections.

• Search 3 - Combined Lens and Source Models: models each component in the system for parameter inference.

Between each search, various aspects of the fit can be altered (e.g. a mask applied to the data can be customized to show only the specific features of interest to each fit). This offers a more efficient lens modelling procedure overall, as the parameter spaces of reduced complexity are sampled faster.

Search chaining uses a technique called ‘prior passing’ to initialize the regions of parameter space that are searched later in the chain. Here, the models inferred in earlier non-linear searches initialize the priors of the more complex models fitted by the searches later on. This ensures the non-linear search samples only the higher likelihood regions of parameter space [see Nightingale et al. (2018)] and therefore reduces the probability that a local maximum is inferred. Prior passing sets the prior of each parameter as a Gaussian. The mean is that parameter’s previous inferred median PDF value, and the width is a value specific to each lens model and parameter. Prior widths have been carefully chosen to ensure they are broad enough not to omit lens model solutions by trimming valid solutions but sufficiently narrow to ensure the lens model does not inadvertently infer local maxima.

The dynesty nested sampling algorithm (Speagle et al. 2019) can also balance efficiency in computation with how thoroughly it explores parameter space. Initial model fits require only a rough estimate of the lens model that provides a reasonably approximate fit to the data. These searches therefore use faster dynesty settings that give a less thorough sampling of parameter space. Our final results require accurate and robust parameter estimates with precise and well-quantified errors. By Search 3, the priors are initialized such that a deeper exploration of the parameter space can be performed more efficiently, ensuring that dynesty does not spend considerable time in regions of parameter space that previous searches tell us do not give a physical lens model. Some basic settings that can be varied to affect the performance of the non-linear search are: (i) number of live points, (ii) number of steps of random walks per iteration, (iii) target acceptance fraction for random walks, (iv) Bayesian evidence tolerance, (v) positions threshold, and (vi) sub-grid size.

The technical details of our modelling method, including data preparation and pipeline, are described more fully for the interested reader in Appendices  A and  B. The sequence of chained searches and specific parameters that are set via prior passing are listed in Table B1, and dynesty settings are tabulated in Table B2. More details on pyautolens’s use of dynesty are provided in (Nightingale et al. in preparation).

4.2 Physically motivated priors

Where possible, we calculate priors using photometric observations from GAMA-DR3 preferentially over a universally applied ‘typical’ value. For either case, it is important not to fix the parameters too restrictively to values from observations.

4.3 Effective radius

In search 1, we initialize the effective radius parameter of the foreground lens light profile with a Gaussian centred at the lesser of two possible radii: (i) effective radius determined from photometric observations from GAMA-DR3 sersiccatSDSS v09 catalog (Kelvin et al. 2012), or (ii) the median SLACS lens effective radius and standard deviation from Auger et al. (2010a; 7 ± 3.3 kpc). We expect the GAMA-DR3 observation to include extended blended light from the source feature, which may result in a higher measured effective radius than would be measured from the foreground galaxy if it were not lensing. In order to assist the search in the task of deblending the lens and source light, we ensure that the prior is not predisposed to unusually large effective radii. Another failure state of early models resulted in unrealistically large source galaxy effective radii. Instead of attributing the extended lens features to lensing of a compact background object (which is most often the case for strong lensing), the model makes up for that extra flux as the physical extent of an extremely large, bright background source galaxy at high redshift. This motivated an upper limit to the effective radius of the source galaxy based on typical disc galaxy properties. We take a rough value of 7.5 ± 2.5 kpc and upper limit of about 15 kpc.

4.3.1 Lens mass-to-light ratio

In order to be implemented as priors in the lens models, these rest-frame constraints must be k-corrected, calibrated to the flux units in which the observed data are given and rewritten in the model mass and intensity units. We use SED-calculated k-corrections from GAMA-DR3 kcorr_auto_z00 v05 (Loveday et al. 2012) for each lensing galaxy to constrain the prior in the observed bandpass. These constraints are converted to angular mass units per eps (electrons per second) with the gain and exposure time of the KiDS observation. These constraints ensure that the model does not attribute mass to a stellar population that is impossible within current stellar evolution models. In cases where the maximum possible M_*/L is fit, the maximum possible stellar mass has also been attributed. In these cases, the model may end up having to compensate with very high amounts of dark matter to account for the lensing potential.

4.3.2 NFW profile scale radius

The scale radius of a dark matter halo is one of the key parameters of the NFW mass density profile for dark matter halos. Gavazzi et al. (2007) modelled 22 SLACS lenses with strong and weak lensing constraints and a two-component mass profile consisting of a de Vaucouleurs stellar component and spherical NFW dark matter component. We adopt their resulting mean scale radius of r_s = 58 ± 8h⁻¹ kpc for our dark matter profile Gaussian prior distributions. These values are converted to arcseconds from angular diameter distances in our assumed cosmology.

4.4 Choice of image bandpass

KiDS observations of each object in different bandpasses are not equal in exposure time, signal-to-noise quality, or PSF. KiDS r-band imaging is the highest quality of the bandpasses, with an exposure time of 1800 s and a mean PSF of 0.65 arcsec. g-band exposure times are 900 s. For each candidate and for each model search, we select the better of r- or g-band images. Search 1 fits the r-band image because it most clearly shows the foreground lens light. If the lens and source are clearly distinguishable in the r-band, then the same image is used for searches 2 and 3. However, the g-band image often most clearly shows the lensed features of the background source; in these cases, the g-band is preferable for search 2. Since search 3 models the entire system, the image that most clearly shows both profiles is used.

Each search assists the subsequent searches to distinguish between the lens and source light, which is one of the more difficult challenges of modelling lenses from images of the resolution and S/N attainable by ground-based observatories. Two effective first solutions are (i) separating the initial search of foreground lens and background source light profiles by colour-band and (ii) masking specific regions of the image. In this case, the sacrifice in quality between the r- and g-bands as a consequence of survey design presents an additional challenge to the fitting process as well as in later analysis. pyautolens uses units of electrons per second, so the effect of the difference in exposure times and S/N between observations with the two bandpasses is minimized. However, measurements taken from models that fit images of the same bandpass are much easier to compare.

5 MODEL QUALITY ASSESSMENT AND GRADING

With models complete, we assess the highest-likelihood models for each candidate. Ideally, a reliable objective figure of merit, such as image χ² or Bayesian evidence, would sufficiently quantify the quality of each lens model fit. Forming robust quantitative goodness-of-fit metrics is currently an open problem in automated lens modelling. Etherington et al. (2022) explored this problem using pyautolens with much higher resolution images and found that none did a particularly satisfactory job. We take the Bayesian evidence to be the reference figure of merit and follow this with a blind visual inspection of the image, fit, and spectrum of each modelled candidate. We inspect the images and spectrum separately in order to isolate some of the failure states that occur for each and have a clear picture of the factors limiting the precision of the models. Three collaborators give a separate score between 0 and 4 for each of the image, fit, and spectrum for each candidate, so that each of the candidates has three scores out of 12 and a total possible score of 36.

The procedure for assigning quality scores is as follows: The collaborator (the ‘scorer’) is shown via a randomized selection either the spectrum or a set of images (observed and model-fit) of a randomly selected candidate. The spectrum and image set are not shown sequentially in order to keep the scores unbiased by each. The spectrum score is based on the detection of redshifted line features that correspond to both the foreground and background redshifts. Wavelength accuracy, strength, and number of detectable line features are considered, in addition to template type and the presence of overlapping line features. The image and fit are scored simultaneously from the set of four observed and model-fit images because the fit score is informed by the image score. The image score is based on two images – (i) the observed image and (ii) the observed image with the model’s lens light Sérsic profile subtracted. The scorer considers how well the two images appear to show a well-defined structure outside the central foreground lens light-profile that could be reasonably described as a lens feature. For the fit score, the scorer compares the lens-subtracted model image to the lens-subtracted observed image and examines model background source-plane image. The fit score is influenced by the image score; the fit score cannot be higher than the image score +1. This means that a poor image that is fit perfectly should not get a high fit score, and a good image that is fit poorly should reflect the failure of the model to adequately attribute the image features to the lensing of a background source.

Following the scoring exercise, we remove any candidate that received a ‘0’ for any of the image, fit, or spectrum scores by any scorer. This removes catastrophic failures and ensures that the final set is reliable for follow-up analysis. The 19 candidate models that remain are assigned a letter grade of A, B, C, or D according to the structure outlined in Table 1. There are 2 A, 3 B, 9 C, and 5 D grades in the scored subsample, described in Table 2. 17 of the graded models are candidates from the LinKS subsample. Two of the three candidates from the ‘Li-BG’ sample were modelled to a level of success that justified presentation alongside the others, though both models are given grades of D. One D-grade model (G419067) with a negative likelihood was a result of high image residuals in the very centre of the lens light profile. Note the asterisk in the ln(evidence) column of Table 2. This model scored well enough for inclusion (spectrum score 4, total score 22) by the blind visual inspection. However, the visual inspection may have removed this candidate with the inclusion of a residual or χ² map in addition to the model images. The Bayesian evidence is therefore a prudent first cut of extremely poor models. Otherwise, as shown in Fig. 3, we find that the quality of fit determined by careful visual inspection is not correlated to the reported Bayesian evidence.

Figure 3.

Total score from visual quality inspection versus the natural log of the Bayesian evidence from model fitting. X-markers are rejected based on visual quality inspection. Square markers are accepted. There is very little correlation between the visual inspection results and the objective quality-of-fit metric.

Open in new tabDownload slide

To the authors’ knowledge, none of these lens candidates have been previously confirmed with high-resolution (HST-quality) imaging or spectroscopy. G250289 was identified in HSC as J083726 + 015639 by Sonnenfeld et al. (2019). Spectrum scores of six or better can be considered to be probable spectroscopic evidence for the lens candidate, and the highest spectrum scores for the two A grades can be considered spectroscopic confirmations. No additional extensive efforts were made to identify another possible background source redshift if the one determined by autoz was deemed unreliable. All scores can be considered useful follow-up evidence for the quality of the candidates, in that the success of a model lends additional confidence to the identifications. However, Petrillo et al. (2019a) and Li et al. (2020) have already provided extensive studies of the quality of their lens identifications, and our image modelling is conducted on the same observations as their analysis. Therefore, any poor model performance here does not contradict a positive identification.

6 MODEL RESULTS

6.1 Extracting best-fitting parameters

For each lens model, the parameter space is sampled over tens of thousands of iterations, estimating the log-likelihood for each sample fit and constructing a probability density function (PDF) for each free parameter listed in Table B1 of the appendix. Mass and light profiles are fully described by the model-fit parameters. The Einstein radius, total lensing ‘Einstein’ mass, mass fractions, luminosity, and mass-to-light ratios are calculated for each sample from the model grid and integrated over the angular area enclosed by the Einstein radius.

Parameters involving the luminosity require k-corrections to rest-frame (see Section 4.3.1). Three of the four models used the g-band image for Search 3, so they are easy to compare. Special attention should be given to G250289, which was instead modelled from its r-band image. In attempting to give the models the best chance to succeed, we were inconsistent in the choice of bandpass for modelling (see Section 4.4). Luminosity and mass-to-light ratios for this model are corrected to the g-band after r-band rest-frame k-correction. This is done by multiplying (or dividing) by an additional factor 10^{−0.4(g − r)}, where (g − r) ∼ 0.285 is the colour difference calculated by integrating the product of each bandpass response function and a template elliptical galaxy spectrum from Kinney et al. (1996) over the bandpass range. This has the effect of scaling luminosity down and M/L up. Given our goal here, which is to explore the methods, this is sufficient for characterizing the differences between models in a consistent parameter space. However, future efforts that intend to approach these measurements more rigorously should approach the initial modelling with more consistency.

We briefly present best-estimate results from the highest-likelihood model fits for the four highest quality models in Table 3, selected primarily by the blind quality scoring described in Section 5. Bayesian evidence reported by pyautolens and the subjective reasonableness of the inferred quantities are also considered in the selection of this small subsample. One B-grade model, G62734, is not included because its central dark matter content is poorly constrained. This and the other 14 lower-grade models are considered to be worth revisiting but were not successful enough to present alongside the cleaner examples we present here.

To discuss inferred quantities, we estimate bivariate PDFs for the quantities using a Gaussian kernel-density estimate from the final 10 000 iterations. The best estimate for each inferred quantity listed in Table 3 is determined at the maximum of one of these bivariate PDFs, where we use uncorrelated values as much as possible. We show the observed image, maximum-likelihood model fit, and spectrum for the highest scoring model in Fig. 4. The other three models listed in Table 3, as well as G62734, are shown and discussed in more detail in Figs C1–C4 of Appendix  C.

Figure 4.

G323152. A-Grade. Upper left: The observed image shows an apparent arc feature above the central lens galaxy light-profile. Upper right: The model image captures the extra light reasonably, though without the exact shape. This could be a result of internal substructure or the impact of shear along the line of sight, both of which are unaccounted for in the model. Lower: The GAMA spectrum shows strong line features at the redshifts at 0.353 and 0.722 identified by autoz. Dotted lines identify foreground lens galaxy absorption features (H, K, Hβ, Mg, and Na) at z = 0.353, and dashed lines show background source emission features (Hβ, [O ii], [O iii]) at z = 0.722.

Open in new tabDownload slide

We are primarily interested in studying the stellar and dark matter content in the central regions of the lens galaxies. Although the mass and light profiles in the models are inferred to a larger radius, mass measurements via strong lensing are the most precise when considering only mass within the Einstein radius of the galaxy. It is important to note that the Einstein radius is a feature individual to each system, so the quantities are not calculated within a uniform radius from the centre of each galaxy. Values for each Einstein radius can be found in Table 3.

6.2 Comparing highest-quality model results

We show the four highest-quality models for comparison. Figs 5–7 show the four models in parameter space of interest to our study. All of the plotted quantities are taken within the Einstein radius (see Table 3). Each model is identified with a different marker, and B-grade group-member galaxy G250289 is indicated with a red marker to remind the reader that the final model fit utilized the r-band. Green and blue contours enclose 1σ (39 per cent) and 2σ (86 per cent) of the two-dimensional PDF, respectively. With the small sample size shown here, we do not intend to address questions of assembly bias and galaxy formation mechanisms. These plots are intended to discuss the cleanest subset of our sample in the context of what can be considered more thoroughly in future work.

Figure 5.

Mass components integrated within the Einstein radius of each of the four best-fitting models described in Section 6 and Table 3. The legend shows the GAMA identifier, quality grade, environment classification, and the SDSS bandpass used for the final model fitting. Green and blue contours about each point enclose 1σ (39 per cent) and 2σ (86 per cent) of the two-dimensional PDF, respectively.

Open in new tabDownload slide

Fig. 5 shows the integrated stellar mass and dark mass enclosed within the Einstein radius of the lens models. These are obtained by integrating over the Sérsic stellar mass profile and elliptical NFW profile. The galaxies have total enclosed Einstein mass values of order M_E ∼ 3−8 × 10¹¹ M_⊙, which is expected since lensing galaxies are typically quite massive. Assembly bias would show itself here as a trend where group central galaxies tend to have higher stellar mass than isolated galaxies at the same dark matter halo mass. The only group-member galaxy, G250289 (marked with a red cross), lies in one of the smaller dark matter haloes and has the highest stellar mass. G323152 (marked with a black circle, not listed in GAMA GroupFinding catalog) has a similar dark mass and about one-fifth of the stellar mass compared to G250289. Conversely, the two isolated galaxies have the highest dark masses and the lowest stellar masses. The higher stellar mass in G250289 could have more to do with effects from the difference in r-band and g-band S/N than a physical interpretation. With so few data, it is difficult to determine how much of an effect this difference has on the results of the models.

The total lensing (Einstein) mass enclosed within the Einstein radius is generally well constrained. We want our models to further constrain the fraction of this lensing mass that is dark matter. The fraction of dark matter is not directly constrained by a model prior but is very sensitive to the assumed forms of mass and light profiles and the constraints placed upon those. Fig. 6 shows the g-band integrated luminosity and dark matter fraction (both calculated within the Einstein radius) for each model. The uncertainties in the fraction of dark matter along the x-axis are relatively small, which is surprising given the inherent degeneracy of stellar and dark mass in lens modelling. The small parameter space explored could indicate a lack of flexibility of the models’ stellar and dark mass profile priors, perhaps an excessive constraint or weighting towards one mass component over the other. As discussed in Section 4.2, the careful selection of priors can be challenging. Compare again G250289 (red cross) and G323152 (black circle), which have similar enclosed dark masses. Even corrected (scaled down) from r-band to g-band, the enclosed luminosity for G250289 is 2–4 times greater than the other three, which could be why the model attributes a higher fraction of the Einstein mass to the stellar component. These models have relatively similar total Einstein masses enclosed within similar Einstein radii around 1.2–1.7 arcsec. G250289 has an Einstein radius of ∼1.5 arcsec that is typical and within the range of the other model values, so additional luminosity and stellar mass is not a result of an overextended radius of integration.

Figure 6.

SDSS g-band luminosity and dark matter fraction integrated within the Einstein radius of each of the four best-fitting models described in Section 6 and Table 3. Legend and marker information are the same as in Fig. 5.

Open in new tabDownload slide

Fig. 7 shows the g-band stellar mass-to-light ratio compared to the enclosed dark mass. Dotted lines at the 2σ contours indicate upper constraints placed on the mass-to-light ratio. This figure completes our discussion of the degeneracy. In summary, the model has two options for attributing the lensing mass: (i) To favour the stellar component, a higher stellar mass can be the result of a heavier stellar population, and (ii) conversely, a very large, centrally concentrated dark matter halo can make up for a lower stellar mass and luminosity. This is all expected. Bounds of integration for luminosity, stellar mass, and dark mass are all dependent on the measure of the Einstein radius, which is in turn dependent on the total mass which the model is attempting to parse into stellar and dark components. It is a complex problem with degenerate variables that is only constrained by meaningful priors. In our cleanest subsample, the most identifiable differences occur for the candidate that was modelled from a different photometric bandpass. This is another experimental design decision based on limitations of the data that significantly affected our ability to analyse the resulting models. Thus, even our best models suffer.

Figure 7.

Stellar mass-to-light ratio (M_*/L_g) in g-band and dark mass integrated within the Einstein radius of each of the four best-fitting models described in Section 6 and Table 3. Dashed lines at the 2σ contours indicate boundaries corresponding to upper constraints placed on the mass-to-light ratio of the models (see Section 4.2). Legend and marker information are the same as in Fig. 5.

Open in new tabDownload slide

7 AUTOZ CONSIDERATIONS POST-MODELLING

Here, we summarize the results of our quality control procedure described in Section 5 and give several recommendations for removing contaminants. Appendix  D gives a more detailed breakdown of the candidates and models in the context of the cases discussed in Section 3.2.

We discovered some failure modes for our method of utilizing the autoz algorithm for background source redshift determination. Revisions to our initial procedure removed about half of the candidates. Few single cases should be excluded without consideration; for example, some redshift determinations were kept even though there appeared to potentially be a falsely attributed line. However, absorption and emission lines must both be checked for overlap regardless of the template configuration type, and template type should be questioned with this information. Several of the lens ELG galaxies turned out to have prominent absorption features at the same redshift. The binary identification of PG and ELG, as we have done here, is an oversimplification that hinders both our classification and understanding of the spectra and expected galaxy properties.

The two A-grade candidate models were acquired with autoz redshift matches that we consider to be cases that deserve extra care (where the foreground galaxy is best matched with an ELG template). This is interesting because one may be tempted to cut these cases entirely in order to obtain a clean and fairly homogeneous sample (i.e. large elliptical galaxy lensing a bright ELG, PG + ELG). We find the more careful consideration of other possible cases to be fruitful. If we had adopted a selection that focused only on configurations where the primary template matched a passive lensing galaxy, then |$\sim 20{{\ \rm per\ cent}}$| of our final selection, including the two highest-scoring models, would have never been considered. Care should be taken in order to reap the benefits of this expanded population while minimizing contamination.

In the big picture, the improvement of spectral quality in wide-field surveys is essential for making this work in an automated way over large sample sizes, but an automated redshift algorithm like autoz could be optimized for background source redshift determination. Our subjective quality scores show little correlation to the autoz output parameters that we used for the initial selection, as shown in Fig. 8. The two axes show the autoz selection parameter space, composed of the second cross-correlation peak σ₂ and the R parameter. Three of the candidates with the highest σ₂ show very poor spectrum scores because they are instances of overlapping emission lines. This suggests that a higher threshold for σ₂ may not actually yield higher-quality background source redshift matches. We now introduce other options for maintaining a clean sample selection with autoz while expanding the sample size in future works.

Figure 8.

Upper: Spectrum quality scores. Lower: Total score shown as scatter plot colour variations scaled with the colourbar on the right of each plot. The axes of the scatter plot are the selection criteria from Fig. 1. Vertical axes are the second-highest cross-correlation peak σ₂, and horizontal axes are the R parameter. Dark purple colours indicate poor scores, with higher scores at orange and yellow. Little correlation can be seen between these parameters and the results of the models or their subjective spectrum scores. In the upper plot, four low-scoring spectra with high σ₂ correspond to overlapping emission lines.

Open in new tabDownload slide

7.1 Recommendations to remove contaminants from AUTOZ selection

One could quite easily remove contaminants during the automated selection. Each of these recommendations pays particular attention to the redshifts of ELGs in both the foreground lens and background source positions. The two most convincing (A-grade) spectra and models were cases of ELG+, so we want to remove as many contaminants as possible without the blanket removal of either of these cases. In order to maintain the applicability of this procedure to an even larger set of data than is considered here, we recommend the following selection criteria be implemented when adopting automated redshift determinations:

• Remove ELG+ELG and PG + PG matches where emission or absorption lines redshift to overlapping wavelengths between foreground lens and background source. One can calculate lens and source redshift combinations that result in overlapping observed emission lines similarly to the procedure described in Holwerda et al. (2015). The following equation defines a region of parameter space between a lower and upper linear function of (1 + z_source) to (1 + z_lens). Within this region, an overlap will occur for a given pair of rest-frame emission line wavelengths:

  $$\begin{eqnarray*} \left| \frac{1+z_{\mathrm{source}}}{1+z_{\mathrm{lens}}} - \frac{\lambda _{r,\mathrm{lens}}}{\lambda _{r,\mathrm{source}}} \right| \lt \frac{A}{R} \end{eqnarray*}$$

  (4)

  where z_source and z_lens are the source and lens redshifts, λ_{r, lens} and λ_{r, source} are the rest-frame wavelengths of emission lines from lens and source, R is the spectral resolution of the instrument, and A is a coefficient that widens the range of exclusion for potential overlapping features. The equation implicitly accounts for the dependence of resolution on observed wavelength. Fig. 9 shows the regions where a redshift combination of foreground lens and background source results in overlapping lines given GAMA’s spectral resolution of R ∼ 1300 and A = 4 (i.e. overlapping emission lines are closer than four times the smallest resolvable wavelength difference). This prescription identifies most of the cases of overlap that were flagged by direct visual inspection of the spectra. We retained one of them with the lowest possible D-grade because its SDSS-BOSS spectrum showed fairly reasonable source Hβ and [O iii] couplet lines at higher-wavelength that were not in the range of the GAMA spectrum.

• Remove + ELG configurations where Hβ and [O iii] couplet emission lines are redshifted beyond the wavelength range of the observation. Table 4 and Fig. 10 show the redshift limits beyond which Hβ and [O iii] λλ4959,5007 lines will be above the survey upper wavelength limit for several optical spectroscopic surveys. Spectroscopy in the 1-|$\mu$|m range is necessary in order to detect the emission lines from sources at z ∼ 1 and will become more important for modelling lenses with foreground lens redshifts higher than z ∼ 0.5. DESI and the upcoming 4MOST have higher resolution and cover extended optical wavelength ranges that correspond to these redshifts, which will reduce the significance of this problem. Background source redshifts could also be assessed by looking for emission lines in very near-infrared (e.g. MOSFIRE Y-band, 0.97–1.12 |$\mu$|m) spectra, where possible. In principle, template-based automated redshift identification could be run for each separately, which would negate some of the difficulties inherent to identifying the two distinct signals within the single observation. The obvious negative to this option is the need for additional observations.

• Remove ELG+ matches with any of the following characteristics: (a) low foreground lens redshift (less than z ∼ 0.2 for our sample), (b) ELG + PG configuration, and (c) primary redshift match to the background source. ELG+ configurations with low foreground lens redshifts suffered from several failure conditions in addition to being a less-likely configuration than PG+. For GAMA, the noisy short-wavelength end of the observed spectral range corresponds to emission lines with z_lens less than ∼0.2, leading to mistaken classification of noise peaks as emission lines. While this is specific in part to GAMA’s wavelength-specific spectral performance, several of these low-redshift foreground lens matches are also ELG+PG configurations, of which almost all were removed in quality scoring. The one remaining candidate had the lowest score of all those accepted. Perhaps even more condemning is the fact that many of these low-redshift ELG+ foreground lens matches were also cases where the background source was assigned the primary redshift. This trio of doubtful cases coincided for several of the candidates that were removed from our sample during quality scoring.

Figure 9.

Combinations of foregrounds lens and background source redshift that result in overlapping important emission or absorption line features. Coloured line-regions are calculated with equation (4) and indicate parameter space where the labelled line features will overlap. Each is labelled in the legend as ‘source feature - lens feature’. Using these functions identifies 19 of the 20 overlaps in the 42 candidates and all 6 that made the final selection of 19. The region in the lower right between the solid and dashed black lines shows the selection criterion utilized in the initial autoz selection, which removed candidates where the source redshift was within 0.1 of the lens redshift.

Open in new tabDownload slide

Figure 10.

Blue, green, and purple lines show the redshift of rest-frame Hβ and [O iii] λλ4959,5007. Dotted vertical lines show the upper wavelength limit of various spectroscopic surveys, and the overlaid arrows point to the redshift at which the first of these three lines will disappear in the survey observations. The red shaded region indicates the very near-IR coverage of MOSFIRE Y-band (0.97–1.12 |$\mu$|m) that could also potentially reveal background source emission lines at around z ∼ 1 and greater.

Open in new tabDownload slide

8 DISCUSSION

8.1 GAMA environment

The sample of 42 KiDS lens candidates and the subsample of 19 with accepted grade A–D models are both close to evenly split between group-member and isolated galaxies according to the GAMA groupfinding metrics. 22 (8 graded) are associated with groups and 19 (9 graded) are isolated. (G323152) is not represented in GAMA group catalogs.

We note that most strong lensing galaxies should be the most massive galaxies in their halo, either in a group or in isolation. Fig. 11 shows the rank of the proximity of the object to the centre of mass of the group relative to other group members, with one indicating the closest or centre-most galaxy. Most of the candidates shown here are group central galaxies, but our models failed for 11 of those. On the other hand, four of the five candidates that are not the central galaxy of their group were accepted and given grades, comprising 40 per cent of the graded group-galaxy members. Low scores for group member galaxies could mean that their identification as lenses is a false-positive given their proximity to other significantly large galaxies. Alternatively, if these are lenses, the modelled mass structure of the lensing galaxy as a single light-mass component and assumption of its presence in the centre of the dark matter halo may not be accurate enough to reproduce the lensing observables. If this were the case, one might expect the group centrals to be more easily modelled than the subdominant or ‘satellite’ galaxies with ranks of 2 or greater, which is not apparent in Fig. 11. The two B-grade group galaxies are ranked 1 and 2, indicating that one of them is a central while the other has at least one companion that is competing for dominance in the group. The rank 2 group member is G62734, which was removed from the final selection because its dark matter content was poorly constrained. This could be a result of this galaxy’s distance from the centre of the group mass.

Figure 11.

Stacked histogram. Quality scores for 41 of 42 candidates in reference to data from GAMA groupfinding catalogs (one candidate does not have environment data). Location on the x-axis distinguishes ‘Isolated’ from group member galaxies, which are further separated by the rank in projected distance from the centre of mass of the group. A rank of 1 indicates that the lens candidate is the central galaxy of the associated group. Colours indicate the quality grade of models, with no colour indicating models that were not accepted.

Open in new tabDownload slide

Compared with the SLACS study in Treu et al. (2009), in which 12 of 70 (17 per cent) were associated with groups, our KiDS/GAMA strong lens candidate sample and selected subsample of 19 models is more highly represented by group-member galaxies. Definitions of group membership based on environmental parameters are not the same between these studies. The nearly 50/50 split between group-member and isolated galaxies in our sample does not necessarily support or dispute a preference for overdense environments by lensing (and all massive) elliptical galaxies. However, the high completeness of GAMA compared with SDSS may instead suggest that our sample minimizes the apparent environmental preference. The distinction of group association here could be affected by selection bias, as those designated as isolated could in fact be groups with satellite members beyond the GAMA flux limit. If this were the case, there should be a systematic bias in isolated galaxies towards higher redshift. Fig. 12 shows that neither subsample of group member or isolated candidates is significantly distinguished in redshift or stellar mass.

Figure 12.

Redshifts and stellar masses of autoz sample from GAMA-DR3 stellarmasseslambdar v20 (Taylor et al. 2011) determined by stellar population and separated by group-member and isolated galaxies according to GAMA team internal groupfinding catalogs (Robotham et al. 2011). There is no clear distinction between the subsamples in either observable.

Open in new tabDownload slide

With more data and better measurements than we have accomplished here, one may be able to compare observations to the scatter in the upper plot of fig. 9 of Zehavi et al. (2018), where for fixed dark halo mass, higher stellar-mass galaxies tend to exist in denser environments. Note that the modelled mass components here are calculated within the Einstein radius and not the full extent of the galaxy. The majority of the dark halo component should extend well beyond the stellar halo, and these high-mass lensing galaxies are more likely to exist in more massive dark matter haloes (⁠|$\mathrm{log}(M_h/h^{-1}\, \mathrm{M}_{\odot })\sim 12\!-\!14$|⁠) where the suggested environmental trend is less supported. The precision and numbers required to test assembly bias will require more refinement of the methods discussed in this study as well as the power of more sophisticated surveys to come.

8.2 Future work: a place for ground-based observations

Realistic lens modelling by fitting mass and light profile parameters is a complex problem with a large number of parameters. With even the highest-quality ground-based imaging offered by the likes of the KiDS, the angular resolution is insufficient to constrain the individual model solutions to levels where one can make strong inferences about the individual lens galaxies. There are simply too many solutions that fit the image to a high probability, which inflates the uncertainty to levels that make it difficult for one to draw conclusions from the inferred quantities. These uncertainties on a single lens can be significantly constrained with the level of imaging afforded by adaptive optics (AO) or space-based instruments. Fig. 13 shows a model solution for one of the lens models after being simulated with the optics for three observatories: (i) VST used for KiDS, which was the instrument that collected the original image, (ii) LSST at the Vera Rubin Observatory (VRO) representing the next generation of ground-based observatories, and (iii) Advanced Camera for Surveys on the Hubble Space Telescope. The same model-fitting procedure applied to HST images or observations with AO of the same lensing galaxies would result in error estimates an order of magnitude better than the results we achieve here. Alternatively, future systematic modelling of orders of magnitude more ground-based, lower-resolution observations (as we expect to achieve with observatories like the VRO) can result in similar precision. Constraints at the population level (made possible with these larger sample sizes) can enhance higher-resolution individual measurements through Bayesian hierarchical frameworks. Our work demonstrates the value of wide-field, lower-resolution surveys as a complementary tool to the expensive and hyper-competitive observing campaigns that are the default for strong lens studies.

Figure 13.

Upper left: Maximum log-likelihood model for G3629152 shown with pixel-scale 0.2 arcsec pixel⁻¹. Other images are simulated with identical background sky and convolved with optics of upper right: VST (r-band PSF 0.65 arcsec, pixel scale 0.2), lower left: LSST at Rubin (PSF 0.5, pixel scale 0.2), and lower right: ACS on HST (PSF 0.1, pixel scale 0.05). Imaging from space-based observatories or AO would allow for better model-fitting and tighter uncertainties for future efforts.

Open in new tabDownload slide

The next generation of spectroscopic surveys is already underway, e.g. the DEVILS deep survey on the AAT (Davies et al. 2018; Holwerda et al. 2021) and the DESI redshift survey. One can expect increased numbers of spectroscopic lensing candidates as well as opportunities to identify the redshift of the potential background source. More comprehensive spectroscopic surveys are being planned with the 4MOST instrument (de Jong et al. 2012; Depagne & 4MOST consortium 2015). These planned surveys include extra-galactic ones such as the two-tiered Wide Area Vista Extragalactic Survey (Driver et al. 2019), the Optical, Radio Continuum and HI Deep Spectroscopic Survey (Duncan et al. in preparation), and a cosmological low-S/N wide-area survey (Richard et al. 2019). These 4MOST surveys are expected to achieve high completeness in their target fields and yield a boon of spectroscopically confirmed strong lensing systems with the same advantages exploited by the procedure outlined here at better spectral resolution and wider fields.

Identifications of strong gravitational lenses through imaging are also expected to increase in the near future with observations by the likes of the VRO, Euclid, and the Roman Space Telescope, in addition to improved machine learning techniques. Following the discussion in Knabel et al. (2020), one expects the selection functions of the spectroscopic surveys and these optical and near-infrared imaging surveys to show a limited improvement in overlap. The analysis presented here, however, shows that in the overlap a useful subset of strong lenses can be utilized for modelling by imaging and spectroscopy combined.

The lenses discussed here, and in future similar ground-based efforts, are also ideal candidates for deeper follow-up observation with higher-resolution imaging and spectroscopy. These follow-up observations could include Integral Field Unit observations to measure the stellar population characteristics across the elliptical, chart the foreground lens galaxy kinematics, as well as study the background source light and stellar population characteristics. These considerations are more aligned with and have been sufficiently described in the existing literature and will not be discussed further here.

9 CONCLUSIONS

We arrive at the following conclusions from our analysis of strong lens candidates in GAMA and KiDS using autoz and pyautolens:

• Meaningful strong-lens studies can be conducted by applying lens-modelling methods such as those we have outlined here to large imaging and spectroscopic surveys.

• Automated template-matching redshift algorithms like autoz can be utilized to determine reliable background source redshifts required for lens modelling. Careful consideration should be taken in cleaning the algorithm’s selection, following the recommendations outlined in Section 7.1.

• Limits of optical resolution in large ground-based surveys present significant challenges to the uniqueness of solutions in our Bayesian modelling of individual strong lenses.

• As sample sizes grow, refinements to these techniques can produce lensing measurements in quantities that will offer considerable statistical power. This approach is complementary to the more detailed modelling of individual lenses that is possible with deeper and higher resolution observations.

ACKNOWLEDGEMENTS

SK acknowledges NASA Kentucky, National Science Foundation (NSF), University of Louisville, and University of California, Los Angeles, for financial and technical support. The material of this work is based upon work supported by NASA Kentucky under NASA award No. 80NSSC20M0047. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. 2021325146. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

AHW is supported by an ERC Consolidator Grant (No. 770935).

This work uses the following software packages: astropy (Astropy Collaboration 2013; Price-Whelan et al. 2018), colossus (Diemer 2018), corner.py (Foreman-Mackey 2016), dynesty (Speagle 2020), matplotlib (Hunter 2007), numba (Lam, Pitrou & Seibert 2015), numpy (van der Walt, Colbert & Varoquaux 2011), pyautofit (Nightingale, Hayes & Griffiths 2021a), pyautolens (Nightingale & Dye 2015; Nightingale et al. 2018, 2021b), python (Van Rossum & Drake 2009), scikit-image (Van der Walt et al. 2014), scikit-learn (Pedregosa et al. 2011), and scipy (Virtanen et al. 2020).

GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/ .

Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples).

Based on data from the STELIB service developed by the Spanish Virtual Observatory in the framework of the IAU Comission G5 Working Group: Spectral Stellar Libraries.

DATA AVAILABILITY

KiDS images and data used in this paper are available from the Astro-WISE Data base Viewer Web Service (dbview.astro-wise.org). LinKS-specific data can be found at the LinKS website (https://www.astro.rug.nl/lensesinkids/). The GAMA autoz catalog is available from GAMA-DR3 website (http://www.gama-survey.org/dr3/, aatspecautozall, specall, lamdarstellarmasses, sersiccatsdsS, and kcorr_auto_z00 catalogs) and the team internal GroupFinding catalog for the full GAMA fields will be made available in GAMA-DR4 (Driver et al. in preparation).

Footnotes

REFERENCES

APPENDIX A: PREPARING DATA FOR MODELLING

Images and weight maps are 101 × 101 pixel (∼20 × 20 arcsec²) cutouts from coadded images of KiDS tile observations acquired from the publicly available Astro-WISE Data base Viewer Web Service.⁴ g- and r-band images are cut out centred on the object’s RA and Dec., recentred to the brightest pixel in the central (lens) galaxy light profile, and converted to eps (electrons per second) for modelling. KiDS image pixel values in the Astro-WISE Data base are given in calibrated flux units relative to the flux corresponding to magnitude 0 and are converted to ‘brightness’ units of electron counts by multiplying by the tile’s average gain, which includes additional factors necessary for this conversion. pyautolens is by default set to be optimally utilized with units of electrons per second (eps), which is acquired by dividing by the exposure time (1800 s for r-band, 900 for g-band).

pyautolens requires input of the PSF and noise map for each image. The inverse square root of the weight map corresponding to the cutout image gives the rms noise, which is converted to electron counts and squared to recreate the background sky. We then add this image to the corresponding cutout image and take the square root to give the noise map, after which we convert to eps. We generate a Gaussian PSF for each image from the average full width at half maximum (FWHM) of the PSF for each image.

We next mark pixel-positions of the distorted images of the lensed background source in each image, when visible, using a GUI distributed with pyautolens. During a lens model fit, pyautolens casts aside all mass models where these image-pixel positions do not trace within a designated threshold of one another in the source plane. This narrows the parameter space that is searched and ensures that the model fits the observed image features of interest.

We generate three masks for the three searches with each candidate: (i) lens mask – a circular aperture tailored to show only the lens galaxy (on the order of but usually slightly less the effective radius, typically around 1–1.3 arcsec); (ii) source mask – a circular annular aperture showing only the light we determine to be the lensed background source features (with inner radius about the size of the circular lens mask and outer radius around 3 arcsec); and (iii) full mask – a circular aperture of typically around 3 arcsec that includes most of the light from the lens and source features and masks as many peripheral contaminants as possible.

APPENDIX B: LENS MODELLING PIPELINE

This section details continue the description of our lens modelling methods summarized in Section 4.1. Through experimentation, we have designed a pipeline composed of a chain of three dynesty searches that we use as a template for fitting each lens. The variety of lensing configurations, image quality, etc. force us to tailor aspects of each model-fit individually, in particular alternating the masks that segment the foreground lens light and distorted background source light. In order to institute the least bias possible, we allow the models to probe a wide range of possible solutions for each parameter. The shape of the prior distribution has significant effects on the performance of the search. We use uniform, log uniform, and Gaussian functions depending on the parameter and informative auxiliary observations. In the following sections, we describe this three-step automated pipeline, where from here on we refer to a ‘search’ as a model-fit performed by the non-linear search dynesty. Each subsequent search in the chain has more complexity in the form of additional parameters, which we balance in computational time by passing priors from previous search outputs. For each non-linear search in the chain, the priors are described in Table B1, and the dynesty settings are given in Table B2.

B1 Search 1 – lens light

For the r-band images, many of the lensed background source’s features are positioned within the lens galaxy’s effective radius. Search 1 therefore struggles to deblend the lens and source light, and the distorted arcs of background source light are attributed to the foreground lens galaxy. In some cases, this leads to a model solution that describes a lens light profile that is very large and very elliptical. To mitigate this systematic effect, we use the aforementioned lens mask for search 1 and constraints on the effective radius to assist the search to focus on fitting the lens light and not the source light. The residuals and uncertainties of this search therefore tend to be quite high. In fact, the residuals often outline the lensed images of the source itself, and the resulting maximum log-likelihood is lower than the value inferred in the second and third searches.

B2 Search 2 – Lens mass and source light

Search 2 focuses on the light from the background source. This component is modelled as a spherical exponential light profile in the source plane defined at the source redshift. The exponential light profile corresponds to the simple n = 1 case of equation (B1) and is parametrized using its (source-plane) centre, effective radius, and intensity. With higher-resolution imaging data, the background source light profile could be fit with a more detailed model. The background source centre coordinates are initialized to values within 2 arcsec of the line of sight of the foreground lens centre. The intensity of the background source light is again set to a wide log-uniform prior distribution, as for the lens light in Search 1. The source effective radius is initialized with a Gaussian distribution around a typical disc galaxy size, as discussed in Section 4.2. To map coordinates to the source plane, the lens galaxy’s total mass is modelled as a singular isothermal elliptical (SIE) profile. The mass profile’s Einstein radius prior is a wide Gaussian distribution centred at 1.0 arcsec with a hard upper limit of 2.5 arcsec. The centre and elliptical components of the SIE are paired with the light profile with the assumption that the ellipses will be aligned. The lens light profile takes prior distributions passed from the results of search 1. By fixing the centre of the lens profiles to the results of search 1, the lens model in this search is reduced to 10 free parameters. This, in addition to taking informed Gaussian priors from search 1 for the lens light, helps the model to focus on solutions that fit the source-light instead of systematic solutions that fit artefacts in the data. We use the annular source mask that removes the lens light from the observed image and therefore further focuses the model on fitting the source light. We also utilize pyautolens’s position resampling functionality, whereby the brightest pixels in the lensed source are marked (via a GUI). Again, the results of this search are passed as priors to Search 3.

B3 Search 3 – Combined lens and source models

Search 3 fits every component of the system. To model the foreground lens, we use a combined elliptical Sérsic mass-light profile for the stellar component and an elliptical NFW profile for the dark matter halo. Background source light is modelled again as a spherical exponential profile. Priors are passed from search 2 (see Table B1), except for the dark matter halo profile and stellar mass-to-light ratio. The prior distributions for these crucial parameters are determined by calculating central and limiting values according to a more careful process, as described in Section 4.2. The full mask including all the lens and source features is used to remove background features and other contaminants that exist in the periphery, which saves computational time. Lensed image positions are again used to discard unphysical mass models. This final search produces a reasonable fit to the complexities introduced by each component and gives uncertainties on each of the inferred quantities. Additional disc and bulge components, cores, and multiple galaxies at different planes along the line of sight can be fitted and would allow more precise and realistic models. These improvements to realism are unhelpful here given the quality of imaging available for the objects in question but would be simply applied in future studies following the same principles outlined in our strategy.