The CluMPR Galaxy Cluster-Finding Algorithm and DESI Legacy Survey Galaxy Cluster Catalogue

Galaxy clusters enable unique opportunities to study cosmology, dark matter, galaxy evolution, and strongly-lensed transients. We here present a new cluster-finding algorithm, CluMPR (Clusters from Masses and Photometric Redshifts), that exploits photometric redshifts (photo-z's) as well as photometric stellar mass measurements. CluMPR uses a 2-dimensional binary search tree to search for overdensities of massive galaxies with similar redshifts on the sky and then probabilistically assigns cluster membership by accounting for photo-z uncertainties. We leverage the deep DESI Legacy Survey grzW1W2 imaging over one-third of the sky to create a catalogue of ~ 300,000 galaxy cluster candidates out to z = 1, including tabulations of member galaxies and estimates of each cluster's total stellar mass. Compared to other methods, CluMPR is particularly effective at identifying clusters at the high end of the redshift range considered (z = 0.75-1), with minimal contamination from low-mass groups. These characteristics make it ideal for identifying strongly lensed high-redshift supernovae and quasars that are powerful probes of cosmology, dark matter, and stellar astrophysics. As an example application of this cluster catalogue, we present a catalogue of candidate wide-angle strongly-lensed quasars in Appendix C. The five best candidates identified from this sample include two known lensed quasar systems and a possible changing-look lensed QSO with SDSS spectroscopy. All code and catalogues produced in this work are publicly available (see Data Availability).


INTRODUCTION
Galaxy clusters are the most massive gravitationally bound structures in the universe, and therefore they are of great use in cosmology and astrophysics.One motivation for the creation of catalogues of clusters is to locate likely areas for observing strongly-lensed time-variable phenomena such as gravitationally lensed supernovae and quasars (Meyers et al. 2009).The observation of the time delays between the lensed images of variable phenomena, coupled with proper modeling of the lensing system, can provide an independent constraint on the Hubble-Lemaitre parameter (Refsdal 1964).In 2015, multiple images of a supernova lensed by a galaxy cluster were observed (Kelly et al. 2015).Recently, the intermediate Palomar Transient Factory and its successor, the Zwicky Transient Facility (ZTF), have discovered several additional gravitationally lensed supernovae (Goobar et al. 2017;Goobar et al. 2023).Meanwhile, programs to observe ★ E-mail: michael.barth@umontreal.edugravitationally lensed quasars are ongoing (Wong et al. 2020;Napier et al. 2023).Further observations of gravitationally lensed transients will lead to improved constraints on the nature of dark matter and the expansion of the universe.
Galaxy cluster catalogues have many other applications.As one example, they can be used to study the impact of clusters on the evolution of galaxies in their vicinity (Dressler 1980).Gravitational magnification of background galaxies by galaxy clusters enables distant galaxies to be detected and studied at high resolution (e.g.Bayliss et al. 2014;Livermore et al. 2015;Johnson et al. 2017;Cornachione et al. 2018;Rivera-Thorsen et al. 2019;Ivison et al. 2020;Khullar et al. 2021;Bezanson et al. 2022).Furthermore, they can enable the detection and analysis of merging galaxy clusters, which are useful for constraining the dark matter scattering cross section (Wittman et al. 2023).Additionally, measurements of the apparent abundance of clusters can be a sensitive probe of cosmology (Haiman et al. 2001;Newman et al. 2002).
A purely photometric galaxy cluster-finding algorithm can take advantage of the large number of distant sources which can be observed using imaging.Conveniently, estimates of redshift and stellar mass can be obtained from the same photometric data.An algorithm which is optimally adapted for photometric observations will be essential for the upcoming generation of optical and infrared surveys at facilities such as Rubin Observatory (LSST Science Collaboration et al. 2009), Euclid Observatory (Euclid Collaboration et al. 2022), and the Nancy Grace Roman Telescope (Spergel et al. 2015), all of which will generate immense amounts of photometric data spanning many passbands.The redMaPPer algorithm (Rykoff et al. 2014) is an example of a well-tested photometric galaxy cluster-finding algorithm.The redMaPPer algorithm relies on a red-sequence model to detect candidate galaxy clusters.Currently redMaPPer has been applied to data from the Sloan Digital Sky Survey Data Release 8 (SDSS DR8) (Aihara et al. 2011) and Dark Energy Survey Science Verification (DES SV) (Dark Energy Survey Collaboration et al. 2016).The redMaP-Per SDSS DR8 catalogue contains 26,311 clusters over the redshift range 0.08<<0.55(Rykoff et al. 2014), while the redMaPPer DES SV contains 786 clusters from 0.2<<0.9(Rykoff et al. 2016).
Another well-tested cluster-finding algorithm was developed by Wen et al. (2009), henceforth WHL.This algorithm has been applied to SDSS (Wen et al. 2012;Wen & Han 2015), Subaru Hyper Suprime-Cam (data: Aihara et al. (2018) and clusters: Wen & Han (2021)), and the Dark Energy Survey (data: The Dark Energy Survey Collaboration (2005) and clusters: Wen & Han (2022)).The WHL algorithm, just like the algorithm described in this paper, uses the "nearest-neighbor" technique to find neighboring galaxies, but it uses a deterministic slice in redshift rather than the probabilistic approach developed in this work.
The current largest galaxy cluster catalogue (Zou et al. 2021) was created by applying the Clustering by Fast Search and Find of Density Peaks (CFSFDP) algorithm (Rodriguez & Laio 2014) to the Dark Energy Spectroscopic Instrument (DESI) Legacy Survey Data Release 8 (Dey et al. 2019).This work identified 540,432 galaxy cluster candidates reaching =1 in the Dark Energy Spectroscopic Instrument Legacy Survey.Zou et al. (2021) use Sunyaev-Zel'dovich and X-ray data to calibrate total galaxy masses from richness-mass relations for this sample.Similarly to WHL, (Zou et al. 2021) use a deterministic redshift slice, rather than a probabilistic approach.
In this paper we present a new galaxy cluster-finding algorithm, which we call CluMPR (Clusters from Masses and Photometric Redshifts).This algorithm is designed to minimize the impact of projection effects caused by uncertainties in the photometric redshifts (photo-'s) of individual galaxies.In this work, we optimize the details of this algorithm to enable the selection of a high-purity sample of high-stellar-mass clusters.These choices aid in selecting those galaxy cluster candidates which are most likely to cause strong lensing and help to optimize searches for transients associated with quiescent galaxies, while simultaneously minimizing contamination of the sample by lower-mass galaxy groups.Our algorithm determines total cluster stellar masses probabilistically by incorporating photometric redshift uncertainties and compensating for redshiftdependent mass incompleteness.We also provide a catalogue of cluster member galaxies and a catalogue of candidate gravitationally lensed quasars.Our catalogues are intended to complement existing galaxy cluster catalogues such as the ones created by Zou et al. (2021) and Yang et al. (2021) using DESI Legacy Survey data, and can also be used to compare the results of applying different algorithms to similar data.
This paper is organised as follows: In Section 2, we describe the data products we used to generate our galaxy cluster catalogue.Sec-tion 3.1 describes the CluMPR algorithm in detail.In Section 3.2, we describe our process for calibrating galaxy cluster parameters and assigning galaxy cluster member galaxies.In Section 4, we provide some summary graphs and data summarizing the characteristics of the CluMPR galaxy cluster catalogue.Finally, in Section 5 we summarize the results, applications, and potential future extensions of our work.The Data Availability section provides the publicly accessible locations of all data and code related to this paper, including the cluster catalogues and candidate lensed quasar catalogues.Appendices A and B describe details regarding the calibration of our algorithm; Appendix C describes a set of candidate gravitationally-lensed quasars we have identified; and Appendix D describes the columns of all catalogues assembled in this work.Throughout this paper, we assume a cosmology with  0 = 70 // , Ω  = 0.286, Ω  = 0.714 (Bennett et al. 2014).

DATA
The DESI Legacy Imaging Surveys (Dey et al. 2019) has produced multi-band images of over a third of the sky (roughly 14,000 deg 2 ) in three optical bands (, , ) and two infrared bands from WISE (1, 2) (Wright et al. 2010).The optical component of this survey is composed of three separate programs: the Beĳing-Arizona Sky Survey (BASS) (Zou et al. 2017), the Dark Energy Camera Legacy Survey (DECaLS) (Flaugher et al. 2015;Dey et al. 2019) with supplementary data from the Dark Energy Survey (The Dark Energy Survey Collaboration 2005), and the Mayall z-band legacy survey (MzLS) (Silva et al. 2016).A full description of the Legacy Imaging Surveys is available at Dey et al. (2019).The primary purpose of the DESI Legacy Imaging Surveys is to obtain targets for spectroscopic follow-up with DESI (Dark Energy Spectroscopic Instrument, DESI Collaboration et al. (2022)).We obtained photometric redshifts for our analysis from Zhou et al. (2021) and stellar mass estimates from Zhou et al. (2023).
We apply several photometric cuts to select galaxies for our analysis with good measurements of photometric redshifts and stellar masses and with as little contamination from stars as possible.We restrict ourselves to objects with (i) z magnitude < 21 (photo-z's are only reliable for z band magnitude < 21, see Zhou et al. (2021)), (ii) Either Gaia photometric mean magnitude > 19, Gaia astrometric excess noise > 10 0.5 , or Gaia astrometric excess noise = 0 (in order to remove stars), (iii) TYPE != PSF (setting TYPE to exclude point spread function objects removes stars and quasars while keeping the vast majority of galaxies, since bright galaxies are able to be resolved) (iv) Photometric redshift > 0.01 (due to photo-z and de-blending unreliability at low z) We restrict ourselves to galaxies with a known stellar mass from Zhou et al. (2023).These are objects which have valid photometry in the , , , W1, W2 bands and satisfy the stellar rejection cut of  − 1 > 1.75 × ( − ) − 1.1.We limit ourselves to galaxies with a DESI legacy Survey DR9 bitmask value 0 or 4096, which correspond respectively to unmasked galaxies and SGA galaxies (Moustakas et al. 2021).The photometric redshifts and the stellar masses were estimated using the Random Forest algorithm trained on legacy survey photometry and a carefully vetted training set of spectroscopic redshifts and stellar population synthesis-based mass estimates.The errors on the photometric redshifts were generated by perturbing the original photometry and using the spread in the distributions of the predicted redshifts as a measure of the error.Detailed descriptions of the methods and catalogues can be found in Zhou et al. (2021Zhou et al. ( , 2023)).
Note: throughout this text, a "sweep" refers to a continuous subsection of DESI Legacy survey photometry data which is contained in a single file.

CluMPR Algorithm Overview
The CluMPR algorithm is comprised of four essential steps: (i) Potential Cluster Center Selection: Select galaxies with large stellar mass.
(ii) Membership Probability Calculation: For each potential cluster center galaxy, evaluate the expected number of neighbouring galaxies within a cylinder oriented lengthwise along the line of sight () direction (see Fig. 1).In the direction of the cylinder's radius, membership is determined deterministically using a binary search tree.Along the z direction, membership is determined probabilistically based on photo- uncertainties.
(iii) Richness Thresholds: Keep potential cluster center galaxies with a significant excess of neighbouring galaxies (this is equivalent to defining a galaxy cluster by setting a minimum richness).
(iv) Aggregation: neighbouring potential cluster center galaxies are aggregated and the galaxy with the highest local stellar mass is chosen as the cluster center.This cylinder-based method is a probabilistic variant of the widelyused (deterministic) "counts-in-cylinders" technique (Hogg et al. 2004;Blanton et al. 2006;Kauffmann et al. 2004;Barton et al. 2007;Berrier et al. 2011).The individual steps of the CluMPR algorithm are explained in greater detail below.

Potential Cluster Center Selection
Potential cluster center galaxies are chosen from the general pool of galaxies based on a minimum stellar mass threshold.We set this minimum stellar mass threshold to be log(stellar mass) > 11.2 [log( ⊙ )].
This mass threshold was chosen based on visual inspection of galaxy clusters (we chose to minimize the fraction of low-mass groups and spurious detections in order to create a high-purity cluster catalogue, but this mass threshold can be lowered if a lower-purity catalogue is desired).

Membership Probability Calculation
For each potential cluster center galaxy chosen in 3.1.1,we count the number of neighbours in a cylinder oriented lengthwise along the  direction and centered on the potential cluster center galaxy.
First, galaxies are selected within a fixed physical radius of the cluster center galaxy (we use three radii: 1 Mpc, 0.5 Mpc, and 0.1 Mpc.The significance of these radii is explained in the subsequent sections).This first step is accomplished using a binary search tree; our implementation uses the cKDTree algorithm from the SciPy Python package (Virtanen et al. 2020).In order to avoid propagating photo- errors into the angular positions of galaxies, we perform this step using angular coordinates.The physical radius is calculated in angular coordinates using the following equation: where R ang is the radius in angular coordinates (radians), R phys is the physical radius in megaparsecs,  is the redshift, and  () is the comoving distance given by where  ( ′ ) is the Hubble-Lemaitre parameter.A graph of Equation 1for  ℎ = 1  is shown in Figure 2.For each potential cluster center galaxy, we calculate the radius using Equation 1 and the photo- of the cluster center galaxy.To ensure equal area projection of angular coordinates, we use a sinusoidal (Sanson-Flamsteed) projection.To avoid issues associated with very small values for angular coordinates, we add an offset of 50 radians to each coordinate.Next, to determine cluster membership along the z-direction (lengthwise along the cylinder), we apply a redshift cut based on the photo- of the cluster center galaxy.The length of this cut (the half-length of the cylinder) is given by where   is the interpolated redshift-dependent binned average of the photo- error for all galaxies.To train our model for   , we used 100 bins between the lowest and highest photo-'s in the training sample; we limited   to 0.1 at high z.We used 20 randomly-selected DESI Legacy Survey sweep files (10 from North (BASS/MzLS) and 10 from South (DECaLS)) as our training samples.The linearlyinterpolated trained results for  , ℎ and  , ℎ are shown in Figure 3.   is significantly larger than the true redshift extent of a galaxy cluster at all redshifts, which means that a probabilistic approach to cluster membership along the z direction is required.We use   to determine the probability that a galaxy within the search cylinder's radius is a true member of the cluster (that is, the probability that it is truly within the cylinder).In order to determine this probability, we use statistical resampling, whereby we create many realizations of each galaxy's potential true redshift; these realizations are drawn from a Gaussian distribution centered on the photo- of the galaxy and with a standard deviation equal to the photo- error for the galaxy.The percentage of these redshift realizations falling within   of the cluster center galaxy's photo- is the membership probability.

Richness Thresholds
In order to identify galaxy clusters among the potential cluster center galaxies selected in Step 3.1.1,we apply a threshold for the minimum number of galaxies that must be located within the galaxy's search radius.This is equivalent to defining what is a galaxy cluster by setting a minimum richness threshold.We apply a redshift-dependent threshold to both the 1-Mpc cylinder and the 0.5-Mpc cylinder; only potential cluster center galaxies that satisfy both thresholds are kept as galaxy clusters.For the 1-Mpc cylinder, the threshold we used is and for the 0.5-Mpc cylinder, the threshold we used is where  is the redshift-binned average number of galaxies within a cylinder of that size, and   is the redshift-binned standard deviation of the number of galaxies within a cylinder of that size.These thresholds were chosen based on visual inspection for purity of resulting galaxy clusters.We assume a Poisson distribution for the number of galaxies within a cylinder, which allows us to define The richness threshold is interpolated linearly as a function of ; the results of this interpolation are shown in Figure 4.The training data for these thresholds are the same as in Section 3.1.2,and the redshift bins are =0.01apart with a width of  = 0.05, ranging from the lowest to the highest photo- in the training samples.

Aggregation
In order to remove duplicate galaxy clusters (several center galaxy candidates for one galaxy cluster), we use a similar method to Sec-tion 3.1.2,except that the radius of the aggregation cylinder is 1.5 Mpc, and the   cutoff is not done probabilistically.First, all cluster central galaxies are sorted by (uncorrected, see 3.2.2) total stellar mass within a 0.5 Mpc radius.Beginning with the cluster central galaxies with the greatest 0.5 Mpc stellar mass, neighbouring cluster central galaxies within R=1.5 Mpc and   are labelled as belonging to the same cluster.The galaxy with the greatest (uncorrected) total stellar mass within a 0.1 Mpc radius is then chosen as the cluster center galaxy.

Data Cleaning
We apply several flags to our galaxy clusters to remove poor objects.First, galaxy clusters within 0.285 degrees of the edge of the DESI Legacy Survey footprint are flagged.We flag clusters near large foreground galaxies from the Siena Galaxy Atlas (SGA) (Moustakas et al. 2021).Our criteria for the SGA galaxies are either (i) D26 > 1.5, (ii) R_MAG_SB24 < 13.5, R_MAG_SB24>0, and Z_LEDA <0.05Some of these flagged objects are false detections (which are mostly due to failed de-blending of foreground galaxies), while others are real galaxy clusters with significant foreground contamination which could affect the photo- and stellar mass estimates.Finally, we remove all clusters located in the isolated strips of the DESI Legacy Survey footprint with declination < -10.3 degrees and 110 < Right ascension (degrees) < 250.We provide a catalogue which includes all flagged objects as a supplement to the main CluMPR catalogue.

Cluster Parameters
Each galaxy cluster in our catalogue is labelled with a unique "gid", which is created from the cluster center galaxy's coordinates using the following formula: round(RA, 6 decimals) * 10 16 + round(DEC + 90, 6 decimals) * 10 6 where RA is right ascension and DEC is declination (both in degrees).Each galaxy cluster has a richness estimate (number of neighbours) for the 1 Mpc, 0.5 Mpc, and 0.1 Mpc search radii.The "gid" is necessary to identify the galaxy cluster members described in Section 3.2.3.Additionally, we provide galaxy cluster parameters for redshift and mass as described below.

Probability and Mass Weighting of Cluster Redshift Parameters
For each galaxy cluster we calculate the following redshift parameters (in addition to providing the photo- information for the cluster center galaxy): (i) Average photo- of all galaxy cluster member galaxies (ii) Average photo- of all galaxy cluster member galaxies (weighted by membership probability) (iii) Average photo- of all galaxy cluster member galaxies (weighted by galaxy stellar mass times membership probability) (iv) Standard error of photo- of all galaxy cluster member galaxies (v) Standard error of photo- of all galaxy cluster member galaxies (weighted by membership probability) (vi) Standard error of photo- of all galaxy cluster member galaxies (weighted by galaxy stellar mass times membership probability) The average photo- estimates for the clusters are biased high, particularly for certain galaxy clusters at low z.This is due to systematic errors caused by the higher probability of a line-of-sight galaxy at high z falling within the cylinder of the cluster (as high-z galaxies have larger photo- errors on average).This bias is further amplified by volume effects: the angular diameter of the cylinder encloses a larger volume at high z.In this work, we do not attempt to correct for this systematic error.For this reason, we recommend using the photo- of the cluster center galaxy as a relatively unbiased estimate of the galaxy cluster redshift.Discrepancies between different estimates for the cluster's average photo- are useful for identifying individual clusters with significant line-of-sight effects which may affect cluster parameters.The standard error estimates are (to first order) independent of the systematic error described above, so they can be used as an approximation for the random error uncertainty in cluster redshift.If more robust estimates (such as Hodges-Lehmann) for cluster redshift and uncertainty are desired, they can be calculated from the cluster membership catalogues (see Section 3.2.3).

Cluster Richness and Stellar Mass Estimates
We provide a background-corrected richness estimate which is background-corrected but not corrected for survey incompleteness.This richness is calculated by summing the membership probabilities of the galaxy cluster member galaxies and then subtracting the average galaxy background from this sum.The calculation of the average galaxy background is described in Appendix A.
We determine galaxy cluster stellar masses by summing the membership probability-weighted stellar masses of all cluster member galaxies within the 1 Mpc, 0.5 Mpc, and 0.1 Mpc search radii.We exclude galaxies with a mass below a threshold described in Appendix B (Equation B3).We then subtract the average galaxy stellar mass background, which is described in Appendix A. Finally, we apply a redshift-dependent correction factor which compensates for redshift-dependent mass incompleteness.This factor is described in Appendix B; the result is an incompleteness-corrected cluster stellar mass containing all cluster members above log( ★ / ⊙ ) = 10.

Cluster Membership
For every galaxy cluster, we have compiled a list of member galaxies (galaxies whose membership probability > 0).For galaxies which are found to belong to more than one galaxy cluster, we assign the galaxy cluster for which there is the highest probability of membership.For each galaxy cluster, the cluster member galaxies can be identified using the "gid" parameter (see Section 3.2).The membership probability for the cluster member galaxies provided in the membership catalogue is not the same as the probability defined in Section 3.1.2.Instead, it is the probability of the member galaxy having a true redshift equal to the photo- of the galaxy cluster's central galaxy, assuming a Gaussian distribution centered on the member galaxy's redshift with a standard deviation equal to that member galaxy's photo- error.

RESULTS
The CluMPR algorithm has returned 309115 candidate galaxy clusters from 0.1<<1.Galaxy clusters with <0.1 are excluded as our algorithm is not optimised for those objects, which are frequently contaminated by poorly deblended foreground galaxies.Figure 5 shows a few examples of galaxy clusters from our catalogue at various redshifts.The distribution of our catalogue's galaxy clusters on the sky is shown in Figure 6. Figure 7 shows the redshift and stellar mass distribution of CluMPR clusters.For 0.1<<0.3, the sample is volume-limited; above =0.3,the sample is limited by incompleteness.The stellar mass correction factor causes the stellar mass distribution to remain relatively flat across z.Low-mass galaxy clusters are excluded by our thresholds and incompleteness at high z, while high mass galaxy clusters are volume-limited.The rise in high-stellar-mass clusters at <0.3 indicates that the catalogue is volume-limited for those redshifts; the decline in high-stellar-mass clusters at z > 0.6 may be due to limits imposed by structure formation (relations between our stellar mass estimates and total halo mass must be established to verify this cause).

Cross-matching with Zou et al. 2021
The galaxy cluster catalogue described in Zou et al. (2021) was compiled using the same DESI Legacy Survey data as our catalogue, so cross-matching the catalogues yields a direct comparison of our implementation of CluMPR with the implementation of the CFSFDP  algorithm by Zou et al. (2021).For cross-matching, we apply a simple 2-dimensional cross-matching which ignores line-of-sight coincidence.The radius for determining a match was 5 arcminutes, which is sufficient to handle most differences in the choice of cluster center galaxy while minimizing spurious matches.We found that this radius worked reasonably well at all redshifts.Overall, Zou et al. (2021) found 85.4% of the galaxy clusters identified by CluMPR above =0.1.Figure 8 shows the stellar mass   2021) missing these objects is that they use a minimum richness of 10 which is not corrected for incompleteness; at high redshift, many clusters will have a low observed richness.This effect is best illustrated by Figure 9, which shows the 1-Mpc CluMPR richness as a function of redshift; color indicates the fraction of CluMPR galaxies found by Zou et al. (2021).At low-, the CluMPR richness threshold for a cluster is relatively high, but it drops rapidly with ; by  0.75, the richness threshold falls below 10, and the fraction of galaxy clusters recovered by Zou et al. ( 2021) falls to 0.5.CluMPR is able to maintain relatively high purity at high z despite the low richness threshold by using high-stellar mass galaxies as a reliable starting point for finding galaxy clusters, which works even in the absence of a large number of cluster member galaxies.
Overall, CluMPR finds 58.9% of the galaxy clusters identified by

Cross-matching with SDSS redMaPPer and SDSS spectroscopy
Both CluMPR and Zou recover most of the objects from the SDSS redMaPPer cluster catalogue.Using the 5 arcminute radius cross-matching described in Section 4.1, CluMPR recovers 95% of redMaPPer clusters, while Zou et al. (2021) recover 99% of redMaP-Per clusters.Figure 11 shows the 1-Mpc cluster stellar mass vs redMaPPer richness.As redMaPPer richness is known to be correlated with total cluster mass (Simet et al. 2017), the observed correltion between CluMPR stellar mass and redMaPPer richness is a strong indicator that CluMPR stellar masses are correlated with total cluster mass.We provide spectroscopic redshifts for galaxies whose spectra Figure 12 shows the residuals between SDSS spectra and three photometric estimates for cluster redshift from CluMPR, as well as the residuals for redshift estimates by redMaPPer and Zou et al. (2021) The bias for average redshifts obtained by CluMPR is caused by line-of-sight galaxies falling within the cluster; this bias is increased when photo- errors are systematically over-estimated, as this creates an artificially long cylinder and increases the likelihood for galaxies to fall within the cylinder.Zhou et al. (2021) indicate that photo- errors are indeed overestimated, so this bias would be reduced with an improved redshift uncertainty quantification.However, our algorithm is still able to reduce the bias by weighting the cluster members by mass*probability (shown, middle right) or just probability (not shown).This indicates that our stellar mass and richness estimates are likely also improved by our probabilistic approach.

CONCLUSIONS
We have used the CluMPR algorithm to identify over 300,000 galaxy clusters using the DESI Legacy Survey imaging data.This algorithm provides incompleteness-corrected stellar masses for each cluster as well as cluster membership probabilities for potential member galaxies.Our catalogue shows good correspondence to previous catalogues such as Zou et al. (2021) and SDSS redMaPPer (Rykoff et al. 2014) in overlapping regimes, and visual inspection of the detected galaxy clusters indicates that our methods have produced a high-purity sample.CluMPR has several advantages over non-probabilistic clusterfinding algorithms such as the WHL algorithm (Wen et al. 2009)  We include a catalogue of the member galaxies for all clusters, which can be used to calculate alternative estimates of the cluster redshift.All plotted redshift estimates underestimate spectroscopic redshift near the high- limit due to sharply decreasing galaxy counts at increasing redshifts (attenuation bias, see Freeman et al. 2009), which leads to a redshift estimate for the cluster central galaxy which is biased low.
et al. (2021).As an example, for strong lensing searches, CluMPR enables us to select a sample of high-stellar-mass clusters that simultaneously has high purity.Higher-redshift lensing searches are further enhanced as CluMPR is able to more easily identify galaxy clusters at higher  ( > 0.75), a regime where smaller number of cluster galaxies will be brighter than magnitude thresholds and hence apparent richnesses will be lower.However, by using the presence of a single massive galaxy (( ★ ) > 11.2) as a signpost to reliably anchor cluster assignments, CluMPR is able to still identify clusters with high purity even when the apparent richness drops.In contrast to red-sequence-based algorithms such as redMaPPer (Rykoff et al. 2014), CluMPR includes all potential galaxy cluster members in its search, which leads to a more complete member sample.This difference will become more relevant for deeper surveys, as we would expect a higher fraction of cluster members which are not on the red sequence at high  (Schechter 1976;Dahle et al. 2013).Finally, our probabilistic approach to cluster membership reduces systematic errors in cluster properties caused by the incorporation of infalling foreground or background galaxies.The CluMPR DESI Legacy Survey catalogue can now be used by astronomers to search for strongly-lensed time-variable phenomena (such as supernovae and quasars) as well as other applications.As an example of how this catalogue may be used, we provide a supplementary catalogue of candidate lensed quasars in Appendix C. The five best candidates identified from this sample include two known systems and a possible changing-look lensed QSO with SDSS spectroscopy.The CluMPR DESI Legacy Survey catalogue can also be used for galaxy evolution studies and to search for merging galaxy clusters.
The CluMPR sample could be made more useful for cosmological studies by calibrating the relationship between the cluster stellar masses we have measured and total cluster masses as derived from Sunyaev-Zel'dovich, X-ray, or weak lensing measurements.Future datasets from Rubin Observatory (LSST Science Collaboration et al. 2009), Euclid Observatory (Euclid Collaboration et al. 2022), the Nancy Grace Roman Telescope (Spergel et al. 2015), and other nextgeneration surveys will provide another application for the CluMPR algorithm.These projects should provide higher-precision photometric redshifts and useful stellar mass measurements for objects several magnitudes fainter than the Legacy Survey objects studied here.As a result, these deep imaging surveys could extend CluMPR's reach to  = 2 and beyond.ACKNOWLEDGEMENTS M. J. Yantovski-Barth would like to acknowledge his PhD advisors Yashar Hezaveh and Laurence Perreault-Levasseur for their patience and advice as he finalised the results presented in this paper.Additionally, we thank Gourav Khullar for a helpful discussion of cluster-lensed quasars.The efforts of M. J. Yantovski-Barth were supported by three grants from the NASA Pennsylvania Space Grant Consortium in 2020 and 2021, as well as by the bourse d'excellence du centenaire from the Fonds du Centenaire of the Department of Physics at the Université de Montréal.The efforts of J. A. Newman were supported by grant DE-SC0007914 from the U.S. Department of Energy Office of Science, Office of High Energy Physics.
The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beĳing-Arizona Sky Survey (BASS; NOAO Prop.ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop.ID #2016A-0453; PI: Arjun Dey).DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF's NOIR-Lab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab.The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.
NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration.Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Minis-terio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l'Espai (IEEC/CSIC), the Institut de Fisica d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF's NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.
BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program "The Emergence of Cosmological Structures" Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance.The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 11433005).
The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology.NEOWISE is funded by the National Aeronautics and Space Administration. The

DATA AVAILABILITY
The CluMPR galaxy cluster catalogue and cluster member galaxy catalogue are currently publicly available at DOI:10.5281/zenodo.8157752.The candidate lensed quasar catalogue is publicly available atDOI:10.5281/zenodo.8164870.Descriptions of each catalogue's columns can be found in Appendix D.
These catalogues are also available to NERSC users at global/cfs/cdirs/desi/users/mjyb16/CLUMPR_DESI_2023.The main catalogue can be found in the folder clean_data, the extended catalogue including flagged objects can be found in the folder unclean_data, catalogues of cluster member galaxies can be found in the folder members, and the candidate lensed quasar catalogues can be found in the folder lensed_qso_candidates.
Code from this project can be found on GitHub. is the number density of galaxies with a mass , and Φ ★ is the normalization factor.We fit this model to the high-mass end of the mass distribution histograms of galaxies observed to be in a cluster.We are able to use this method since we observe the histograms at all redshifts to roughly follow the Schechter function at high masses, though the number count approaches 0 for low mass galaxies (an example of a stellar mass histogram is shown in Figure B2).As the first step of building our model, we use the Schechter function to approximate the true mass distribution of galaxies, given by the double Schechter function (Davidzon et al. 2017) . This approximation is necessary in order to reduce the number of parameters to be fitted.A demonstration of the sufficiency of a single Schechter function to fit to a double Schechter function above () = 10 is shown in Figure B1.We fit  in our model to the double Schechter for the values shown in Table B1. ★ and Φ() in our model are fitted to the stellar mass distribution of galaxies in our clusters.
We fit  ★ and Φ ★ to the high-mass end of the galaxy mass distribution histogram for several redshift bins (as an example, see Figure B2).We use 22 redshift bins between 0.025<<1.025with bin width =0.05.The high-mass ends of each histogram are located by fitting a piecewise linear function of the form to the peaks of the histograms, where z is redshift bin.The maximum limit of (  ) = 11.2 is required in order to ensure that the mass of the cluster central galaxy is not excluded from the mass calculation in the main algorithm (see section 3.2).We find the bestfit values to be a = 1.36202 and b = 9.96855.This model, along with an example of actual histogram peak locations for a small sample of galaxies, is shown in Figure B3.
The contribution of missing low-mass galaxies to the total stellar mass is determined by taking the ratio between the integral of Equation B1 from () = 10 to infinity and the integral of Equation B1 from the mass limit B3 to infinity.The latter integral represents the incomplete stellar mass calculated by the main algorithm.The ratio between the integrals is thus where   is the total stellar mass above () = 10 and  is the incomplete (observed) stellar mass.Multiplying the observed galaxy cluster stellar masses by this redshift-dependent mass correction factor allows us to estimate the true stellar masses above () = 10 within the 1 Mpc, 0.5 Mpc, and 0.1 Mpc search radii.We find that a quadratic equation is a good fit for the observed log(stellar mass) correction factor as a function of redshift: where a = 0.19347, b = 0.23392, and c = 0.00694.A plot of the stellar mass correction factor as a function of redshift is shown in Figure B4; we trained this model on a dataset consisting of 150 sweep files randomly chosen from DESI Legacy Survey North and South.The stellar mass correction factor is provided for all galaxy clusters in our catalogue so that the incomplete stellar mass can be easily retrieved.

APPENDIX C: GRAVITATIONALLY LENSED QUASAR CANDIDATES
As a supplement to our galaxy cluster catalogue and galaxy cluster member catalogue, we provide a catalogue of candidate gravitationally lensed quasars.Gravitationally lensed quasars are valuable as probes of the Hubble parameter for the expansion of the universe (Suyu et al. 2017).As data, we used quasars from the DESI Bright and Dark survey target selection (Myers et al. 2023).We compiled the candidate lensed quasar catalogue by counting the number of quasars within the Einstein radius of each galaxy cluster identified by CluMPR.We use two Einstein radii: one corresponds to wide angle lensing by the entire cluster ( = 10 15  ⊙ ), and the other corresponds to lensing by the core of the cluster ( = 0.25 * 10 15  ⊙ ).Since there is contamination of QSO targets by stars, particularly in dense star fields, we remove galaxy clusters which have counts within 3 Einstein radii >= 6 times the counts within one Einstein radius.We  B4) is the ratio between two integrals along the green line (Schechter function): the total mass is the integral above log(mass) = 10, and the incomplete mass is the integral above the peak of the histogram (in this case, above log(mass) = 11).This ratio is designed to compensate for the missing galaxy counts for galaxies with masses below log(mass) = 11.rate candidates based on how many of the quasars have similar colors in bands g-r, g-z, and r-W1.If at least two quasars are within 1 magnitude of each other in all 3 colors, we rate the candidate at Grade C. If at least three quasars are within 1 magnitude of each other in all 3 colors, we rate the the candidate at Grade B. If either 4 quasars or two combinations of 3 quasars are within 1 magnitude of each other in all 3 colors, we rate the candidate as Grade A. Spectroscopic follow-up of the lensed quasar candidates and cluster lens modeling will enable confirmation of true lensed quasars.

C1 Recovery of known cluster-lensed quasars
To date, six widely-separated gravitationally lensed quasars (corresponding to a cluster-scale lens mass) have been discovered (Inada et al. 2003(Inada et al. , 2006;;Dahle et al. 2013;Shu et al. 2018Shu et al. , 2019;;Martinez et al. 2023).Of the six known cluster-lensed quasars, we recover three: SDSS J1004+4112 (Grade A core-lensed), SDSS J1029+2623 (Grade C core-lensed), and COOL J0542-2125 (Grade A corelensed).We do not recover SDSSJ2222+2745 or SDSS J1326+4806 due to only one of the quasar images being in the DESI target catalogues as a QSO, and we do not recover SDSS J0909+4449 due to the lensing cluster being absent from both the CluMPR and Zou et al. ( 2021) catalogues (it is more accurately described as a galaxy group).B4) as a function of redshift.Mass incompleteness increases with redshift (higher mass limit), so the correction factor increases correspondingly.

C2 New Grade A core-lensed candidates
Of the five core-lensed candidates ranked Grade A, two are known lensed quasars.We discuss the remaining Grade A core-lensed candidates below.

C2.1 CluMPR J1278885160133433856: A candidate lensed changing-look quasar
The candidate lensed quasars located near the cluster CluMPR J1278885160133433856 are shown in Figure C1.The galaxy cluster has an SDSS spectroscopic redshift of =0.525, and the cluster center galaxy is located at RA = 127.8885,DEC = 43.4339(J2000).Two of the candidate lensed images, QSO 1 and QSO 2, are confirmed quasars with SDSS redshifts of  = 1.3056 ± 0.0002 and  = 1.3021 ± 0.0004, respectively.The spectra of these quasars are shown in Figure C2, with zoomed-in views of the C III] and Mg II broad emission lines (Figures C3 and C4,respectively).The alignment between emission lines in the two spectra indicates that QSO 1 and QSO 2 are at the same or nearly the same redshift ( ≈ 1.3), but the noise in the spectra and the absence of narrow emisison lines makes evaluating the accuracy of the assigned SDSS spectroscopic redshifts challenging.There are notable differences between the spectra: QSO 2 appears brighter and has a redder color (negative slope) compared to QSO 1, and the broad emission lines appear more prominent in QSO 2 than in QSO 1.We offer two potential explanations for the differences in the spectra: either these are two distinct quasars that happen to be at very similar redshifts, or these are two images of the same highly-variable ("changing-look") quasar.Changing-look quasars are quasars whose luminosity and colors (i.e., continuum slopes) change quickly by a large factor over the course of hundreds of days to several years (e.g.Bian et al. 2012;MacLeod et al. 2016MacLeod et al. , 2019)).The luminosity changes are accompanied by the appearance and disappearance of broad emission lines (e.g.LaMassa et al. 2015;Runnoe et al. 2016;Ruan et al. 2016).Changing-look quasars have recently been discovered at redshifts beyond 1 (Ross et al. 2020;Guo et al. 2020).The variability timescale of changing look quasars is on the same order as the time delays between lensed images for a cluster-scale lens, so it is possible that the in spectra between QSO 1 and QSO 2 is due to a gravitational time delay.
To evaluate the changing look quasar interpretation further, we consider the following: • According to the formal uncertainties, the redshifts of the two quasars (∼ 30-50 km s −1 ) differ by ≈ 8.However, the detailed discussion of redshift determination methods used by the SDSS by Lyke et al. (2020, see their Section 4.6 and Figure A1) suggests a much larger redshift uncertainty of ∼ 300 km s −1 .Viewed in this light, the quasar redshifts differ only by 1.5.
• We measured the rest-frame equivalent widths of the broad emission lines and found them to be EW(Mg II) = 34 ± 2 Å and EW(C III]) = 33 ± 3 Å in QSO 1 and EW(Mg II) = 26 ± 1 Å and EW(C III]) = 25 ± 1 Å in QSO 2. These measurements reflect the Baldwin effect (anti-correlation between EW and luminosity; e.g., Green et al. 2001, and references therein) and follow the behavior seen in extremely variable quasars by Ren et al. (2022).
• A caveat to the above is that, the relative luminosities of the two quasars cannot be readily evaluated as the observed brightness of quasar images depends on the magnification of the lens (Bauer et al. 2011), so the true luminosity of the quasar images cannot be estimated without a well-constrained mass model for the lens.The above considerations suggest that the the lensed changinglook quasar hypothesis remains plausible for QSO 1 and QSO 2. To test this hypothesis further we suggest the following.
• Improving the redshift determination.-This can be done by obtaining spectra with higher signal-to-noise ratio in the near-IR or sub-mm bands to detect narrow emission lines.Spectra in the wavelength range 7800-8600 Å can catch the narrow [Ne V] 3346,3426 and [O II] 3726,3729 doublets.The wavelengths of these lines and their relative strengths (the [Ne V]/[O II] ratio) can test the hypothesis that QSO 1 and QSO 2 are lensed images of the same quasar.In the same spirit, spectroscopy in the J band can catch the H line and [O III] 4959,5007 doublet and allow a similar test.Alternatively, spectroscopy with ALMA may detect the CO(2-1) line from the quasar host galaxies (rest frequency of 230.5 GHz) that will allow for an exquisite redshift determination.
• Monitoring for spectroscopic variability.-If the spectra of QSO 1 and QSO 2 do indeed reflect a change in the state of a single quasar, we should expect that one of the two spectra will change to match the other in the next few years.Therefore, the hypothesis can be tested by repeated spectroscopic observations of the two quasars.
There are two more possible quasar images (C1 and C2) which have been selected for DESI spectroscopy.Finally, there is a dim object resembling a lensing arc.Further spectroscopy and deeper imaging are necessary to determine whether these objects are related to QSO 1 and QSO 2.

