MIFAL: fully automated Multiple-Image Finder ALgorithm for strong-lens modelling -- proof of concept

TL;DR

This paper presents MIFAL, an automated method for identifying multiple images in strong lensing clusters, combining arc detection, photometric redshifts, and initial lens models to improve mass distribution analysis.

Contribution

The paper introduces a fully automated procedure for finding multiple images in strong lensing clusters, reducing reliance on manual identification and enabling large-scale analyses.

Findings

Successfully recovered all known systems in the first cluster

Recovered about half of the known systems in the second and third clusters

Identified some false images depending on parameter thresholds

Abstract

We outline a simple procedure designed for \emph{automatically} finding sets of multiple images in strong lensing (SL) clusters. We show that by combining (a) an arc-finding (or source extracting) program, (b) photometric redshift measurements, and (c) a preliminary light-traces-mass lens model, multiple-image systems can be identified in a fully automated (`blind') manner. The presented procedure yields an assessment of the likelihood of each arc to belong to one of the multiple-image systems, as well as the preferred redshift for the different systems. These could be then used to automatically constrain and refine the initial lens model for an accurate mass distribution. We apply this procedure to \emph{Cluster Lensing And Supernova with Hubble} observations of three galaxy clusters, MACS J0329.6-0211, MACS J1720.2+3536, and MACS J1931.8-2635, comparing the results to published SL analyses where multiple images were verified by eye on a particular basis. In the first cluster all originally identified systems are recovered by the automated procedure, and in the second and third clusters about half are recovered. Other known systems are not picked up, in part due to a crude choice of parameters, ambiguous photometric redshifts, or inaccuracy of the initial lens model. On top of real systems recovered, some false images are also mistakenly identified by the procedure, depending on the thresholds used. While further improvements to the procedure and a more thorough scrutinisation of its performance are warranted, the work constitutes another important step toward fully automatising SL analyses for studying mass distributions of large cluster samples.

Tables

Table 1. Table 1: MIFAL results for M0329

Id	R.A.	Dec	$z_{ima}$	$z_{sys}$	$z_{norm}$	$R_{find}$	${\chi }_{col}^2$	${\chi }_z^2$	${\chi }_{R_{find}}^2$	${\chi }_{SB}^2$	${\chi }_{ima}^2$	${\chi }_{sys}^2$	Comment
	${(}^{\circ })$	${(}^{\circ })$				${(}^{\prime \prime })$
1.1	52.424900	-2.187690	3.25	3.26	1.46	2.3	0.60	0.00	2.74	2.30	5.63	23.74	Sys. c5 in Z15
1.2	52.415100	-2.196250	3.26	3.26	1.46	2.3	13.08	0.00	2.69	0.00	15.78	23.74	”
1.3	52.416100	-2.194830	3.26	3.26	1.46	0.4	2.24	0.00	0.08	0.00	2.32	23.74	”
2.1	52.412700	-2.197070	2.87	2.83	1.46	1.6	2.16	0.02	1.33	1.71	5.23	40.05	Sys. 3 in Z15
2.2	52.422000	-2.187520	2.79	2.83	1.46	1.9	17.78	0.11	1.92	8.93	28.74	40.05	”
2.3	52.417500	-2.190630	2.83	2.83	1.46	2.3	2.91	0.03	2.67	0.47	6.07	40.05	”
3.1	52.425725	-2.190436	2.13	2.19	1.46	3.0	0.23	0.34	4.49	0.25	5.31	34.48	Sys. 2 in Z15
3.2	52.426393	-2.198453	2.24	2.19	1.46	2.3	3.83	1.76	2.66	0.19	8.44	34.48	”
3.3	52.415078	-2.200176	2.26	2.19	1.46	2.7	1.87	3.02	3.66	1.12	9.68	34.48	”
3.4	52.421065	-2.191334	2.13	2.19	1.46	1.8	1.65	0.34	1.65	7.41	11.05	34.48	”
4.1	52.416600	-2.204170	3.20	3.06	1.46	4.7	14.00	3.02	11.39	16.13	44.54	51.35	False detection
4.2	52.430903	-2.199237	2.93	3.06	1.46	0.2	5.06	0.01	0.02	1.72	6.81	51.35	”
5.1	52.417500	-2.195560	5.79	6.06	1.46	1.9	1.97	4.16	1.89	10.41	18.42	48.62	Sys. 1 in Z15
5.2	52.429900	-2.188120	6.08	6.06	1.46	3.9	0.66	0.01	7.87	1.43	9.98	48.62	”
5.3	52.421900	-2.201400	6.05	6.06	1.46	2.5	0.89	0.01	3.12	0.47	4.49	48.62	”
5.4	52.416800	-2.198060	6.31	6.06	1.46	4.9	0.14	2.66	12.39	0.55	15.73	48.62	”

Table 2. Table 2: MIFAL results for M1720

Id	R.A.	Dec	$z_{ima}$	$z_{sys}$	$z_{norm}$	$R_{find}$	${\chi }_{col}^2$	${\chi }_z^2$	${\chi }_{R_{find}}^2$	${\chi }_{SB}^2$	${\chi }_{ima}^2$	${\chi }_{sys}^2$	Comment
	${(}^{\circ })$	${(}^{\circ })$				${(}^{\prime \prime })$
1.1	260.067820	35.601192	2.47	2.39	3.85	5.9	22.94	0.65	17.69	0.39	41.67	65.64	False identification
1.2	260.062460	35.612745	2.31	2.39	3.85	0.2	22.94	0.38	0.02	0.63	23.97	65.64	”
2.1	260.076530	35.609749	3.22	2.88	3.85	3.0	0.47	0.85	4.51	0.18	6.01	79.88	Sys. 1 in Z15
2.2	260.073320	35.612099	3.11	2.88	3.85	4.4	1.17	0.16	10.06	0.10	11.49	79.88	”
2.3	260.059560	35.611026	2.66	2.88	3.85	5.7	3.06	3.54	16.85	0.31	23.75	79.88	Likely
2.4	260.066420	35.604602	2.65	2.88	3.85	7.7	4.68	3.79	30.07	0.10	38.64	79.88	Uncertain
3.1	260.078200	35.604918	1.68	1.68	3.85	7.7	12.69	0.50	30.26	3.68	47.13	50.84	Sys. 3 in Z15
3.2	260.067580	35.609483	1.68	1.68	3.85	0.2	3.18	0.50	0.02	0.00	3.70	50.84	”
4.1	260.067140	35.604648	0.98	0.86	3.85	2.5	2.51	14.06	3.11	1.66	21.35	71.33	Sys. 2 in Z15
4.2	260.074950	35.605803	0.84	0.86	3.85	5.3	23.21	0.05	14.48	5.31	43.04	71.33	”
4.3	260.065810	35.605677	0.86	0.86	3.85	2.7	1.10	0.04	3.77	2.03	6.94	71.33	”

Table 3. Table 3: MIFAL results for M1931

Id	R.A.	Dec	$z_{ima}$	$z_{sys}$	$z_{norm}$	$R_{find}$	${\chi }_{col}^2$	${\chi }_z^2$	${\chi }_{R_{find}}^2$	${\chi }_{SB}^2$	${\chi }_{ima}^2$	${\chi }_{sys}^2$	Comment
	${(}^{\circ })$	${(}^{\circ })$				${(}^{\prime \prime })$
1.1	292.960600	-26.569164	1.80	1.80	2.03	1.1	1.31	0.03	0.58	0.01	1.94	32.90	Sys. 1 in Z15
1.2	292.957990	-26.568611	1.41	1.80	2.03	0.8	2.44	2.35	0.32	2.99	8.10	32.90	”
1.3	292.949590	-26.570605	2.46	1.80	2.03	3.3	0.53	9.78	5.45	7.09	22.86	32.90	”
2.1	292.963980	-26.570404	2.22	2.23	2.03	3.9	20.38	0.00	7.70	0.87	28.95	72.29	False identification
2.2	292.966320	-26.571364	2.28	2.23	2.03	5.1	6.55	0.03	13.22	16.48	36.28	72.29	False identification
2.3	292.951220	-26.569077	2.23	2.23	2.03	0.2	6.69	0.00	0.02	0.34	7.05	72.29	Image 2.3 in Z15
3.1	292.955220	-26.583679	2.90	2.92	2.03	2.3	5.10	0.09	2.65	1.93	9.76	74.10	Sys. 4 in Z15
3.2	292.966700	-26.579142	3.06	2.92	2.03	3.7	2.88	0.86	6.85	0.17	10.75	74.10	Differs from Z15 but possible
3.3	292.967010	-26.578874	2.81	2.92	2.03	5.7	7.92	0.08	16.85	0.17	25.02	74.10	Nearby but likely unrelated
3.4	292.952110	-26.582897	2.95	2.92	2.03	2.7	22.12	0.22	3.72	2.51	28.57	74.10	Sys. 4 in Z15

Table 4. Table 4: Parameters used to construct the preliminary lens models

Cluster	$q$	$S$	$K$	$K_{gal}$	${\gamma }_{ex}$	${\varphi }_{ex}$	BCG1	BCG1	BCG1	BCG1	BCG1	BCG1	BCG2	BCG2	BCG2	BCG2	BCG2	BCG2
name		(pixels)				${(}^{\circ })$	RA ${(}^{\circ })$	DEC ${(}^{\circ })$	weight	$r_{core}\left[″\right]$	ellipticity	P.A. ${(}^{\circ })$	RA	DEC	weight	$r_{core}\left[″\right]$	ellipticity	P.A. ${(}^{\circ })$
M0329	1.3	150	1.5	0.07	0	0	52.423219	-2.1962308	5	0	0	0	52.413118	-2.1914246	3	0	0	0
M1720	1.3	150	1.5	0.07	0	0	260.06980	35.607312	5	0	0.123	87.2	—	—	—	—	—	—
M1931	1.3	150	1.2	0.07	0.1	82	292.95678	-26.57572	5	0	0.5a	82a	—	—	—	—	—	—

Table 5. Table 5: Thresholds used for running MIFAL

Cluster	Cluster	$\Delta r$	$\Delta z$	$\Delta SB$	min.	max.	${\chi }_{ima}^2$
name	redshift	(″)			images	images	cut
M0329	0.450	5.2	0.5	$300\%$	2	5	50
M1720	0.387	7.8	0.5	$300\%$	2	5	50
M1931	0.352	6.5	0.7	$50\%$	2	5	50

Equations

$$Q_{ij}=\frac{{\int }_A(x_i-{\bar{x}}_i)(x_j-{\bar{x}}_j)q(I(\vec{x}))d^{2}x}{{\int }_{A}q(I(\vec{x}))d^{2}x}.$$

$${\chi }_{ima}^2={\chi }_r^2+{\chi }_z^2+{\chi }_{colour}^2+{\chi }_{SB}^2,$$

$${\chi }_r^2=((x^{\prime }-x{)}^2+(y^{\prime }-y{)}^2)\text{leavevmode}/\text{leavevmode}{\sigma }_r^2,$$

$${\chi }_z^2=((z^{\prime }-\overline{z}{)}^2)\text{leavevmode}/\text{leavevmode}{\sigma }_z^2,$$

$${\chi }_{colour}^2=\sum _i(colou{r}_i^{\prime }-{\overline{colour}}_i{)}^2\text{leavevmode}/\text{leavevmode}{\sigma }_{colou{r}_i}^2,$$

$${\chi }_{SB}^2=((S{B}^{\prime }-SB{)}^2)\text{leavevmode}/\text{leavevmode}{\sigma }_{SB}^2,$$

$$\alpha (\theta )=\frac{4G\frac{2\pi {\Sigma }_{q0}}{2-q}{d}_l^{\text{leavevmode}1-q}}{c^2}\frac{d_{ls}}{d_s}{\theta }^{1-q},$$

$${\vec{\alpha }}_{gal}(\vec{\theta })\propto \sum _i{F}_i|\vec{\theta }-{\vec{\theta }}_i{|}^{1-q}\frac{\vec{\theta }-{\vec{\theta }}_i}{|\vec{\theta }-{\vec{\theta }}_i|}.$$

$${\alpha }_{gal,x}({\vec{\theta }}_m)=\sum _i{F}_i\frac{\Delta x_{mi}}{[(\Delta x_{mi}{)}^2+(\Delta y_{mi}{)}^2{]}^{q/2}},$$

$${\alpha }_{gal,y}({\vec{\theta }}_m)=\sum _i{F}_i\frac{\Delta y_{mi}}{[(\Delta x_{mi}{)}^2+(\Delta y_{mi}{)}^2{]}^{q/2}},$$

$${\alpha }_{DM,x}({\vec{\theta }}_m)=\sum _i{P}_i\frac{\Delta x_{mi}}{[(\Delta x_{mi}{)}^2+(\Delta y_{mi}{)}^2]},$$

$${\alpha }_{DM,y}({\vec{\theta }}_m)=\sum _i{P}_i\frac{\Delta y_{mi}}{[(\Delta x_{mi}{)}^2+(\Delta y_{mi}{)}^2]},$$

$${\alpha }_{ex,x}({\vec{\theta }}_m)=\gamma \cos (2\varphi )\Delta x_m+\gamma \sin (2\varphi )\Delta y_m$$

$${\alpha }_{ex,y}({\vec{\theta }}_m)=\gamma \sin (2\varphi )\Delta x_m-\gamma \cos (2\varphi )\Delta y_m,$$

$${\vec{\alpha }}_T(\vec{\theta })=K\text{leavevmode}[N_1{K}_{gal}{\vec{\alpha }}_{gal}(\vec{\theta })+N_2(1-K_{gal}){\vec{\alpha }}_{DM}(\vec{\theta })]+{\vec{\alpha }}_{ex}(\vec{\theta }),$$

Full text

††pagerange: MIFAL: fully automated Multiple-Image Finder ALgorithm for strong-lens modelling – proof of concept–B††pubyear: 2013

1 Introduction

The matter in the Universe is thought to consist of mostly an unseen component, known as Dark Matter (DM; e.g. Komatsu et al. 2011; Planck Collaboration et al. 2018). Although DM has still not been detected directly, evidence for it existence is manifested in many astrophysical phenomena and measurements on various scales (Einasto, 2009, for a review).

One of the most direct ways to map the unseen DM distribution in massive objects such as galaxy clusters is through the analysis of gravitational lensing signatures. Due to the high projected mass density of clusters, background sources are lensed, typically forming magnified and distorted images. Close to the cluster centre, in the strong-lensing (SL) regime where the mass density is higher than the critical value for lensing (see Bartelmann 2010; Kneib & Natarajan 2011), there are often seen multiple images of the same background source. As lensing is coupled to the deflecting mass distribution, the location and redshift of the multiple images supply constraints to map the DM distribution. Matching up multiple-image families, often acknowledged as the most complicated and time-consuming step in a SL analysis (see, e.g. Jullo et al., 2007; Bradač et al., 2008), is thus a crucial step for the modelling of the inner cluster (unseen) mass distribution.

Methods for analysing strong-lenses have been evolving over the past decades in response to higher quality data, mainly from space, and increasing computational power (e.g., Kneib et al. 1993; Keeton 2001; Broadhurst et al. 2005; Halkola et al. 2006; Liesenborgs et al. 2006; Jullo et al. 2007; Coe et al. 2008, see also Kneib & Natarajan 2011). Moreover, extensive lensing surveys recently undertaken with the Hubble Space Telescope (HST) such as the Cluster Lensing and Supernova Survey with Hubble (CLASH; Postman et al. 2012), the Reionization Lensing Cluster Survey (RELICS; Coe et al. 2019) and the ultra-deep Hubble Frontier Fields (Lotz et al., 2017), have allowed for high-end modelling using large numbers of multiple images (e.g. Jauzac et al., 2015; Caminha et al., 2017), and for the comparison of different lens modelling techniques (Treu et al., 2016), including with simulations (Meneghetti et al., 2017). Methods for SL analyses generally divide into two categories, namely parametric and so-called non-parametric or free-form methods. Free-form techniques (e.g., Diego et al. 2005; Liesenborgs et al. 2007; Coe et al. 2010; Williams & Saha 2011) usually make no assumptions regarding the underlying distribution of matter, and directly constrain the mass distribution from the location and redshift of multiple-images. Although this procedure is appealing due to its model-independency, the number of multiple-images readily identified in lensing clusters is often too small to obtain a high-resolution result (Ponente & Diego, 2011). For this reason, most free-form techniques often do not possess the capability to predict many multiple images in terms of location, exact shape, internal details and orientation (although few can also be successful in doing so, for example see Liesenborgs et al. 2008 or when combining the free-form method with the smaller-scale parametric galaxy contributions, e.g., Diego et al. 2015; Diego et al. 2016). Parametric techniques (e.g. Jullo et al., 2007; Johnson et al., 2014; Oguri et al., 2013; Grillo et al., 2015), on the other hand, make certain well-motivated assumptions regarding the underlying mass distribution and use physical parametrisations to represent it, and when well constructed, could be used to successfully find multiple-images and accurately predict the location and shape of other counter images (e.g. Treu et al., 2016; Meneghetti et al., 2017).

In fact, this prediction power is the key point for the work presented here. Based on previous studies in which we have analysed a few dozen galaxy clusters (e.g. Zitrin et al., 2009, 2015), we noticed that by constructing a simple preliminary yet well-guessed lens model, multiple-images can be readily found using the model, by sending manually each arc-like image to the source plane and back through the lens to the image plane, and then looking for similar looking (and similarly redshifted) objects where the model predicts them. Thus, our goal here is to simply automatise this procedure, so that clusters could be analysed “blindly”, in an automated manner.

Despite the various methods available for SL analyses or mass modelling, only several dozen clusters have been analysed to date with regards to their SL features (albeit with rapidly increasing numbers), with the vast majority only in the last decade or so (e.g. Smith et al., 2005; Richard et al., 2010; Zitrin et al., 2015; Kawamata et al., 2016; Cerny et al., 2018; Acebron et al., 2018). One of the main reasons, apart from the finite amount of high-resolution data particularly from space, is that nearly all methods are often time-consuming and the analysis has to be performed on a particular basis with significant intervention and fine-tuning by the user. Ongoing large sky-surveys and the amount of data taken by high-resolution space missions such as the HST, the upcoming launch of the James Webb Space Telescope (JWST), and especially the Euclid and WFIRST missions, call for a more (statistically) efficient and automated method which could be operated on a broad database.

Broad studies of galaxy clusters are of high importance, as they will be able to probe and put independent constraints on the acceptable $\Lambda$CDM paradigm (e.g. Wambsganss et al., 1995; Dalal et al., 2004; Oguri & Blandford, 2009), which is often confronted by various studies (e.g. Hennawi et al., 2007; Broadhurst et al., 2008; Broadhurst & Barkana, 2008; Sadeh & Rephaeli, 2008; Puchwein & Hilbert, 2009; Sereno et al., 2010; Gralla et al., 2011; Meneghetti et al., 2010; Umetsu et al., 2011, see also Kroupa 2012; although most show a mild discrepancy with debatable consequences), and on related and important observables such as the geometry of the Universe (e.g. Jullo et al., 2010), the mass-concentration relation of clusters (e.g. Oguri et al., 2012; Merten et al., 2015; Umetsu et al., 2014, see also von der Linden et al. 2014; Hoekstra et al. 2015; Okabe & Smith 2016), the density profiles in the inner core (Navarro et al., 1996; Newman et al., 2013), the high-end of the mass function and its evolvement in redshift, or the (in)homogeneity of the Universe. This is especially relevant given the recent tension between the abundance of massive clusters found by the Planck satellite, and the Planck CMB cosmology (Planck Collaboration et al., 2015). As signatures of lensing are coupled, via the mass distribution, to the gravitational potential, it is also one of the very few phenomena that might, conceivably, hold the key for independently testing the nature of DM, or its alternatives such as MOND or TeVes/MOG (e.g., Milgrom 1983; Bekenstein 2012, see also Amendola et al. 2008; Tian et al. 2013), Early Dark Energy propositions (see Fedeli & Bartelmann, 2007) and other innovative ideas such as the large local void in a Lemaître-Tolman-Bondi (LTB) universe (see Redlich et al., 2014), or non-Gaussianity of early perturbations (Chongchitnan & Silk, 2012), and free-streaming neutrinos (Emami et al., 2017).

In addition, since galaxy clusters magnify background sources due to the lensing effect, lens models and their magnification maps are central in the discovery and higher-resolution studies of early, high-$z$ galaxies (e.g. Bradley et al., 2012; Bouwens et al., 2009; Zheng et al., 2009; Zackrisson et al., 2012; Bradač et al., 2012; Ishigaki et al., 2015; Hashimoto et al., 2018). In particular, recent deep observations along with successful mass modelling, enable also the detection of multiply-imaged, higher-redshift galaxies (e.g. Franx et al., 1997; Kneib et al., 2004; Richard et al., 2011; Coe et al., 2013).

Here, we combine three successful (and independent) tools, with the goal of “paving the road” towards a completely automated and innovative (albeit simple) multiple-image finder, or equivalently, SL analyser. We begin with a sophisticated arc-finder (Seidel & Bartelmann, 2007) to detect stretched and arc-shaped objects across the cluster central field (that are likely lensed and, potentially, multiply imaged). The photometry of these arcs could then either be input into a photometric-redshift code (such as the Bayesian Photometric Redshift, BPZ; Benítez 2000; Benítez et al. 2004; Coe et al. 2006) in order to assess the corresponding source redshift, or in our case, cross-matched with objects in the public CLASH catalogues, which include multi-band photometry and photometric redshift estimates for each object. We construct an initial mass model using a revised version of the (semi-)parametric light-traces-mass (LTM) SL modelling method of Zitrin et al. (2009), which has a minimum of free parameters and is initially well constrained even without using multiple-images as input, and is therefore optimal for finding and confirming new multiple-images (e.g. Zitrin et al., 2015). We make use of the initial mass model to automatically lens each arc candidate back to the source plane via the lens equation, with the measured photometric redshift, and relens it through the lens to predict the appearance of other counter images of the same background source. In each of the predicted locations in the image-plane, the arc catalogue is searched for close-by objects with similar photometric redshifts, surface brightness, and colours. The new multiple-image families thus found can then, possibly iteratively, be used in order to refine the mass model. Although it constitutes another important step towards automated lensing analyses potentially applicable to high-end data or large sky surveys (e.g. Zitrin et al., 2012a; Wong et al., 2012; Stapelberg et al., 2019), the outlined procedure is clearly simple in essence. Also, each of the three independent ingredients of the algorithm could be in principle easily replaced with other alternatives, or future improvements.

We apply here this method to the X-ray luminous (Ebeling et al., 2010) clusters MACS J0329.6-0211, MACS J1720.2+3536, and MACS J1931.8-2635 (hereafter M0329, M1720, and M1931, respectively), in observations taken with the HST as part of the CLASH treasury program (Postman et al., 2012). The clusters were imaged to a depth of $\sim 27$ AB in 16 filters with the Wide Field Camera 3 (WFC3) UVIS and IR cameras, and the Advanced Camera for Surveys (ACS) WFC. The images were reduced and mosaiced ($0.065\arcsec/pix$) using standard techniques implemented in the MosaicDrizzle pipeline (Koekemoer et al., 2002, 2011). CLASH also delivered HST photometry for all identified objects performed automatically by SExtractor (Bertin & Arnouts, 1996), and accompanying photometric redshifts estimated for all the galaxies in the field using the full 16-band UVIS/ACS/WFC3-IR photometry via the BPZ (Benítez, 2000; Benítez et al., 2004; Coe et al., 2006) program (Postman et al. 2012; Molino et al. 2017 for additional details).

The paper is organised as follows: In §2 we detail the multiple-image finder and SL analysis algorithms. In §3 we report the results of applying this procedure to M0329, M1720, and M1931, and compare the results with previous analyses. The results are discussed in §4 and the work is summarised in §5. The code to construct an initial blind model is given in Appendix A, and the relevant parameters that were used are given in Appendix B. Throughout this paper we adopt a concordance $\Lambda$CDM cosmology with ($\Omega_{\rm m0}=0.3$, $\Omega_{\Lambda 0}=0.7$, $H_{0}=100$ $h$ km s*-1Mpc-1*, with $h=0.7$).

2 The Multiple-Image Finder ALgorithm

We present here the Multiple-Image Finder ALgorithm (MIFAL)111Based on work presented in M. Carrasco’s PhD Dissertation. The general idea is as follows: given (i) a list of arclet candidates near the cluster core, (ii) their photometry and photometric redshifts, and (iii) a preliminary lens model, run over all arclets, send each one to the source plane and back to obtain predicted locations for counter images, and search around these positions for objects with similar photometric redshift, Spectral Energy Distribution (SED) or colours, and surface brightness. Rank the potential candidates for each lensed arc by similarity, and select the best combinations. We now detail the different ingredients of the algorithm.

2.1 The arc candidate catalogue

The first step we apply here for automatically searching for multiple images is to primarily identify as many lensed galaxy candidates as possible in the cluster lens field. Since lensed features are typically sheared, the first step is to build an arc catalogue, which constitutes a list of likely lensed candidates. This is the first ingredient needed to run MIFAL.

The arc catalogue is constructed here using the Arcfinder algorithm (Seidel & Bartelmann, 2007), which employs the first and second brightness moments for the detection of elongated objects across the frame. As detailed in Seidel & Bartelmann (2007), the algorithm is sufficiently robust to detect such features even if their surface brightness is near the pixel noise of the image, yet the amount of spurious detections remains (relatively) low.

The algorithm subdivides the image into a grid of overlapping, circular cells (see Fig. 1a) which are iteratively shifted towards their local centre of brightness in their immediate neighborhood (see Fig. 1b to 1e). The centre of brightness is defined as the first moment $\int_{A}\vec{x}d^{2}x=\bar{\vec{x}}$, where $A$ is the cell area and $q(I)$ is a weighting function that depends on the image intensity $I(\vec{x})$. The code then computes the cell ellipticity $(Q_{11}-Q_{22}+2iQ_{12})/(Q_{11}+Q_{22})$ using their second brightness moments

[TABLE]

The ellipticity of each cell provides a natural measure of its orientation in the image and allows us to compute their angular separation. Cells oriented in the same direction, and spatially aligned, are then combined into initial objects using a simple coherence measure (see Fig. 1f); essentially the product of the cosine of the angle between cell orientations and a (scaled) perpendicular cell separation, with a friends-of-friends type algorithm. At this point, the objects are nothing more than sets of correlated cells (see 1g) where one cannot directly infer much morphological information, and there are still many spurious detections, e.g. for spatially connected point sources. A large number of spurious detections is removed with a filter interpreting the brightness distribution in each cell to remove those which are unlikely to belong to an elongated object. For objects still containing a sufficient number of valid cells, several consecutive steps compute isophote contours which allow the computation of basic properties like length, length-to-width and signal-to-noise. These are then used to choose the most likely arc-shaped objects, and further reduce the number of spurious detections. Some further filtering and noise cleaning procedures take place, for example to remove elongated objects which are clearly not gravitational arcs (such as diffraction spikes, spiral galaxy arms, edge-on galaxies, etc.), resulting in a final catalogue of arc candidates (see 1i).

To obtain an arc catalogue as complete as possible, we run the Arcfinder on a detection image constructed from all optical ACS and WFC3-IR data. It should also be noted that the arc detection part takes only several seconds on a standard CPU.

Two main caveats follow the use of an automated arcfinder compared to examining each arc candidate by eye. First, some spurious detections of elongated, un-lensed galaxies may appear in the catalogue. These, however, will be mostly filtered out at a later stage based on their photometric redshifts and on the mass model. Second, images with a feeble elongation or distortion may not be detected even if they are multiply lensed (see example the multiply-imaged spiral galaxy in MACS J1149.5+2223, Zitrin & Broadhurst 2009). Our goal throughout, though, is not to detect all previously known multiple-images, but rather, enough of them so that a proper mass model could be automatically constructed.

2.2 Photometry and photometric redshifts

The second step in the procedure is to obtain, for each arc candidate, fluxes in each band and a photometric redshift estimate. In principle, this can be done by measuring the fluxes using the Arcfinder contours directly, although here we take on an alternative route, matching each arc in the arc-candidate catalogue to an object from the public CLASH catalogues by Molino et al. (2017), which include aperture-matched PSF-corrected multi-band photometry, and enhanced Bayesian Photometric Redshifts (BPZ2.0; Benítez 2000; Benítez et al. 2004; Coe et al. 2006; Molino et al. 2017), based on the multi-drizzled (0.065 $\arcsec/pix$) science images produced by CLASH (Koekemoer et al., 2011). These redshifts will later be used to compute the (relative) lensing distances to re-lens the candidate arcs through the mass model and potentially, further constrain the model in the cases where multiple image systems were found. After matching arcs to their photometry and photometric redshifts, the arc catalogue is then cleaned from objects whose photometric redshift is lower than $z_{l}$ + 0.3, where $z_{l}$ is the redshift of the cluster, and from objects lying outside the field-of-view (FOV) of the initial lens model.

2.3 The preliminary lens model

The third and last ingredient needed to run MIFAL is a lens model. The preliminary mass model we use here is based on the semi-parametric method of (Zitrin et al., 2009, 2015, see also Broadhurst et al. 2005), which is based on the assumption that light traces mass (hereafter LTM). The mass model can be quickly constructed by adopting a simple representation of the cluster member galaxies and the underlying DM component, together with a prior guess of the ratio between the two components, and a certain (effective) mass-to-light ratio ($M/L$) to account for the overall normalisation. Details of this mass modelling method and it prediction power can be found in, e.g., Zitrin et al. (2009, 2012a, 2015). Explicit details regarding the code and its implementation are given in Appendix A, and here we give only a brief overview.

The mass model consists of three components. The first component are the cluster galaxies, whose surface mass density is modelled as a power-law of slope $q$, scaled by their luminosity. The superposition of the mass density contributions from all red-sequence galaxies (brighter than a certain threshold), represents the lumpy, galaxy component. All galaxies are modelled as circular, aside from one or two Brightest Cluster Galaxies, for which ellipticity, and potentially a core, can be assigned. This mass distribution is then smoothed using a 2D Gaussian kernel in Fourier space, to provide a model for the DM distribution in the cluster halo, which is the second component of the model. The deflection fields from these two components are then normalised onto a similar arbitrary scale, and added together with a relative weight ($K_{gal}$). The resulting deflection field is calculated over the input image grid (or a lower resolution version of it) and scaled as desired (with parameter $K$). Finally, the third component is an external shear which can be added to the deflection field for further flexibility. In all, there are six free parameters: the power law of the galaxies’ density profile, the degree of smoothing (Gaussian width), the ratio between the galaxy and DM contributions, the overall normalisation (i.e., the mass scaling of both the smooth and galaxy components), and the strength and direction of the external shear. This simple procedure seems to produce relatively well-guessed initial mass models that can be used to identify multiple images, and in turn, be iteratively refined. This makes the procedure ideal for our purpose.

Our goal here is to first construct a reliable initial mass model to allow for the automated identification of multiple images. For that purpose, as we now detail, we fix all free parameters to either averaged values, or to values indicated by the input galaxy catalogue. As we have mentioned in previous work (e.g. Zitrin et al., 2009), the choice of the power-law ($q$) and smoothing Gaussian width ($S$) only mildly affects the reproduction of multiple images. These parameters mainly affect the mass profile, so that the assumption is that nearly all systems can, in principle, be reproduced with any reasonable choice of $q$ and $S$. Therefore we set these two parameters to typical values deduced from our previous analyses (see also Zitrin et al. 2012a), adopting here a power-law $q$=1.3, and a smoothing Gaussian width of $S$=150 pixels (0.065 ″/pix). Similarly, for the relative galaxy to DM ratio, we use a typical value of $K_{gal}=0.07$. The value of the overall normalisation, $K$, basically determines to which source redshift the model is scaled, and is unknown. Instead we can adopt a typical value, and iterate on different cosmological distances ($D_{ls}/D_{s}$ normalisations) searching for multiple images given each normalisation, and then use the multiple images found to normalise the model. Moreover, the overall normalisation $K$ correlates empirically with the size of the lens (Taiber et al., in preparation). We use a typical value of $K=1.5$ , which – based on previous analyses – corresponds roughly to an Einstein radius of $\sim 20\arcsec$. A moderate external shear ($\gamma_{ex}=0.1$) is used along the BCG direction, in cases where the BCG is highly elliptical; otherwise no external shear is used in the preliminary model. If external shear as such is used, which adds to the deflection field, then a smaller $K=1.2$ is adopted (after rescaling of the galaxy and DM deflection fields such that the mean of their absolute value across the field of view is 200 pixels, equal to 13″; see Appendix A for more details).

As the input red-sequence list, including the BCG ellipticity parameters, we use the same member-galaxy catalogues used by Zitrin et al. (2015), based mainly on measurements with SExtractor (Bertin & Arnouts, 1996) and a simple colour-magnitude red-sequence cut. In addition, we note that the flux measurements of BCGs, which are often much more diffuse than other cluster members, do not seem to trace well the encompassed mass of the BCGs (or, perhaps, BCGs may lie effectively on a different mass-to-light ratio compared to other galaxies, at least in our modelling scheme). We have found empirically in our previous modelling using this method, that the BCG $M/L$ ratio needs to be multiplied by a factor $\sim 3-5$, which we adopt here as well. We thus have now a single, preliminary mass model for each cluster (see Fig. 3) that can be used to find multiple images. It is worth mentioning that the construction of this initial mass model takes only a few seconds on a typical computer. Given it is key to the success of the automated multiple image finding, we make our code for constructing the initial lens-models publicly available222\colorbluehttps://github.com/adizitrin/initialLTM, and detail its implementation in Appendix A below. Table 4 lists the parameter values used to construct each of the three models used here.

2.4 Finding multiple images

We describe now the simple algorithm used by the code for finding multiple images. Since the redshift, or $D_{ls}/D_{s}$ normalisation to which the preliminary mass model is scaled is not known a priori, to explore the relevant range the code probes 25 discrete values of $D_{ls}/D_{s}(z_{s})$ normalisations, from source redshifts $z_{s}=z_{l}+0.3$ to the maximum source redshift given in the arc catalogue, where $z_{l}$ is fixed to the cluster redshift and $z_{s}$ increases such that the increments in $D_{ls}/D_{s}(z_{s})$ are linear. In practice the 25 realisations are grouped in five principal $D_{ls}/D_{s}$ bins, each bin contains five sub-bins, corresponding to five different source redshifts, so that systems from the different sub-bins are merged together. The process now described refers to one $D_{ls}/D_{s}$ normalisation, and hence is repeated five times, once for each principal normalisation.

2.4.1 Preliminary multiple image file

The preliminary lens model is used to automatically project each candidate in the arc-candidate catalogue to the source plane and then back, to predict the location of possible counter images. This is done using the best-fit photometric redshift of the re-lensed arc as input. The program then runs over all the candidates in the catalogue searching for objects near the predicted locations and with similar photometric redshift and surface brightness as the re-lensed candidate. These will later be ranked according to their similarity to the re-lensed arc (§2.4.3). The process is repeated for the five sub-bins of the $D_{ls}/D_{s}$ normalisation adopted. To reduce spurious detections, thresholds are defined past which potential counter image candidates are not considered. The typical thresholds we probe are (in different iterations and for the different clusters – see Table 5 for the final parameters for each cluster) $\Delta z=0.2$ to $\Delta z=0.7$ for the redshift difference, $\Delta r=5\arcsec-10\arcsec$ for the search radius, and a relative $\Delta SB=20\%$ to $\Delta SB=300\%$, compared to the input re-lensed arc.

Once MIFAL has listed all the possible multiple images for the re-lensed candidate, for all candidates in the file and for the five normalisation sub-bins, it takes on a few organising steps – such as removal of single images, combination of repeated images (i.e., if these correspond to the same arc), and coaddition of the lists from the five different sub-bins – towards creating a preliminary list of potential multiple image systems for the adopted $D_{ls}/D_{s}$ normalisation. It then calculates the median redshift of each temporary system, $z_{sys}$, using the redshift of all members of that system, and registers the number of images in each system. At the end of this stage, each candidate in our catalogue has its own candidate multiple-image system and its respective $z_{sys}$. The file now is a preliminary list of the potential multiple image systems. Note that at this point this file may include some repetitions of systems, in the sense that if arc $i$ was relensed, and arcs $j$ and $k$ are listed as possible counter images so that these make a multiple image system, then the same system may appear again when listing candidate counter images for arcs $j$, and $k$. Also possible is the repetition of systems that are reproduced in more than one normalisation sub-bin.

2.4.2 The semifinal list

The program now makes a few more reorganising steps to the preliminary multiple image list. First, it assigns a temporary index to each system. The index assignment is made following a simple count of the different systems, with the exclusion that systems originating from the same relensed arc, and that have similar median redshifts ($<\Delta z$), are identified as the same system and given the same system number. This allows for grouping the same systems together (from different sub-bins, for example), which is the second organising step. Then, the code merges images that repeat in the same system, averaging their properties where differences exist (for example, the $D_{ls}/D_{s}$ may differ if images originate in different sub-bins), and recalculates the median redshift of each system.

The next stage is to obtain the so-called semifinal list of multiple images (per adopted $D_{ls}/D_{s}$ normalisation), merging together systems that have at least 2 images in common and similar redshift ($z_{sys}\pm\Delta z$), merging again images that repeat in each new system, and recalculating the system’s properties such as median redshift and number of images.

2.4.3 Final multiple image lists

To obtain a final list of multiple images, MIFAL now assigns to each multiple-image candidate a $\chi^{2}$ grade to estimate its likelihood of being a member of a certain family. The $\chi^{2}$ for each image is defined by:

[TABLE]

where the first component is the $\chi^{2}$ of the image location with respect to the expected position. This is defined by:

[TABLE]

with [x*′,y′*] and [x,y] being the position of the candidate, and the model-predicted multiple-image, respectively, and $\sigma_{r}$ the expected image-plane reproduction uncertainty. Motivated by studies that probed the contribution from large-scale structure and matter along the line of sight to the lensing signal (e.g. Host, 2012), we set $\sigma_{r}=1.4^{\prime\prime}$ throughout. The second component is defined as:

[TABLE]

where $z^{\prime}$, $\overline{z}$, and $\sigma^{2}_{z}$, represent the measured photometric redshift, the median system redshift, and the dispersion of photometric redshift measurements within the system, respectively. The third component is defined as:

[TABLE]

where, $colour^{\prime}_{i}$, $\overline{colour}_{i}$, and $\sigma^{2}_{colour_{i}}$ are the $i^{\prime}th$ colour of the candidate in question, the median value of this colour in the candidate system (from the other images found for that system excluding the image in question), and its dispersion, respectively. While any combination of colours can be chosen, and while CLASH has more filters available, in practice we limit our case study to using only the following colours (over which the sum is performed): F475W – F435W, F625W – F475W, F814W – F625W, F814W – F475W, F105W – F814W, F125W – F105W, and F140W – F125W. If there are non detections in some bands, or if the candidate arc is a dropout, i.e., a Lyman break galaxy, for example, then the relevant colours are ignored. The fourth component is defined as:

[TABLE]

where $SB^{\prime}$ is a surface brightness measure of the candidate that was lensed (in practice we use the F814W flux with a fixed, arbitrary zero point, divided by the area used for the flux measurement of each candidate). $SB$ is the surface brightness measure of the candidate counter image, and $\sigma^{2}_{SB}$ the respective $1-\sigma$ error given by the dispersion of surface brightness measures within the system.

After MIFAL assigns a $\chi^{2}$ grade to each candidate, it selects the so called final, most likely images per system that are below a user-defined $\chi^{2}_{ima}$ threshold (which we take throughout as 50). In practice, the allowed number of images in a system is typically being limited to at least 2 (to minimise spurious detections, a minimum of 3 images can be used). The maximum number of images in a system we use is 4. Thus, if the semifinal system has between two and four images below the threshold, the program now adopts it as is. If it has five images or more (some of the semifinal systems may have them), the algorithm selects the four images with the lowest $\chi^{2}_{ima}$ that meet the predefined $\chi^{2}_{ima}$ threshold. Then, it calculates the total $\chi^{2}_{sys}$ of the system, simply adding the $\chi^{2}_{ima}$ of the selected images. For now, we allow also for multiple images to be considered in more than one system (with different $\chi^{2}_{ima}$ weights), although this can clearly be modified to choose only a single system that corresponds to the best $\chi^{2}$. After each system has been finalised, the program recalculates again its properties such as median redshift and updated number of images.

2.4.4 Repetition and normalisation

Once the program has run through all semifinal systems, filtering and updating them using the said $\chi^{2}$, we obtain the final multiple-image systems catalogue. Since, as mentioned, the code repeats the process described above for 5 principal $D_{ls}/D_{s}(z_{s})$ normalisation bins, we end up with one multiple image catalogue for each $D_{ls}/D_{s}$ bin. This is the final stage of the current process. However, since multiple images have now been found these can be used to refine the model (e.g., using a quick Monte-Carlo Markov Chain (MCMC) minimisation), hence calibrating it to a specific chosen redshift. This calibrated model could now, in principle (i.e., we do not perform this additional step here) be used to repeat the MIFAL process, potentially refining the list of multiple images.

2.4.5 Adding lens configuration and parity information

The simple $\chi^{2}$ ranking uses only the metric defined above for grading each multiple image and system: the similarity in terms of colours, redshift, position, and surface brightness. In addition, typical lensing configurations, as well as the parity of the images, can be used to improve further the selection of multiple image families. For completeness, we briefly mention this possibility (partly incorporated in the code at present), although we do not present results based on these features in the current work, which is aimed merely to constitute a first proof-of-concept for the automated finding of multiple images.

If one assumes that most clusters would show, roughly, lensing configurations as those given by smooth elliptical mass distributions, extra information can be incorporated to improve the multiple-image list: for example, to exclude two nearby images that are on the same side of the critical curve (and so the method could choose the preferred one based on its $\chi^{2}$ ranking). Some of the most typical lens configurations are the Einstein cross, cusp arc, and fold arc systems, which have five multiple-images each (four distributed around the lens and one, often smaller and fainter than the rest, close to the cluster centre). As is common in massive elongated clusters, cusp like systems often show only three highly magnified close images; a configuration that could be considered independently, as well. Hence the method can use these four lens configurations to add position constraints in the selection process of the final multiple images. In practice, this means the program would take each semifinal system and check if there are image combinations that match with some of these four lens configurations, including the parity, and save or prioritise the successful image combinations. The program would assign a $\chi^{2}$ to each successful combination, that is combined with the previous $\chi^{2}$ ranking.

3 Results

As the first proof-of-concept for the presented method and the feasibility for detecting multiple images in a “blind” fashion, we analyse three CLASH clusters: M0329, M1720, and MACS1931. These clusters were previously analysed in the framework of the CLASH program by Zitrin et al. (2015), who found several multiple imaged systems in each of them (see also Caminha et al. 2019). The first two clusters, M0329 and M1720, show a complex, unrelaxed, or substructured central mass distribution, whereas the last one, M1931, shows quite an elongated, elliptical-like shape.

As detailed in §2, we run the Arcfinder on all three cluster detection images, and then cross-match the arcfinder’s results to the CLASH catalogues to obtain the photometry and photometric redshift of each object. We then clean the catalogue from objects lying outside the model’s FOV, or that do not lie far enough in redshift beyond the cluster’s redshift. We construct preliminary mass models for each cluster, as detailed in §2.3. The preliminary, “blind” mass models of these clusters, constructed without any input regarding the multiple images, are shown in terms of their magnification maps in Fig. 3, and the parameters used for constructing these are given in Appendix B. We run MIFAL on these three clusters, and detail here the results. To probe the method with respect to the input thresholds, we start with the same thresholds for all three clusters [$\Delta z=0.5$, $\Delta r=6.5\arcsec$, relative $\Delta SB=50\%$] and slightly refine them per cluster, to probe the effect on the results. The most likely multiple-image systems found by the automated procedure are reported in Tables 1, 2, and 3.

MACS J0329.6-0211: Zitrin et al. (see also 2012b, 2015) reported six sets of likely lensed systems in M0329, although only three systems they considered secure and were used to constrain their model. For comparison, Fig. 4 shows the multiple-image systems found by the automated MIFAL for a normalisation redshift of $z_{norm}=1.46$, and the following thresholds [$\Delta z=0.5$, $\Delta r=5.2\arcsec$, relative $\Delta SB=300\%$]. Without any multiple image constraints as input, MIFAL fully recovers the three multiple-image systems, including the quadruply imaged $z\sim 6$ dropout galaxy. Another candidate system proposed by Zitrin et al. (2015) is also fully recovered. In addition, there is a seemingly false two-image system also proposed (system 4 here), which can be discarded by either accepting only systems that have three and above multiple images, or by applying thresholds to the different $\chi^{2}$’s. Note also that with the default thresholds MIFAL initially suggested two other (most likely unreal) systems, that were either excluded by the code based on their very high $\chi^{2}$ – that is 1-2 orders of magnitude higher than the secure systems – or were not detected when revising the thresholds as above. Especially, lowering $\Delta r$ to $5.2\arcsec$ excluded some other suggested counter image candidates. We also note that when $\Delta SB$ was set initially to $50\%$, the $z\sim 6$ dropout system went missing (because its F814W flux is noisy) - but was recovered within the relative $\Delta SB=50\%$ threshold if we instead calculated the SB based on a redder band (F125W). In that case all MIFAL systems seen in Fig. 4 are discovered, aside from image 2.1, and with no spurious detections. We thus conclude that of the probed combinations, that is the most successful set of thresholds for M0329. The systems are also listed in Table 1. These are sufficient to fully constrain a proper lens model for the cluster.

MACS J1720+3536: Zitrin et al. (2015) reported 5 secure systems in M1720, and one additional candidate system (system 5 in their notation). Fig. 5 shows the most likely multiple-image systems found by MIFAL for a normalisation redshift of $z_{norm}=3.85$, and the following thresholds [$\Delta z=0.5$, $\Delta r=7.8\arcsec$, relative $\Delta SB=300\%$]. Our automated algorithm recovers three of the five secure systems by Zitrin et al. (2015), including a distinct, bright two-image red system. The two remaining systems from Zitrin et al. (2015) we do not expect to be recovered by MIFAL, since many of their images are not present in the automatic arc detection+photometric redshift catalogue. Among the three systems that are recovered by MIFAL, two (systems 3 & 4 here) are exactly reproduced, while for system 2 three images seem to be well matched but an additional image is also included (image 2.4), which may not be related. Additionally, Zitrin et al. (2015) considered another possible counter image as image 2.3, but MIFAL’s choice seems reasonable. The proposed two-image system number 1, seems to be falsely detected altogether. If we use the default thresholds [$\Delta z=0.5$, $\Delta r=6.5\arcsec$, relative $\Delta SB=50\%$], then only system 2 and false system 1 are recovered. Both the SB threshold and the radius have to be increased to detect the remaining systems. The multiple-image systems and their properties, for the chosen thresholds, can be found in Table 2. Therefore, in this case our algorithm was able to recover $\sim 2.5$ out of 5 previously detected systems found by aid of a visual inspection. Still, the procedure supplies enough secure constraints to revise the initial lens model.

MACS J1931-2635: Zitrin et al. (2015) reported 4 secure systems in M1931 (and three other potential candidate systems, that we disregard). Fig. 6 shows the most likely multiple-image systems found by MIFAL for a normalisation redshift of $z_{norm}=2.03$, and the following thresholds [$\Delta z=0.7$, $\Delta r=6.5\arcsec$, relative $\Delta SB=50\%$]. It is first noted that MIFAL recovers well system 1. System 1 shows a large scatter in photometric redshifts, probably due to contaminated photometry by nearby stars or cluster members, and so the initial default $\Delta z$ had to be raised to $0.7$ to recover it. MIFAL also recovers a second system, namely two images of system 3 here (3.1 & 3.4; images 4.1 & 4.2 in the notation of Zitrin et al. 2015). As the counter image for this system expected on the other side of the cluster MIFAL assigns a different image than the one Zitrin et al. (2015) found, but appears to be similar-looking as well. In fact MIFAL lists two nearby images as the counter image but if these are not a single object, only one of them can be the true counter image. Next, we do not expect MIFAL to reproduce system 2 from Zitrin et al. (2015) since two of its images do not appear in the arc+photometry catalogue, and indeed MIFAL assigns two other (false) counter images to image 2.3 (image 2.3 also in their notation), which is the only one from this system that appears in the catalogue. Similarly, we do not expect MIFAL to reproduce system 3 from Zitrin et al. (2015), because it is only partly listed in the arc+photometry catalogue and hence also absent from the final multiple-image catalogue. The multiple image systems are listed in Table 3. Overall MIFAL does well in reproducing the systems that appear in the catalogue for this cluster, with some differences in the image designated as 3.3.

4 Discussion

We have probed the feasibility of finding multiple images automatically, on three massive clusters imaged with the HST. The outlined procedure was operated on the three clusters and compared to previous analyses. The method seems to be quite successful: without any multiple images used as input, it recovers all secure known multiple images in M0329, and several secure multiple-image systems in M1720, and M1931. Given the trade off between completeness of true images and spurious detections, some erroneous multiple images were identified as well, and thresholds had to be refined to maximise the true-to-false identification ratio. The number of identified secure systems is sufficient to fully constrain, or refine, the initial lens models. In fact, the output multiple image list can be used as input for any lens modelling technique and so the final solution does not necessarily have to be limited to the LTM technique.

The method we presented for automatically finding multiple images consists of three components: the arc candidate catalogue, a photometry and photometric redshift catalogue, and the initial lens model. Each of these constitute a crucial ingredient for the success of the method. The arc-finder and the photometry are known to operate very well, yielding stable results. The promise of the process is perhaps driven by the initial lens model – such blind automated lens modeling for SL clusters are not yet too common – constructed following a revised version of the LTM methodology (Zitrin et al., 2009, 2015, see details of the LTM code in Appendix A).

The LTM is very promising for the automated analysis of SL clusters, given that the initially guessed model, based only on the distribution of cluster members and their luminosity, is already reasonable enough so that multiple images can be found using that model (see as a few examples Zitrin et al. 2009, 2015, and references therein). However, the LTM initial-guess model, even while evidently successful, cannot be expected to trace the mass distribution exactly, and hence its success could differ from cluster to cluster. In addition, compared to analytic models such as elliptical or pseudo elliptical mass distributions, the LTM methodology has a few limitations. First, the mass distribution is often rounder, where the required ellipticity of the critical curves often needed, is introduced as an external shear contribution. This results in a limited overall elongation of the model. Second, the mass profile is often shallower, which results in higher magnifications compared to common analytical models. Finally, it is susceptible to galaxies that lie outside the adopted M/L ratio, which can significantly bias the resulting model.

This range of success levels is only partly manifested in the analysis we perform here, where the initial model for the three clusters seems to be successful in finding all or most multiple images using the automated procedure, at least, those that are included in the arc-photometry catalogue. MIFAL has a harder time reproducing some images but this, however, seems to be in large part a more direct result of noisy or contaminated photometry, and erroneous photometric redshifts, in some cases.

The current work aims to constitute a proof-of-concept for the automated finding of multiple images. Clearly, there are many avenues for improving and further scrutinising the presented procedure. Such improvements, for example, include the use of additional information such as known lensing configurations, parity and/or shear or magnification of the images, using the full photometric-redshift probability distribution of each arc (rather then just the best-fit value), or explicit use of the structure of the images (i.e., internal details) through a full-arc reconstruction comparison for each image. In addition, ways to improve the initial model guess, including the usage of other alternative parametrisations, should be probed as well. Similarly, one major improvement could be to find the effective scaling relation between the luminosity and mass that corresponds to the current methodology, as was done, for example, by (Zitrin et al., 2012a). It might also be beneficial to examine whether the weighting of the different $\chi^{2}$’s should remain as presented or perhaps modified, to give more weight to certain, more stable measurements.

While the success of the method is promising, there is now a need to continue operating the method on more clusters, to continue examining its performance in different configurations.

Finally, it should also be mentioned that the current procedure is one example of many possible procedures, including, conceivably, those that combine machine learning algorithms for image identification.

5 Summary

In this work we presented an innovative Multiple-Image Finder ALgorithm (MIFAL) designed to automatically find multiple images in galaxy cluster lenses, so that their mass model could be efficiently and automatically constrained. We combine an arcfinder algorithm with CLASH photometric redshift measurements, along with a preliminary mass model, to physically match together multiple-image systems in an automated (“blind”) manner. We obtain a robust assessment of the likelihood of each arc to belong to one of the multiple-image systems, as well as the preferred redshift for the different systems. MIFAL then selects the most likely multiple images for each system based on the assigned grades, to finally construct a catalogue of multiple-image systems.

We applied MIFAL to the recently studied galaxy clusters M0329, M1720 and M1931 in deep CLASH/HST images. We compared the results of our automated procedure with the results by the conventional SL analysis where multiple images were verified by eye and using parametric lens reconstruction methods. Our automated SL analysis by MIFAL blindly recovered all three known systems in M0329, and about two or three, or half of the systems in M1720 and M1931. In cases where multiple images and/or systems were not picked by the automated procedure, this was mainly due the arcs being faint and missing from the original arc catalogue. Some false identification also occurred, mostly due to ambiguous or erroneous photometric redshift measurements, or in part, inaccurate initial lens model.

Although more comparisons are required for a more robust assessment of the proposed algorithm accuracy, and although the basic proof-of-concept procedure we outline here is perhaps not yet fully competitive with regular SL analyses, the potential for expanding this to completely automated high-end SL analyses is evident. We hope the present work will constitute another step towards fully automatising SL analyses, as a standard tool for studying cluster mass distributions of very large samples, with the HST and successors such as JWST, and particularly, Euclid and WFIRST.

Finally we wish to highlight that the multiple-image automated finding process takes no more than a few minutes on a normal CPU. A second iteration that includes a refinement of the first lens model with the newly found multiple images, can take up to several hours due to the multi-dimensional minimisation.

acknowledgements

The authors are grateful for useful discussions with Matthias Bartelmann. This work was supported in part by the transregional collaborative research centre TR 33 “The Dark Universe” of the German Science Foundation.

Appendix A The LTM model

We outline our revised version of the LTM formalism (Broadhurst et al., 2005; Zitrin et al., 2009, 2012a) and its implementation for constructing an initial lens model. This is the same formalism recently used to analyse some RELICS (Coe et al., 2019) galaxy clusters (e.g. Acebron et al., 2018). We assume that the mass distribution of each cluster member galaxy can be described by a surface-density profile, $\Sigma(r)=\Sigma_{q0}r^{-q}$, with r being the radius from the centre of the galaxy and $\Sigma_{q0}$ some normalisation factor. This profile can be integrated to give the interior mass, $M(<\theta)=\frac{2\pi\Sigma_{q0}}{2-q}(d_{l}\theta)^{2-q}$, which results in a deflection angle, due to a single galaxy, of:

[TABLE]

where $\theta$ is the angular distance from the galaxy’s centre, $G$ the gravitational constant, $c$ the speed of light, and $d_{l}$, $d_{s}$, $d_{ls}$ the angular diameter distances to the lens, to the source, and between the lens and the source. Simply put, $\alpha(\theta)\propto\theta^{1-q}$, per given lens and source redshifts (or angular diameter distances). Note that for all galaxies we assume circular symmetry, aside from one or two BCGs to which ellipticity can be incorporated. To these galaxies a core can also be added.

The deflection angle at a certain point $\vec{\theta}$ due to the lumpy galaxy components is simply a linear superposition of all galaxy contributions scaled by their luminosities, or fluxes $F_{i}$ (times some normalisation factor):

[TABLE]

In practice we use a discretised version of equation 8, over a 2D square grid $\vec{\theta}_{m}$ of $N\times N$ pixels, given by:

[TABLE]

[TABLE]

where $(\Delta x_{mi},\Delta y_{mi})$ is the displacement vector $\vec{\theta}_{m}-\vec{\theta}_{i}$ of the $m$th pixel point, with respect to the $i$th galaxy position $\vec{\theta}_{i}$. Note that there should be some normalisation factor in the expression of the deflections fields, which we initially set to unity.

In the code the sum is not performed explicitly as above, but instead is performed in Fourier space for speed-up purposes. The Fourier transform kernel is the same for all galaxies (substantially speeding up the calculation), and differs only for the independently modelled BCGs. The contribution of the BCGs can also be independently scaled.

The mass density distribution from the galaxy component is then given by the divergence of this deflection field (times half).

The second component of the model, the DM distribution, is assumed to be a Gaussian-smoothed version of the galaxy component. Therefore if $P$ represents the mass surface density (say, $\kappa$) of the smoothed map, then the deflection field from the DM component is simply:

[TABLE]

[TABLE]

where $P_{i}$ represents the (unnormalised) mass density value in the $i$th pixel of the smooth component as obtained by the smoothing procedure.

Note that in practice, instead of smoothing the galaxy-component mass density distribution directly, the convolution with a Gaussian kernel and the relevant summation to a deflection field are performed in one quick step in Fourier space.

The third component of the model is an external shear, given by:

[TABLE]

[TABLE]

where $\gamma$ is the strength of the external shear ($\gamma\geqslant 0$), and $\phi$ its position angle.

The total deflection field is obtained by combining the components as follows:

[TABLE]

where $N_{1}$ and $N_{2}$ are normalisation factors placing the galaxy and DM deflection fields onto the same (arbitrary) scale, $K_{gal}$ is the relative contribution of the galaxy component to the deflection field, and K is the overall normalisation. Note that the above procedure is typically performed with a $4\times 4$ lower resolution than the input HST images, and then the deflection fields are interpolated to the full HST resolution.

A MATLAB implementation of the code is available online333\colorbluehttps://github.com/adizitrin/initialLTM. Typically, in our various analyses with this implementation we find that $q\sim 1.3$, the Gaussian smoothing length is roughly $\sim 100-200$ pixels, and $K_{gal}\sim 0.05-0.1$. The external shear is about $\sim 0.1$ on average for elliptical lenses, aligned typically along the main axis of the cluster – which is usually also the BCG elongation direction. With these, it is easy to construct an initial lens model for any given cluster (see more details in 2.3). The only missing factor is the overall scaling. Since the overall normalisation is directly related to the size of the lens, or Einstein radius, a good initial guess can be such that a typical Einstein radius size is obtained (say $\sim 15-20\arcsec$). In our implementation we use $N_{1}$ and $N_{2}$ such that the mean of the absolute value of each deflection field is scaled to an arbitrary value of 200 (pixels), and then $K$ varies between $K\sim 0.5$ for the small lenses and $K\sim 3$ for large lenses. In addition, a scaling between the galaxy luminosities and the normalisation factor $K$, as was done by (Zitrin et al., 2012a) for their older-version code, remains for future work (see also Stapelberg et al. 2019).

Appendix B Parameters used to run MIFAL

(see Tables on next page)