C2.2 CluMPR J2439482900096790528
The candidate lensed quasars located near the cluster CluMPR J2439482900096790528 are shown in Figure C5.The galaxy cluster's central galaxy is located at the photometric redshift  = 0.333 with RA = 243.9483,DEC = 6.7904.There are no publicly available spectra for the candidate quasars at this time, but they are part of the DESI dark time targets.Spectroscopy and deeper imaging are necessary to determine whether this system is a quasar-lensing cluster.

C2.3 CluMPR J2420498300113853952
The candidate lensed quasars located near the cluster CluMPR J2420498300113853952 are shown in Figure C6.The galaxy cluster's central galaxy is located at the photometric redshift  = 0.750 with RA = 242.0498,DEC = 23.8542.One of the candidate quasars (QSO 1) has SDSS spectroscopy with a spectroscopic redshift of 2.205.There are no publicly available spectra for the other candidate quasars at this time, but they are part of the DESI dark time targets.Spectroscopy and deeper imaging are necessary to determine whether this system is a quasar-lensing cluster.

Figure C6
. Candidate gravitationally lensed quasar system near galaxy cluster CluMPR J2420498300113853952 (cluster center galaxy in green).Objects in white circles are QSO targets for DESI, and there is a possibility that some of them are images of the same gravitationally lensed quasar (it is relatively unlikely that they are all the same quasar due to their configuration).QSO 1 is a confirmed quasar with SDSS spectroscopic redshift of 2.205.

APPENDIX D: CATALOGUE DESCRIPTIONS
Below we provide descriptions of the column names and contents for the CluMPR DESI Legacy Survey catalogue, associated cluster member catalogue, and the candidate gravitationally lensed quasar catalogues.The public location of the catalogues can be found in the Data Availability) section.Log(stellar mass) of the cluster within a 1-Mpc physical radius, corrected for incompleteness and galaxy background cluster_mass_halfmpc Log(stellar mass) of the cluster within a 0.5-Mpc physical radius, corrected for incompleteness and galaxy background cluster_mass_tenthmpc Log(stellar mass) of the cluster within a 0.1-Mpc physical radius, corrected for incompleteness and galaxy background richness_onempc Expectation number of observed galaxies in cluster within a 1-Mpc physical radius, corrected for galaxy background richness_halfmpc Expectation number of observed galaxies in cluster within 0.5-Mpc physical radius, corrected for galaxy background richness_tenthmpc Expectation number of observed galaxies in cluster within 0.1-Mpc physical radius, corrected for galaxy background Log(stellar mass) of the cluster within a 1-Mpc physical radius, corrected for incompleteness and galaxy background cluster_mass_halfmpc Log(stellar mass) of the cluster within a 0.5-Mpc physical radius, corrected for incompleteness and galaxy background cluster_mass_tenthmpc Log(stellar mass) of the cluster within a 0.1-Mpc physical radius, corrected for incompleteness and galaxy background mass_bkgd_onempc Mass background (which has been subtracted from cluster_mass_onempc) for one-mpc radius mass_bkgd_halfmpc Mass background (which has been subtracted from cluster_mass_onempc) for one-mpc radius mass_bkgd_tenthmpc Mass background (which has been subtracted from cluster_mass_onempc) for one-mpc radius correction_factor Correction factor for mass incompleteness neighbours_onempc Expectation number of observed galaxies in cluster within a 1-Mpc physical radius, not corrected for galaxy background neighbours_halfmpc Expectation number of observed galaxies in cluster within a 0.5-Mpc physical radius, not corrected for galaxy background neighbours_tenthmpc Expectation number of observed galaxies in cluster within a 0.1-Mpc physical radius, not corrected for galaxy background richness_onempc Expectation number of observed galaxies in cluster within a 1-Mpc physical radius, corrected for galaxy background richness_halfmpc Expectation number of observed galaxies in cluster within 0.5-Mpc physical radius, corrected for galaxy background richness_tenthmpc Expectation number of observed galaxies in cluster within 0.1-Mpc physical radius, corrected for galaxy background flag_foreground 1 indicates the presence of a major foreground galaxy (cluster parameters may be unreliable in such cases).Value 0 should indicate good clusters.edge_mask 0 indicates that the cluster is located close enough to the edge of the footprint that parameters may be affected.
Value 1 should indicate good clusters.flag_footprint 1 indicates a galaxy cluster which is located in an isolated area of the footprint (these are excluded from the main dataset).

Figure 1 .
Figure 1.Search cylinder for CluMPR algorithm, oriented lengthwise along redshift (line-of-sight) direction and centered on the potential cluster center galaxy.

Figure 2 .
Figure 2. Search cylinder radius in angular coordinates (radians) as a function of redshift.Obtained from Equation 1.

Figure 3 .
Figure3.Interpolated redshift-dependent binned average of the photo- error as a function of redshift.This is the half-length of the search cylinder.Note that   is limited to 0.2.South refers to the DECaLS survey data, while north refers to the data from BASS and MzLS.

Figure 4 .
Figure 4. Interpolated redshift-dependent thresholds for 1-Mpc and 0.5-Mpc cylinders.South refers to the DECaLS survey data, while north refers to the data from BASS and MzLS.

Figure 5 .
Figure 5. Examples of galaxy clusters found by the CluMPR algorithm in the DESI Legacy Survey imaging.The redshifts of these clusters range from  = 0.305 (left column, third row) to  = 0.707 (right column, second row).

Figure 6 .
Figure 6.Density of galaxy clusters found by the CluMPR algorithm in the DESI Legacy Survey on the sky.

Figure 7 .
Figure 7. 1-Mpc total stellar mass vs photo- for galaxy clusters found by the CluMPR algorithm in the DESI Legacy Survey imaging.Cluster stellar masses have been background-subtracted and corrected for incompleteness with Stellar Mass Correction Factor.The top histogram shows the redshift distribution of CluMPR clusters with the redshift distributions of Zou et al. (2021) (black) and redMaPPer (red) overplotted.The right side histogram shows the stellar mass distribution of CluMPR clusters.
of CluMPR clusters detected byZou et al.

Figure 8 .
Figure 8. 1-Mpc total stellar mass vs photo- for galaxy clusters found by the CluMPR algorithm in the DESI Legacy Survey imaging (same as Figure 7).Color represents the fraction of CluMPR clusters found by Zou et al. (2021) The catalogue by Zou et al. (2021) has broad agreement on high mass clusters.
Zou et al. (2021); however, this percentage increases to 87.1% if one limits theZou et al. (2021) sample to clusters with richness > 30. Figure 10 shows 1 Mpc richness for Zou et al. (2021) galaxy clusters as a function of redshift; color indicates the fraction of CluMPR clusters found by Zou et al. (2021) As in Figure 8, there is broad agreement on high-richness objects.At all redshifts, CluMPR does not find most of the low-richness galaxy clusters identified by Zou et al. (2021) Thus, it is likely that the CluMPR sample is better optimised for detecting high-mass clusters with high purity.An implementation of CluMPR with lower richness thresholds would likely yield results more similar (for  < 0.75) to Zou et al. (2021).
of CluMPR clusters detected byZou et al.

Figure 9 .
Figure 9. 1-Mpc richness (expected number of neighbours) vs photo- for galaxy clusters found by the CluMPR algorithm in the DESI Legacy Survey imaging.Color represents the fraction of CluMPR clusters found by Zou et al. (2021).The lower bound on richness is set by the thresholds in Section 3.1.3(compare Figure 4).
et al. clusters detected by CluMPR

Figure 10 .
Figure 10.1-Mpc richness (expected number of neighbours) vs photo- for galaxy clusters found by Zou et al. (2021) in the DESI Legacy Survey imaging.Color represents the fraction of CluMPR clusters found by Zou et al. (2021).

Figure 11
Figure 11.1-Mpc cluster stellar mass vs redMaPPer richness.The redMaPPer catalogue does not have clusters below richness = 20.The cluster stellar mass has been background-subtracted and corrected for incompleteness with Stellar Mass Correction Factor.There is a clear correlation between CluMPR cluster stellar mass and redMaPPer richness, which indicates that CluMPR stellar masses are likely correlated with total cluster mass as well.

Figure 12 .
Figure12.Residuals between SDSS spectroscopic redshifts and various estimates of cluster redshift from photo-'s.Plotted intervals contain 68.2% of clusters.Our catalogue offers several redshift estimates; the photo- of the cluster central galaxy is least biased (bottom right).The plotted weighted average (middle right) uses mass*probability weighting, which reduces bias compared to an unweighted average redshift (top right).We include a catalogue of the member galaxies for all clusters, which can be used to calculate alternative estimates of the cluster redshift.All plotted redshift estimates underestimate spectroscopic redshift near the high- limit due to sharply decreasing galaxy counts at increasing redshifts (attenuation bias, seeFreeman et al. 2009), which leads to a redshift estimate for the cluster central galaxy which is biased low.
Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.The Photometric Redshifts for the Legacy Surveys (PRLS) catalogue used in this paper was produced thanks to funding from the U.S. Department of Energy Office of Science, Office of High Energy Physics via grant DE-SC0007914.The Siena Galaxy Atlas was made possible by funding support from the U.S. Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0020086 and from the National Science Foundation under grant AST-1616414.Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah.The SDSS web site is www.sdss4.org.SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für As-tronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Figure A3 .Figure A4 .
Figure A3.Median background galaxy counts within 1 Mpc vs photometric redshift.The differences in backgrounds between North and South are due to differing survey characteristics and different models for redshift for each survey.

Figure B2 .
Figure B2.An example of fitting a Schechter function to the observed galaxy stellar mass histogram for the redshift bin centered on z = 0.66.The histogram follows a Schechter function above log(mass) = 11, but below this mass the number counts decrease due to incompleteness.The Stellar Mass Correction Factor (EquationB4) is the ratio between two integrals along the green line (Schechter function): the total mass is the integral above log(mass) = 10, and the incomplete mass is the integral above the peak of the histogram (in this case, above log(mass) = 11).This ratio is designed to compensate for the missing galaxy counts for galaxies with masses below log(mass) = 11.

Figure B3 .Figure B4 .
Figure B3.The histogram peak-finding function used to determine the log(mass) above which stellar mass histograms should roughly follow a Schechter function.The red line is the location of histogram peaks for one example sweep; there is good correspondence between the model and the example.In practice, we use many sweeps simultaneously when training the Stellar Mass Correction Factor (Equation B4), which reduces the random variation in histogram peak location.

Figure C1 .
Figure C1.Candidate gravitationally lensed quasar system near galaxy cluster CluMPR J1278885160133433856 (cluster center galaxy in green).QSO 1 and QSO 2 are confirmed quasars with SDSS redshifts of 1.306 and 1.302, respectively.Objects C1 and C2 are QSO targets for DESI.C3 is a very faint object which resembles a lensing arc.

Figure C2 .Figure C3 .
Figure C2.SDSS spectra for QSO 1 and QSO 2 in the candidate gravitationally lensed quasar system near galaxy cluster CluMPR J1278885160133433856.The broad emission lines (CIII] and MgII) are redshifted to very similar values near  = 1.3, but the slopes of the spectra are very different.This indicates that these are either two different quasars located at the same redshift, or two images of the same highly variable (changinglook) quasar.

Figure C5 .
Figure C5.Candidate gravitationally lensed quasar system near galaxy cluster CluMPR J2439482900096790528 (cluster center galaxy in green).Objects in white circles are QSO targets for DESI, and there is a possibility that they are all images of the same gravitationally lensed quasar.

Table B1 .
Schechter parameters used to derive Stellar Mass Correction Factor

Table D1 .
Main CluMPR Galaxy Cluster catalogue

Table D2 .
Extended CluMPR Galaxy Cluster catalogue of all cluster galaxies, weighted by membership probability z_average_mass_prob Mean photometric redshift of all cluster galaxies, weighted by membership probability times mass z_std_central Std.deviation of photometric redshift of central cluster galaxy z_std_no_wt Std.deviation of photometric redshift of all cluster galaxies z_std_prob Std.deviation of photometric redshift of all cluster galaxies, weighted by membership probability z_std_mass_prob Std.deviation of photometric redshift of all cluster galaxies, weighted by membership probability times mass RELEASE Integer denoting the camera and filter set used for the central galaxy, which will be unique for a given processing run of the data, from Tractor catalogues BRICKID A unique Brick ID for the central galaxy, from Tractor catalogues OBJID catalogue object number within this brick for the central galaxy, from Tractor catalogues MASKBITS Bitwise mask indicating that the central galaxy touches a pixel in the coadd/*/*/*maskbits* maps, from Tractor catalogues gid Unique identifier of the central cluster galaxy.See Section 3.2 for definition.mass_central Log(stellar mass) of the central cluster galaxy cluster_mass_onempc

Table D3 .
CluMPR Galaxy Cluster Membership catalogue Columns Description Galaxy Unique galaxy id, constructed as gid above.See Section 3.2 for definition Cluster Unique galaxy id of cluster central galaxy Galaxy Mass Stellar mass of galaxy Galaxy Redshift Photometric redshift of galaxy Cluster Redshift Median photometric redshift of cluster central galaxy Galaxy Redshift Uncertainty Uncertainty in the photometric redshift of this galaxy Cluster Membership Probability Cluster membership probability of this galaxy

Table D4 .
Candidate Lensed Quasar catalogues Columns Description lens_RA_central Right ascension of central cluster galaxy (deg) lens_DEC_central Declination of central cluster galaxy (deg) lens_z_spec Spectroscopic redshift of central cluster galaxy (-1 if unavailable) lens_z_spec_err Estimated error in spectroscopic redshift of central cluster galaxy (-1 if unavailable) lens_z_median_central Median photometric redshift of central cluster galaxy cluster_mass_onempc Log(stellar mass) of the cluster within a 1-Mpc physical radius, corrected for incompleteness and galaxy background lens_gid Unique identifier of the central cluster galaxy.See Section 3.2 for definition.