LoVoCCS. I. Survey Introduction, Data Processing Pipeline, and Early Science Results

We present the Local Volume Complete Cluster Survey (LoVoCCS; we pronounce it as “low-vox” or “law-vox,” with stress on the second syllable), an NSF’s National Optical-Infrared Astronomy Research Laboratory survey program that uses the Dark Energy Camera to map the dark matter distribution and galaxy population in 107 nearby (0.03 < z < 0.12) X-ray luminous ([0.1–2.4 keV] L X500 > 1044 erg s−1) galaxy clusters that are not obscured by the Milky Way. The survey will reach Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) Year 1–2 depth (for galaxies r = 24.5, i = 24.0, signal-to-noise ratio (S/N) > 20; u = 24.7, g = 25.3, z = 23.8, S/N > 10) and conclude in ∼2023 (coincident with the beginning of LSST science operations), and will serve as a zeroth-year template for LSST transient studies. We process the data using the LSST Science Pipelines that include state-of-the-art algorithms and analyze the results using our own pipelines, and therefore the catalogs and analysis tools will be compatible with the LSST. We demonstrate the use and performance of our pipeline using three X-ray luminous and observation-time complete LoVoCCS clusters: A3911, A3921, and A85. A3911 and A3921 have not been well studied previously by weak lensing, and we obtain similar lensing analysis results for A85 to previous studies. (We mainly use A3911 to show our pipeline and give more examples in the Appendix.)


Introduction
According to the standard ΛCDM cosmological model, galaxy clusters are the most massive objects and the largest gravitationally bound/collapsed structures that assemble in the late universe. The mass build-up of clusters is a consequence of continual merging and gravitational aggregation of matter clumps in the characteristic bottom-up hierarchical scenario per the ΛCDM cosmogony.
Galaxy clusters form in the high-density knots of the large-scale cosmic web, and their abundance at different evolutionary stages traces the high-mass end of the halo mass function. The mass function describes the mass number density in the universe as a function of mass and redshift, and it has the highest sensitivity to cosmological parameters Ω m and σ 8 (Allen et al. 2011) at the high-mass end (LSST Science Collaboration et al. 2009).
Galaxy clusters are comprised primarily of dark matter and hot gas. The gas constitutes only ∼15% of the cluster mass for massive relaxed systems (Mantz et al. 2014); a fraction that is consistent with the overall cosmological baryon fraction.
The total mass of a cluster is not a direct observable, and therefore it has to be inferred from the lensing effects produced (Bartelmann & Schneider 2001) or from other observables such as optical richness (Rykoff et al. 2014); the Sunyaev-Zel'dovich effect (SZ;Sunyaev & Zeldovich 1972;Carlstrom et al. 2002); and X-ray luminosity (Sarazin 1986). There are considerable uncertainities in the mass determinations via each of these methods (e.g., Nagai et al. 2007;Stopyra et al. 2021). For instance, the scatter in the lensing mass reconstructions ranges from ∼10%-20% (Becker & Kravtsov 2011;Bahe et al. 2012) driven mainly by projection effects.
Scaling relations (e.g., Mantz et al. 2016 and references therein) connect the observables and the (true/"expected") underlying halo mass, and make it possible to deduce the massabundance and hence obtain constraints on the cosmological parameters from observations. The scaling relations have much smaller scatter than the measured quantities described above because they are the statistical result obtained from an ensemble of a large number of clusters.
However, to build a robust scaling relation, we require the accurate measurements of both the cluster mass and the observable with careful attention paid to selection biases. Compared with other methods (e.g., member galaxy velocity dispersions (Evrard et al. 2008) and X-ray emission), gravitational lensing provides a direct measurement of the mass of a cluster without any additional assumptions about its intrinsic dynamics or hydrostatic equilibrium-this is especially helpful for merging and unrelaxed clusters. Lensing can then be used to calibrate the scaling relation since the projection effects can be averaged out with ensemble statistics.
On smaller spatial scales, comparing the mass distribution inside a cluster with different observational features, for instance, the galaxy population/luminosity function, or intracluster light (ICL; Montes & Trujillo 2019), or the gas distribution inferred from X-ray emission, reveals how baryons evolve with dark matter. Clusters also permit more fundamental tests for the nature of dark matter. Mapping the spatial distributions of these various components and comparing them allows us to also set the limits on the self-interaction cross section for dark matter at galaxy-group scales, which bridges the gap between cluster-scale (Clowe et al. 2004;Randall et al. 2008;Kim et al. 2017;Golovich et al. 2019;Robertson et al. 2019) and galaxy-scale (Relatores et al. 2019;Sameie et al. 2021) studies.
Previous studies, based on their research goals, often distinguished and selected cluster samples on the basis of their dynamical states. For example, the relationship between gas and dark matter is typically probed in the relaxed clusters under the assumption of hydrostatic equilibrium (Mantz et al. 2016), and in merging systems, to study their interactions (Harvey et al. 2015).
However, in order to fully understand the cluster assembly and mitigate selection biases, a well-defined complete sample of representative clusters is needed for study. In particular, given that clusters have increasingly gained currency as cosmological probes with the derivation of calibrated massobservable scaling relations across dynamical states, we need a standardized framework that allows us to analyze the full diversity of relaxed and perturbed clusters.
In the universe, low-redshift (low-z) galaxy clusters have undergone a long period of growth and evolution; some appear to have reached dynamical equilibrium and are detected as virialized/relaxed structures. Low-z clusters have high spatialresolution observations due to their proximity, and are naturally useful to calibrate the scaling relations in the local universe using their lensing masses. As noted by Stopyra et al. (2021), weak lensing (WL) offers the best prospects as the most robust mass estimator, provided that the accurate redshifts to a large number of local galaxies can be obtained.
Principal uncertainties in the weak-lensing mass measurements include the following: uncertainties in photometric redshifts (photoz) that are used to tag background galaxies; cluster triaxiality and projection effects, and shear calibration that calibrates the measured lensing reduced shear derived from the shapes of gravitationally distorted source galaxies to true reduced shear (McClintock et al. 2019). Detection of multiple strong-lensing image systems with the measured spectroscopic redshifts alleviates these uncertainties. In general, the photoz error of a source worsens as its true redshift increases, but it does not greatly affect the mass measurements of low-z clusters since the lensing distance ratio D ls /D s is almost constant for middle to high-z background source galaxies. Similarly, the projection bias due to the presence of uncorrelated line-of-sight structures can be addressed with detailed modeling of contaminating background clusters with the help of photometry and color information (e.g., Rykoff et al. 2014); also the bias caused by halo triaxiality can be estimated from the shape of the brightest cluster galaxy (BCG; Herbonnet et al. 2019). The projection biases can be reduced by bin-averaging when building scaling relations (von der Linden et al. 2014a;McClintock et al. 2019) as well. Finally, the (reduced) shear measurement biases of weak (and medium) lensing signals can be well calibrated using the cluster-lensing image simulations (Massey et al. 2007;Mandelbaum et al. 2015;Liu 2020), or direct shear response algorithms such as metacalibration (Sheldon & Huff 2017;Sheldon et al. 2020) instead. Low-z clusters also allow for better mass modeling of dark matter substructures, whose persistence is predicted by ΛCDM (Gao et al. 2012), because of the high spatial resolution of observational data.
At smaller length scales, low-z cluster samples significantly improve the studies of the comparisons between the lensing mass distributions and baryonic maps (as mentioned earlier) by taking advantage of high-resolution light profiles, spatial distribution of member galaxies that can be derived observationally to infer the presence of substructures, and features in the gas distribution. Since low-z cluster galaxies are visually larger and brighter, they also allow for higher-quality imaging and spectroscopic studies of galaxy morphologies as a population (Pranger et al. 2013;van der Burg et al. 2016;Sohn et al. 2019); star formation patterns (e.g., in jellyfish galaxies; Moretti et al. 2018). In addition, nearby massive clusters are generally complete in multiwavelength all-sky surveys (e.g., ROSAT; Voges et al. 1999;and Planck SZ;Planck Collaboration et al. 2015), which mitigates selection biases.
A few recent ground-based optical surveys have studied cluster samples including Weighing the Giants (WtG) (von der Linden et al. 2014b), CCCP (Hoekstra et al. 2012), MENeaCS , LoCCuS (Smith et al. 2016), RCS-2 (Gilbank et al. 2011;Hildebrandt et al. 2016), MACS (Ebeling et al. 2001), etc. However, they either have mainly focused on intermediate-redshift clusters, lack broad wavelength coverage, do not have a wide enough field of view (FOV), or lack deep imaging. Many of these samples also comprise preferentially high X-ray luminosity clusters. The future wide and deep LSST (Ivezić et al. 2019) survey will overcome these issues, but there are still several years before it will be completed. Moreover, many of these previous surveys are plagued by selection biases, and therefore a complete cluster survey like the one that we present here is urgently needed.
We note that MENeaCS studied X-ray-selected massive clusters that are in a similar redshift range under similar observation conditions compared to LoVoCCS. As we show below, however, LoVoCCS (1) focuses on the opposite hemisphere and doubles the sample size by extending the X-ray limit, (2) doubles the wavelength range for better galaxy type and redshift determination, (3) is ∼1 mag deeper in r band with doubled FOV, and (4) uses different instruments and pipelines. In addition to mass measurements, LoVoCCS also studies the mass distribution inside the clusters to compare with baryons and to locate substructures.  2016). The typical seeing is ∼0″.9 in r, i, z, Y bands, which is sufficiently small for lensing studies, and ∼1″.1 in other bands, which corresponds to a resolution of ∼2 kpc at z ∼0.1 that can easily resolve low-z galaxies. DECam has 62 high-resolution imaging/science CCDs (2k×4k pixels) with a pixel scale of ∼0″.263 per pixel and well-samples its point-spread function (PSF). It has the widest FOV (2°.2, 3 deg 2 ) among similar currently available systems and is mounted on a telescope with a large aperture, which enables fast, consistent, and comprehensive observation of both cluster member galaxies and background source galaxies across the whole virial region of a typical low-z cluster, and leads to a considerable lensing signalto-noise ratio (S/N) because the low lensing signal due to the factor D ls D l in the critical surface density is compensated by the huge number of resolved background galaxies within the FOV. It also contains u band, which is important for photoz determination (in particular for low-z galaxies due to the 4000 Å and Balmer break); and the u band also improves the photozs of high-z galaxies by detecting the Lyman break (Schmidt & Thorman 2013;Sawicki et al. 2019). The range of observational filters is particularly useful for robust cluster member selection, star formation history analysis, and metallicity-sensitive stellar population studies.
In order to compare and combine the analysis results of different galaxy clusters in our sample, we require a pipeline that can consistently and uniformly process the observational data. In general, each telescope system has its own processing pipeline, which prevents direct one-to-one comparison between the processing steps and algorithm implementation details in different pipelines. The LSST Science Pipelines software (LSP; Jurić et al. 2017) 24,25 is being developed to use the state-of-theart algorithms to process and measure the future LSST images, but is also applicable to other telescope systems (e.g., DECam, Canada-France-Hawaii Telescope, hereafter CFHT; Mega-Cam; Boulade et al. 2003), Hyper-Suprime Cam (HSC; Miyazaki et al. 2012Miyazaki et al. , 2018, Sloan Digital Sky Survey (SDSS; York 2000) once the user provides the characteristics of the system at the beginning of the processing. After running similar calibrations on input images, the software performs the same processing and analysis procedures on the calibrated images regardless of what telescope it comes from. The performance of the LSP has been well tested and verified on HSC imaging data (Bosch et al. 2018;Aihara et al. 2018Aihara et al. , 2019. The data products of the LSP have the same format and are compatible with the LSST, which allows for direct comparison among the outputs of different optical systems and simplifies future comparative studies. Given the accessible instruments and analysis tools, we have designed the Local Volume Complete Cluster Survey (LoVoCCS; PI: Ian Dell'Antonio)-an on-going National Optical-Infrared Astronomy Research Laboratory (NOIRLab) survey. 26 The survey observes a complete, volume-limited sample of 107 nearby X-ray luminous galaxy clusters (Section 2.1) that are not obscured by the Milky Way (MW) and are observable with DECam. We use the LSP and our independently developed semiautomatic pipelines to process the data. Our key science goals are as follows: 1. Map the dark matter distribution in the surveyed clusters and detect additional associated structures like the filaments and voids between halos via weak gravitational lensing; 2. Measure and quantify the effects of substructure on cluster mass measurements; 3. Determine the properties and evolution of cluster member galaxies as a function of cluster mass and substructure; 4. Provide a calibration of cluster-scaling relations at low redshift for systems where individual structures can easily be resolved and compare the lensing derived distribution of the total mass to X-ray emitting gas/the SZ effect; 5. Provide a first epoch baseline for the study of lensed and cluster transients with the LSST; 6. Establish the low-redshift anchor for the study of clusters and their galaxies with the Nancy Grace Roman Space Telescope.
The LoVoCCS survey was started in Fall 2019 and is expected to be complete (including archiving data) in ∼2023. The complete set of observations for each cluster in the LoVoCCS survey reaches Year 1-2 LSST depth (in ∼2024-2025) but concludes 1 yr before the LSST. The LoVoCCS survey is 1-2 mag deeper than DES, and has more stringent selection criteria in terms of seeing before coaddition of individual exposures. Once finished, LoVoCCS will be the largest uniformly selected and individually analyzed sample of nearby massive galaxy clusters. It will permit the determination of the spatial distribution and properties of dark matter in clusters on the scale of 100 kpc ( 1 ¢ and ∼4 ×10 13 M e ), and allow the studies of the properties and interactions of galaxies that are 4 mag below L å and down to a surface brightness limit of μ = 28 mag arcsec −2 . This depth allows for a proper census of low-surface-brightness (LSB) mergers/tidal features in the clusters and in the backgrounds (Adams et al. 2012); estimates of their galaxy merger rates (Drlica-Wagner et al. 2021); a census of LSB galaxies in clusters (e.g., Impey & Bothun 1997;Roman & Trujillo 2017;van der Burg et al. 2017); and a detailed study of the ICL as well (Montes 2019) in a much larger and complete sample. LoVoCCS will constrain the mass uncertainty in scaling relations to <10% (Herbonnet et al. 2020), comparable to the stacked-cluster result of DES (McClintock et al. 2019) but using individual mass profiles with detailed properties derived from a complete cluster sample.
The depth and large FOV of LoVoCCS provides the opportunity to study low-z lensed transients in member galaxies as well as the intra-cluster environment (Sand et al. 2011) down to faint absolute magnitudes to help establish the synergy of these transient observations with other surveys. In addition to its high resolution, LoVoCCS has the ability to capture the strong-lensing features near cluster galaxies, and will trigger follow-up groundbased spectroscopic measurements to refine the constraints on the granularity of the cluster mass distribution (Gladders et al. 2003;Meneghetti et al. 2020). Because of the wealth of available ancillary X-ray and near-IR data, LoVoCCS also provides various opportunities to study active galaxies, star formation, and feedback processes in nearby clusters. We will analyze the LoVoCCS fields that were observed by the Hubble Space Telescope (HST) with an independent pipeline (e.g., Clowe et al. 2006) to calibrate the shear measurements of the LSP. We also plan to compare the observational results of LoVoCCS with cosmological simulations to verify the performance of the LSP. Comparison of the LoVoCCS sample with simulated analogs will enable us to obtain constraints on the nature of dark matter as well.
The study of LoVoCCS clusters will be published in a series of papers. In this (first) paper, we present the survey with details of the observations and the data processing pipeline. We present the data of three clusters (two of them have not been well studied previously via WL) as examples of the use of our pipeline. In Section 2, we describe the LoVoCCS sample and our observational strategy. Then in Section 3, we present the data processing pipeline, including the LSP steps, our data quality checks, and the post-LSP procedure, with examples of processing outputs. We show data products and results in Section 4. In Section 5, we discuss the sources of uncertainties and some known issues of our pipeline, and finally we summarize the work in Section 6. Throughout this work, we adopt the cosmological parameters used in Child et al. (2018) to be consistent with the simulations, where we obtained the cluster concentration-mass relation: flat ΛCDM H 0 = 71 km s −1 Mpc −1 , Ω m = 0.2648, Ω b = 0.0448. The magnitudes are in the AB system (Oke & Gunn 1983).
The LoVoCCS survey images and catalogs will be made regularly available through yearly releases on the NOIRLab archive. Our analysis framework and tools will aid the development of the cluster-lensing pipeline in the LSST Dark Energy Science Collaboration (LSST Dark Energy Science Collaboration 2012).

LoVoCCS Sample
To build a complete sample of Local Clusters accessible to DECam, we start from the NASA HEASARC Meta-Catalog of X-Ray Detected Clusters of Galaxies (MCXC; Piffaretti et al. 2011), 27 which was constructed on the basis of 7 ROSAT All Sky Survey-based and 5 serendipitous cluster catalogs (1743 clusters in total) and is complete at low redshift and high X-ray luminosity (and thus high mass). We select the X-ray luminous clusters [0.1-2.4 keV] L X500 > 10 44 erg s −1 (M 500  2 × 10 14 M e ), 28 where 500 means the total value within a radius, inside which the average mass overdensity of the cluster is 500 times the critical density of the universe at the cluster redshift. For the redshift range, we set a lower limit at z > 0.03 to make sure that the FOV of DECam covers the whole virial region in the cluster rest frame (∼3-4 Mpc), and we set an upper limit of z < 0.12 where the catalog is nearly complete (Böhringer et al. 2000Schuecker et al. 2001). We select r band for (weak) lensing measurements because of its balance between seeing and sky brightness (in a few cases, i band performs even better), and we also include u, g, z bands (and Y band if sufficient archival data exist) to facilitate the derivation of photometric redshifts. Therefore, we require that the galactic extinction from the MW in the SDSS r band be lower than 0.5 mag at the center of the target. Similarly, in order to reduce the contamination of foreground starlight, we select the cluster fields that contain no more than 30 MW stars that are brighter than 13 in G-magnitude within 15 arcminutes of the cluster center (using Gaia Data Release 2, 2018;Gaia Collaboration et al. 2018. Additionally, we require that our targets can be observed at airmass 1.5, and this leads to a cut of the sample at decl. δ < 20°. 29 The reason is that DECam lacks a compensator to correct for atmospheric dispersion, which distorts PSF and deflects sources, and hence limits the photometry, astrometry, as well as lensing measurements. All of these thresholds result in a final sample of 107 clusters (Figure 1; Appendix A). Moreover, we have started using HSC to observe a similar sample of clusters in the northern sky 30 (PI: Hironao Miyatake). A large portion of the LoVoCCS clusters have been observed by space telescopes in the optical/ IR (e.g., HST/Wide-field Infrared Survey Explorer), highquality X-ray (e.g., Chandra/XMM-Newton), and microwave (e.g., Planck) wavelengths as well.

Observing Strategy
We consider typical colors of the galaxies in our sample: nearby ellipticals (Strateva et al. 2001;Baldry et al. 2004;Chang et al. 2006;Cappellari et al. 2011;Jin et al. 2014) and background spirals (LSST Science Collaboration et al. 2009); and require that at r = 24.5, i = 24.0 the galaxies reach S/N = 20 (for lensing measurements), and reach S/N = 10 at u = 24.7, g = 25.3, z = 23.8 (for photoz and cluster-dwarf-galaxy detection). This is based upon the LSST gold sample as i < 25.3 (S/N  20) in the 10 yr catalog (LSST Science Collaboration et al. 2009). For the LSST Year 1 catalog, the galaxies with the same S/N have magnitudes i  24. Our strategy yields a galaxy number density of ∼20 per arcmin −2 (i < 24 or r < 24.5), which enables the measurement of ∼200k background sources within 1 deg (approximately the virial region) of each cluster and produces the lensing S/N  5, and allows the study of member galaxies ∼9 mags fainter than Lå in z band. As a result, our observations 27 https://heasarc.gsfc.nasa.gov/W3Browse/rosat/mcxc.html 28 We note that this L X500 lower end of the sample is expected to be further adjusted due to (1) the inconsistent measurement of L X500 from different input catalogs of MCXC and (2) the steepness of the luminosity function, and we will seek to undertake supplementary observations of those low X-ray luminosity clusters in the future. We have been granted time to observe a sample of low-SZ-mass clusters to extend the SZ-mass-completeness of LoVoCCS (Prop. ID 2022A-658443). 29 A few clusters are close to the south pole and their minimum airmasses are ∼1.7; we still include them in the LoVoCCS sample to maintain completeness. 30 https://www.naoj.org/Observing/Schedule/s20a.html are able to simultaneously well detect the background sources for a lensing analysis and the cluster galaxies for population studies including faint and LSB ones (∼28.5 mag per arcsec 2 in r band; adequate for the detection of the galaxies in Dragonfly-like surveys; Cohen et al. 2018), and is comparable with the LSST for transient/variable source studies. Using the DECam Exposure Time Calculator 31 and assuming a typical moon day at 3, we find the total exposure time needed at u, g, r, i, z are 5100 s, 2300 s, 4200 s, 3400 s, and 2300 s respectively (in u, r, i, we preferentially consider background sources because of photoz and lensing measurements), and the depth is between the planned LSST one-year and two-year observations (u band is slightly deeper).
Each LoVoCCS cluster requires 5.2 hr of DECam exposure including overheads. However, ∼60% of the clusters were partially observed by DES or other community projects before LoVoCCS. NOIRLab granted LoVoCCS 40 nights, and in order to effectively use the allocated time, we observe in r band (sometimes i band) when airmass <1.4 and seeing 0″.9 for getting the high-quality shape information, and otherwise switch to the other bands to obtain the photometry, as long as airmass 1.7 and seeing 1″.5. Because the bluer bands are more sensitive to moon light, we observe clusters, e.g., in u, g before moon rise and g, r, i, z afterwards. We use dither patterns to fill CCD gaps and clean CCD defects/cosmic rays, and use the short exposures (∼100-200 s) to avoid the saturation of reference stars that will be compared with external catalogs (Section 3.1). So far, ∼60% of the data has been taken (including archival observations), and 49 clusters have been observation-time complete, even though the observation was delayed due to COVID-19 in 2020-2021. We have processed the data for ∼20 clusters and obtained preliminary results (Section 3). We will continue observing the individual clusters in 2021-2022, and will make complementary observations in 2023 after evaluating the depth and completeness of processed clusters.

Data Reduction
We created a pipeline named run_steps that incorporates the LSP and our analysis methods to semiautomatically process the LoVoCCS observational data and consistently generate science results (flow chart: Figure 2). For each cluster, we use the LSP to detrend and calibrate the raw DECam data, select the high-quality CCD images, stack the images, and measure the final coadd images (Section 3.1). Then we calibrate the absolute-magnitude zero-points of the LSP output catalog, measure the photometric redshifts, analyze the lensing signal, and derive a cluster mass (Section 3.2). We choose the LSP for DECam data processing rather than DES Data Management (Morganson et al. 2018) or Community Pipeline (Valdes & Gruendl 2014). Our choice is motivated by the following: (1) they use different flagging/masking (e.g., for cosmic rays and satellite trails) and measurement methods; (2) the LSP uses the more recently developed algorithms (or algorithms that are being developed); (3) we aim to compare our results with and apply our tools to future LSST data products, and therefore we require a similar processing framework, consistent methods, and compatible catalog structure; (4) we seek to compare the DECam results with CFHT and HSC/Subaru whose data can be processed with the LSP under the same framework (obs_decam versus obs_cfht and obs_subaru/ hscPipe 32 ; Bosch et al. 2018), in order to calibrate the processing and establish the synergies; and (5) the LSP has a well-developed parallel processing mechanism (per CCD or patch) that saves CPU time. The run_steps pipeline consists of several steps as summarized above, and we will present the details of each step below. We use the LSP version 19.0.0 and

Raw-Data Processing Using the LSP and Check-Visit
This section explains our procedure for retrieving the DECam data and using the LSP to detrend, calibrate, select, coadd, and measure the data in run_steps. The size of the processing outputs of each cluster is about 1.5 TB; the time cost is about 3-5 days with 20 cores (Intel Broadwell/Skylake/ Cascade), and the memory cost is about 200 GB.
download_raw-We obtain the raw DECam data of our observations and the archival data that cover our targets using either NOIRLab API 34 or NOAO Science Archive (NOIRLab Astro Data Archive). 35 initialize_DATA-The initialize_data stage sets up the data repository that will hold processed science images, and directs the LSP to process our DECam data with obs_decam using the previously ingested calibration files. The calibration files include the following: the (nightly) master calibration files (MasterCals) 36 of biases and dome flats; "defect" images produced by obs_decam; z-and Y-band fringe images 37 ; and external reference catalogs for astrometric and photometric calibration.
The MasterCals are obtained using the same method mentioned in download_raw. We find that they are stable within a few weeks, especially after 2013. Therefore, we select the high-quality MasterCals every 1-2 months (per-pixel differences are mostly 1% in flats and 1 ADU in biases between consecutive selected ones) and use them instead of the MasterCals from individual observation nights in the data processing. This saves disk space and produces consistent calibration, and is similar to the strategy in Morganson et al. (2018).
We use Gaia Data Release 2 (DR2; Gaia Collaboration et al. 2018, 2016) for astrometric calibration because it provides the high-precision star positions (in general uncertainty 2 milliarcseconds; mas) across the entire sky. We select the stars that have standard errors < 25 mas in R.A. and Decl., and G-band S/N > 10. As Gaia uses the idiosyncratic passbands (Jordi et al. 2010) that have strong degeneracy in color-terms if they are converted to DECam filters, we use the Pan-STARRS1 Because there is no all-sky deep survey that provides a homogeneous photometry with Sloanlike filter system, we use these three catalogs to cover the footprint of LoVoCCS. The differences between the surveys that we use and their most recent data releases are small for the calibration.
PS1 reaches 5σ depths in g, r, i, z, y < 23.3, 23.2, 23.1, 22.3, 21.3; it covers δ > −30°and lacks u band. We select the stellar objects using the difference between Kron and PSF magnitudes, the number of detections in each band 2, and S/N > 5. SDSS reaches 5σ depth at u < 22, and its footprint is mostly in the northern sky. We select the clean, detected, observed, 41 and deblended stars with S/N > 5 from the catalog.
We build the reference catalogs in individual cluster fields in order to save disk space and increase computational efficiency. Noting that, in δ < −30°, SkyMapper is the only homogeneous survey with Sloan-like bands but not as deep as PS1 and SDSS, we preferentially choose PS1 and SDSS (or SkyMapper-v if no SDSS data exist) for targets in δ > −30°, and use SkyMapper in δ < −30°. In the future, we will seek to use DES as the reference catalog for the targets in the DES footprint, but this may lose the calibration consistency among the clusters.
ingest_raw-This step ingests the raws into the repository, and also adds the information of airmass, which is required by jointcal (see below) to account for refraction, into the header if it is missing (∼year 2013), by using the zenith distance.
processCcd-In this step, the LSP performs instrumental signal removal (overscan correction, crosstalk correction, bias correction, linearization, masking, interpolation, fringe correction, and flat correction), image characterization (iterations of cosmic-ray cleaning, background subtraction, source detection, PSF measurement by PSFEx from Bertin 2013, photometry, and aperture correction), and calibration (of astrometry and photometry using reference catalogs) on the raw images (Figures 3 and 26). This generates calibrated-exposure images (calexp) and source catalogs (src) of single CCDs.
The Bright-Fatter (B/F) effect is that brighter stars could appear larger on the exposure image because of the accumulated charges in the CCD pixels. The correction to the B/F effect is currently not included, but its percent-level bias on shear measurement (Gruen et al. 2015) is much lower than the uncertainties of individual cluster mass measurements, which are dominated by shape noise. The B/F effect will be significant and addressed at ensemble studies, e.g., calibrating the scaling relation by the whole LoVoCCS sample in the future.
Star flats (using the same stars that are observed at different positions of the detector in different exposures to calibrate the camera response) and illumination correction (a correction to the star flats) 42 are not included as well, but their functionality will be achieved by jointcal (see below).
We assume that the small (stellar) color-terms in the mapping between the bands of DECam and the reference 33 https://pipelines.lsst.io 34 https://github.com/NOAO/nat-nb 35 https://astroarchive.noirlab.edu 36 https://noirlab.edu/science/documents/scidoc1203 37 https://noirlab.edu/science/programs/ctio/instruments/Dark-Energy-Camera/Calibration-Files 38 PS1 DR1 is slightly different from DR2 at coordinate ( <0″.01) and magnitude uncertainty ( <0.02 mag) with mean ∼0. 39 The differences between SDSS DR15 (Aguado et al. 2019, recent) and DR12: coordinate <6 × 10 −7 deg (0″.002), magnitude <0.02, magnitude uncertainty <0.01. They have the same ID and the same number of objects. 40 SkyMapper DR2 is deeper than DR1 by ∼1 mag but only became public in late 2020; DR1 was the one that we had access to when we started to process LoVoCCS data in 2019. Nevertheless, we use DR2 in the absolute-magnitude zero-point calibration after the LSP because of its smaller scattering. 41 Here "detected" means the number of detections of this object (nDetect) > 0 in SDSS DR12; "observed" means the number of observations of this object (nObserve) > 0. The observation is related to the runs, while the detection is related to the catalogs. We use both cuts to ensure that the object does not appear temporarily and is useful for the photometric calibration. 42 https://noirlab.edu/science/documents/scidoc1203 catalogs are zero during the processing and correct them later at the coadd catalog level because (1) they will not affect the coaddition, which is mainly weighted by sky background variance (the mean variance of all pixels over a CCD-scale area), and (2) the published color-terms strongly depend on their observation data sets and thresholds. Ignoring the colorterms leads to a small constant offset on the magnitude zeropoint in the LSP output catalog.
We skip CCD 2/S30, 31/S7, 61/N30, and 62/N31 because of their functional problems 43 and data quality issues. In general, the processing of one CCD will stop if it meets any errors, and the calibrated images or catalogs will not be produced. Because PS1 is slightly deeper than SkyMapper, we find it allows 4% more CCDs to pass this step.
check_visit-This stage selects the high-quality calibrated CCD images and eliminates the bad observations from the coaddition. Some of our observation and archival data may contain unsatisfactory local features. Therefore, before conducting any further calibration or coaddition in the LSP, we make a quality cut outside of the LSP in order to improve the measurement on the coadd images. We pick out stars used for PSF-modeling (flag calib_psf_used) from the calibrated exposures and fit them with 2D Moffat (Moffat 1969) functions to calculate seeing (FWHM), and we determine the shape (ellipticity) of PSF by using the SDSS (weighted and PSFuncorrected) second moments of starlight profiles. We take their median values on each CCD as quality indicators (Figures 4 and 27) and identify images where the stars are too blurry or too elliptical due to weather or instrument issues. Typical cuts would be |e| 1 and seeing 1  ¢¢ in the lensing shape measurement band and  1″.5 in photometry bands. We set the thresholds at r band as |e| < 0.13; u, g, i, z, Y bands as | e| < 0.33; r band as FWHM < 4.4 pixel (or 1″.16); and u, g, i, z, Y bands as FWHM < 6.6 pixel (or 1″.74); to allow a small amount of fluctuation because the (weighted) coaddition could reduce the ellipticity and seeing. We also select CCDs whose centers are within 1°.5 ( the virial region) from the cluster center to ensure the coadd depth. Moreover, we create mosaic images of the whole FOV for each visit (exposure) to simplify identifying contamination such as optical defects and satellite trails.
skymap-We feed the calibrated CCD images that passed check_visit into the sky map construction process of the LSP. The LSP distributes the CCD images onto a grid (12 × 12 because of our radial distance cut) of patches with tangent plane projection (0″.263 per pixel)-4200 × 4200 pixels with a 100-pixel overlapping region between neighbors. Each patch is regarded as a flat area, and the difference between the World Coordinate Systems (WCSs) of the patches is a linear translation.
jointcal-This LSP step works on multiple CCD catalogs from the different visits of the same band. It minimizes the total differences between the measurement results and true values (external or to be fitted) to solve for transformations, because the same (stellar and isolated) source should have consistent positions and fluxes-this improves the astrometric and photometric calibrations on the overlapping sky region. It uses the reference catalog (usually shallower) to carry out the absolute calibrations, and we still use the zero color-terms as in processCcd and correct them later in coadd catalogs. When comparing our results with DES (g, r, i, z, Y) and SkyMapper (v), we find that in general setting the orders of the photometric fit photometryVisitOrder to smaller values, e.g., 2, 2, 0, 2, 2, 0 rather than the default 7 on u, g, r, i, z, Y bands, gives smaller scatters. We also filter out CCDs that have a  significantly small number of stars, or visits that contains very few CCDs, because both cases can weaken the statistics and lead to fitting issues.
makeCoaddTempExp-The LSP warps/resamples/reprojects/remaps the processed images onto the patches to make temporary exposures for coaddition. We let the LSP apply the jointcal results to the calexp images, and convolve the temporary exposures with kernels to match their PSFs in order to find the satellite or asteroid trails.
assembleCoadd-The LSP identifies and flags the transient artifacts by comparing the convolved temporary exposures from different visits, and then coadds/stacks the unconvolved temporary exposures (weighted by the inverse of their mean variances) and cleans the artifacts.
detectCoaddSources-After coaddition, the LSP detects the footprints (regions that are above the detection threshold) and their peaks on the final coadd images in each band.
mergeCoaddDetections-The LSP combines the footprints and peaks from different bands. They are not necessarily spatially matched due to their physical spectral energy distribution (SED) features and redshifts. However, the peaks from different bands help with reducing the degeneracy in deblending a merged footprint into individual child objects.
deblendCoaddSources-The LSP distributes ("deblends") the light of each merged footprint ("parent" object) into child objects. The version of the LSP that we use runs deblending in individual bands (the algorithm is similar to SDSS; Lupton et al. 2001); a newer deblender that uses multiband color information to break the degeneracy is being developed (SCARLET; Melchior et al. 2018), and we plan to test its performance on the LoVoCCS data in the future.
measureCoaddSources-In each band, the LSP measures the physical quantities, including the position, shape, and flux, of each object in the last step. The effective coadd PSF at one position is generated by combining the PSF models at that point in the coadded calexp images. The output catalogs are identically identified among the bands.
mergeCoaddMeasurements-The LSP finds a "reference" band (best measurement among different bands) for each object. The reference band will be used in forced photometry (forcedPhotCoadd).
forcedPhotCoadd-In each band, the LSP measures the photometry of the objects using their reference-band positions and shapes.
read_catalog_all-Finally, we extract a catalog of the coordinates, magnitudes, shape quantities of all deblended and well-measured 44 coadd objects over the cluster field using the functions provided by the LSP. The magnitudes are converted from forcedPhotCoadd fluxes. The shape information comes from measureCoaddSources. Repeated objects in overlapping areas of adjacent patches are filtered out.

Post-LSP Data Analysis
We conduct the data analysis, which takes a few hours on ∼20 cores to finish, on the LSP processing results.
combine_patch_color-To simplify and aid the visualization of our deep coadd images, we connect the coadd images of single patches and write new WCS into the mosaic image. Then we create a color RGB image by scaling and combining the mosaic images from three different bands. In the scaling process, the pixel values in each band are (1) linearly normalized into [0, 1] by upper and lower limits, and then (2) stretched by an arcsinh function, in order to display the faint LSB features while maintaining the details at bright regions (Lupton et al. 1999). Though the mapping is unique, it changes the colors (different from the method in Lupton et al. 2004) and reduces the (color) saturation. However, it prevents ordinary objects from appearing too red or too blue, and thus it distinguishes the quasars, star-forming regions, and high-z objects; this helps, for instance, the visual identification of background clusters that need to be considered in the foreground cluster-lensing mass fitting. (Figure 5) photometric_correction-In this step, we conduct MW galactic extinction correction and absolute-magnitude zero-point calibration.  first separate stars and galaxies in the coadd catalog by extendedness from the LSP-a value showing the difference between CModel and PSF magnitudes (Bosch et al. 2018, and references therein). We use the extendedness value in r band because it is the band for shape measurement and has the greatest depth. The CModel photometry fits a source using PSF-convolved index n = 1 and n = 4 Sérsic profiles (Sérsic 1963) to model the source by a bulge-disk decomposition, while the PSF photometry computes the inner product of the source image and the PSF model at that position (therefore it only covers the central flux of a resolved extended source).
Next, we spatially match our stellar catalog with the photometric reference catalog mentioned in processCcd and jointcal. We assume the photometry of the reference catalog has been well calibrated. In each band, after the extinction correction, we use the theoretical color-terms built from model stars to transform the dereddened reference-catalog color of each matched star to "true" DECam magnitude, and compare that value with its measured DECam magnitude (e.g., the left panels of Figures 6 and 28). Then we use the median of the differences as a correction of the absolute-magnitude zeropoint of DECam; the potential DECam-saturated stars are cleaned out by magnitude cuts beforehand (g, r, i, z > 16; u, Y > 14 in corresponding reference-catalog bands). More details are given below.
The model stars are selected from F-and G-type mainsequence and subgiant stars (V and IV) in the synthetic stellar SED library (Pickles 1998) and Calspec standard star library (Bohlin et al. 2014), because these stars have low scatter in color-terms due to metallicity and surface gravity (Wolf et al. 2018;Sawicki et al. 2019). We calculate the theoretical AB magnitude of a model star in band α by Equation (1) The theoretical color-terms are then calculated by fitting the differences between the theoretical DECam magnitudes and reference-system magnitudes of those model stars in one band by a color of the reference; we use g − i and a third-order polynomial as in Schlafly et al. (2018). We only use the reference stars within 0 < g − i < 0.7 in the zero-point calibration because of the color range in the models. Star distances and SED amplitudes are canceled out in the calculation.
The magnitude zero-point calibration varies in bands and cluster fields, and is generally 0.15 mag in g, r, i, z, Y and 0.3 mag in u. As a final check, we compare the stellar loci of the calibrated stars and model stars (e.g., the right panels of Figures 6 and 28), and we find that they are consistent except a 0.1 offset in u (with cuts g < 18 and 0.81 < u − g < 1.6 for good photometry and reducing scatter in colors), which is then subtracted. 49 Finally, the corrected LoVoCCS photometry closely matches DES and is close to DECam Legacy Survey (DECaLS; Dey et al. 2019). We present the comparisons in Appendix C.1 (Figures 29, 30, and 31).
photometric_redshift-We use the framework of Bayesian Photometric Redshift (BPZ; Benítez 2000; Benítez et al. 2004;Coe et al. 2006) to estimate the photozs of source galaxies. We refactored the public code 50 to improve its performance (e.g., using Python broadcasting, multiprocessing, and reducing repeated calculations) and to more effectively control the priors and script details. 51 We briefly summarize our implementation here.
Equation (2) shows how well a template SED fits an observed galaxy. f α and f tα (z) are the fluxes of the observed galaxy and the redshifted template SED in band α respectively; a is a scaling factor that absorbs the coefficients and minimizes χ 2 at a ; ¢ f s a is the measurement error. a F F ot tt To be specific, from Equation (1), we have f C 10 m for a redshifted template SED (Blanton & Roweis 2007), C α ∝ ∫T(λ)λ −1 dλ, σ fα ∝ f α dm α , where dm α is the measurement uncertainty. The amplitudes of the fluxes are reduced during the calculation; only the colors (c) matter.
Once the a a 2 ( ) c = ¢ values on a grid of (z, t) are obtained, we calculate a likelihood by P c z t , e 2 2 2 c¢ is the mininum on the grid-this avoids numerical overflow issues. The probability distribution P(z) of BPZ is then calculated by P(z|c, m 0 ) ∝ ∑ t P(z|t, m 0 )P(t|m 0 )P(c|z, t), where m 0 is the magnitude in a reference band.
We use traditional metrics (Benítez et al. 2004 We use the default priors derived from HST because they yield the best performance on our data over the other published ones. We experimented with adding a Gaussian peak in the redshift prior at the cluster redshift, and found that it hardly affected the photozs of bright objects (Figure 7) because they have small magnitude errors, and their likelihood functions are sharp, nor did it greatly affect cluster mass measurements (only percent-level) because we use a large cut on odds. In the future, we will calibrate the priors using all galaxies that have spectroscopic redshift (spec-z or z s ) records in the footprint of LoVoCCS and seek to obtain additional spectroscopy. Similarly, we adopt the default CWWSB templates (8 model SEDs including elliptical, spiral, irregular, and starburst galaxies; Coe et al. 2006 and references therein) after comparison with other public libraries and adopt a two-step interpolation between consecutive templates after tests. We set the maximal redshift at 1.5 based on the observation depth; the correction for the Lyman series and neutral hydrogen absorptions 52 in Madau (1995) is currently not included (it has no effect on the galaxies at z < 1.5 observed by DECam).
After the correction of photometry, we measure the photoz of source galaxies by the methods above. We include the uncertainties from the corrections, and also a noise ∼0.002 to compensate possible underestimated uncertainties and achieve better performance. For galaxies that have spec-z records, we evaluate the performance of our code by these metrics (Ilbert et al. 2006): (1) normalized-median absolute deviation (NMAD) 1.48 × median{|z s − z b |/(1 + z s )} for the scatter; (2) the percentage (η) of catastrophic error |z s − z b |/ (1 + z s ) > 0.15. We find NMAD 0.05 and η < 10% in our data, and they drop as the proportion of background galaxies (to foreground cluster galaxies) increases and are then comparable to the results in Ilbert et al. (2006), which is ∼0.6 mag deeper than LoVoCCS especially at u band, when the foreground galaxies are filtered out.
In Figure 7, we show a comparison between our photoz results and archival high-quality spec-z data 53 (i  22) in several cluster fields with cuts z s < 1.5, odds > 0.95, and 1 mod 2 c < (i  21.5); we have NMAD=0.050 and η = 6.8%. For this selected sample of galaxies, we find the odds values are concentrated near 1, which makes it difficult to select background sources using a probability cut. On the other hand, the galaxies near the diagonal tend to have smaller ; mod 2 c we have NMAD = 0.052 and η = 8.1% when the mod 2 c cut increases to 4. Cluster member galaxies concentrate near the cluster spec-zs; they and some low-z field galaxies can have high photozs, but most of them are bright and can be filtered out by magnitude. The photoz measurement has better performance on background galaxies at higher redshifts for lensing studies.
We visually inspected the outliers (the catastrophic photoz errors) in Figure 7 and found that they were caused by different effects. For those with small spec-z but large photoz (A), some (∼30%) are isolated galaxies, most are blended objects. Those isolated galaxies have special colors (because of attenuation, star formation history, merger, AGN, etc.), or are too far away from the cluster center and not covered by u band in our observation (but covered by archival data in other bands). The blending cases have two types: (1) a redder neighbor, and (2) a bluer neighbor that caused the primary red galaxy to be recognized as a blue galaxy in high redshift. Outliers with small photoz but large spec-z (B) were caused by similar effects. Figure 7. Comparison between the photoz (z p ) and spec-z (z s ) in all three cluster fields. Left: the blue points, orange triangles, and green squares represent the data in the A85, A3921, and A3911 fields respectively. The dotted, dashed, and solid lines denote 0.15(1 + z s ), 0.05(1 + z s ), and 0 biases respectively. The outliers at upper left (A) have a large fraction that do not have u-band coverage (gray markers), which is much less frequent in our final catalog for lensing analysis. The outliers at bottom right (B) can mostly be filtered out by a photoz cut. Right: the solid histogram shows the distribution of (z s − z p )/(1 + z s ). The orange dotted Gaussian curve ( 0, NMAD ( ) ~) is for reference only and is normalized to the galaxy count; in theory it would match the histogram if the photozs have small biases. Both A3911 and A3921 fields were observed in u, g, r, i, z, Y; the Abell 85 (A85) field was observed in u, g, r, i, z, and processed using PS1+SDSS as reference catalogs. In the A3911 field, the cluster member galaxies (z s ∼ 0.1) dominate the sample of the galaxies that have spec-zs. Cluster member galaxies and some low-z galaxies (z s  0.2) cause the distribution of the sample to shift to higher photozs, and as a consequence, there is a discrepancy between the histogram and the Gaussian reference.
In Figure 8, we consider a clean sample of galaxies that have S/N > 10 in all five bands, and we make cuts at r > 17, z b > 0.15, and r-band blendedness < 0.42 to select out objects that are likely in the background. We plot the number density of the objects along cluster radial distances, and we stack the results from all three clusters for better statistics. The number density of the "background" ones (the green curve with triangles) shows no clear sign of cluster-galaxy contamination -it is almost flat in 1 < R < 3 Mpc, and it decreases near the cluster centers due to the obscuration and blending of cluster galaxies, magnification (Medezinski et al. 2010;Herbonnet et al. 2020), etc. The drop in R > 3 Mpc is mainly caused by the relatively small FOV of DECam compared with the low redshift of A85. Adding cuts on BPZ-quality metrics, galaxy shapes, and S/N do not greatly change the shape (but amplitude) of the curve for the background objects. In the next step, we use the selections on magnitude, BPZ, and blendedness in addition to other cuts to obtain a sample of sources for the lensing analysis; we use a lower cut of S/N to improve the statistics. A small fraction of cluster galaxies could still exist in that sample. They are not expected to affect the mass S/N map peaks because their shapes are averaged out, but they could bias the measured-cluster mass low. We estimate the effect of the remaining cluster galaxies on the mass measurement in Section 5.
mass_map-In this step, we analyze lensing signal by establishing aperture mass maps. We use the HSM (Hirata & Seljak 2003;Mandelbaum et al. 2005) algorithm integrated in the LSP to measure galaxy shapes. We adopt the following quality cuts for galaxies in our shear catalog; these are modeled in part on thresholds used to make the HSC shear catalog (Mandelbaum et al. 2018a).
(1) HSM PSF-corrected distortion |e| <4: this filters out measured shapes that are too large and allows some |e| > 1 shapes caused by noise.
(2) HSM resolution > 0.3: the resolution compares the size of a source and PSF; an unresolved source has resolution 0.
(4) r-band blendedness <0.42: the blendedness from the LSP estimates the fraction of parent-image light in the neighbor of a child image-an isolated source has blendedness = 0; this avoids deblender errors. (5) r-band CModel magnitude between 17 and 26, and S/N > 5: the lower end of the magnitude cut is close to the saturation limit, and is used to remove low-z field galaxies and the cluster member galaxies; the higher end is based on the LoVoCCS depth. (6) BPZ 0.15 < z b < 1.4: this avoids the contamination by cluster member galaxies and eliminates those high-z galaxies (only occupy a small percentage) that would otherwise require correction for hydrogen absorption features. (7) BPZ odds >0.95: the probability distribution peak is highly concentrated. (8) BPZ 4 mod 2 c < : the observed color is well fitted by a template.
(9) We only select the central 6 × 6 patches ( ∼ 1°.8×1°.8) that contain the high-quality deep-coadd results. (10) We use the shear calibration script 54 for the HSM algorithm to convert the galaxy shapes to per-object galaxy reduced shear estimate g 1,2 (Mandelbaum et al. 2018a) and weights, and similarly, we add a cut g 2 |ˆ| < .
The shear calibration parameters are derived from the HSC image simulations (Mandelbaum et al. 2018b). The LoVoCCS survey has similar instruments, depth, and observing conditions to the HSC survey. Ideally, the shear calibration parameters should be estimated by the image simulations and measurement pipeline of the same survey, but this is beyond the scope of this work, and the correction is much smaller than the mass measurement uncertainty of individual clusters, where galaxy shape noise dominates. We are studying recalibrating the measured shapes to true shear by processing simulated DECam images using the LSP; a new calibration based on source galaxy properties (e.g., S/N and size) will be presented in the future when we analyze the whole LoVoCCS sample.
After the photometric redshift measurement, we use aperture mass statistics (Schneider 1996) to locate lensing mass peaks. The aperture mass is a weighted integral (convolution) of the convergence κ within a finite circular region (aperture/ window) and can be used to detect mass clumps in a κ map. However, κ is not observable and is small in WL. Instead, based on the relationship between κ and shear γ 1,2 , the aperture mass can be rewritten as an integral of tangential shear weighted by another filter function Q. The S/N of the aperture mass is maximized when the filter Q has the same angular scale as the mass peak's tangential shear.
The mass_map stage employs the Schirmer filter for aperture mass maps (Schirmer 2004;Schirmer et al. 2004;Hetterscheidt et al. 2005). This filter is defined in Equation (3); x is dimensionless, R is a plane distance, R ap is the aperture radius, and the selection window Q(x) peaks at x ∼ 0.1. The Schirmer filter gives small weights at both very small and large radii to reduce noise and fluctuation. It is useful for isolating Navarro-Frenk-White (NFW)-like substructure at length Figure 8. Number density of objects as a function of radial distance R. The blue curve with circles represents the clean sample of galaxies (stacked and binned from the three clusters and S/N > 10). The green curve with triangle markers shows background galaxies selected by r magnitude, photoz, and blendedness. The orange curve with square markers gives the opposite of the selection (i.e., the difference between the blue and green curves). 54 https://github.com/PrincetonUniversity/hsc-y1-shear-calib. In fact, we get close results for the mass maps and lensing masses if we only consider a fixedrepresentative shear responsivity The S/N map can be obtained using Equation (4) (Fong et al. 2019), where R and R¢ are plane coordinates (pixels on the image or angular positions) under small angle approximation, ò describes the galaxy shape, 55 and j = T, X (tangential and cross components) correspond to E-and B-mode values respectively.
Ideally, lensing signals vanish in the B-mode, and thus in practice the B-mode serves as an estimation of the noise. Before the calculation, we bin the data using the weighted means in 100 × 100 pixels to reduce noise. Since the calculation at each element is independent, we speed up the process by parallel computation. Our method recovers the result of the significance map approach that uses catalog shuffling in McCleary et al. (2020). The maps give the lensing S/N at each binned "superpixel" rather than the detection S/N of the cluster; a better estimator for the cluster-lensing S/N will be given in the step of global_mass_fitting.
We then record the mass peak coordinates, which will be fed into global_mass_fitting, when R ap maximizes the S/N of the central cluster (within 0°.5) at some cluster-scale size. WCS is also added into the map. Figures 9, 10, and 34 show examples of the S/N maps where R ap maximizes the E-mode S/N. These S/N-maximized maps will be used in the next step to be combined with the optical images and X-ray maps.
combine_map_optical-In this step, we cover the optical image generated by combine_patch_color (or in a single band) with the (S/N-maximized) aperture mass S/N map and X-ray map contours to display their correlations (Figures 11, 12, and 35, where the X-ray data all come from Chandra).
We notice significant offsets between the mass S/N map, Chandra X-ray map, and optical BCG peaks in a number of LoVoCCS cluster fields, and some of the discrepancies have not been well studied. We will give more details in Section 4.3 and follow-up papers.
global_mass_fitting-In this step, we jointly fit the shapes of all selected background source galaxies with the lensing distortion produced by a NFW profile (Navarro et al. 1997;Equation (5)). In the equation, ρ crit = 3H 2 (z)/(8πG) is the critical density of the universe at redshift z, r s is the scale The lensing critical density Σ crit (Equation (6)), convergence κ (Equation (7)), and shear γ (Equation (8)) of a NFW mass halo are given in Wright & Brainerd (2000) (f 1,2 are monotonically decreasing functions of x(>0); x is dimensionless). D s , D l , and D ls are the angular diameter distances between the observer and the source, the observer and the lens, and the lens and the source respectively. v c is the speed of light.
The redshift information comes from NED (for the cluster lens) and the photoz measurement (for the background sources). We take the mass peak from mass_map as the lensing center, and we adopt the concentration-mass (c-M) relation (Equation (9)) in Child et al. (2018), which was derived by fitting both the relaxed and unrelaxed individual halos in cosmological n-body simulations using NFW profiles.
The best-fit mass is obtained by minimizing the total weighted difference between the per-object galaxy reduced shear estimate and the model-predicted reduced shear. The shape noise is assumed to be random. The coordinate transformation is given by Equation (10), where j is the position angle (counterclockwise) of the source toward the lensing center.
e e e e e e cos 2 sin 2 , cos 2 sin 2 10 If multiple high-S/N mass peaks exist (due to, e.g., merging clusters or high-z background clusters), we first search for their optical counterparts; without a concentration of galaxies, it is unlikely that the mass detection is valid. Next, we obtain their redshifts from the NED extragalactic database. In some special cases, NED does not include the redshift records of some neighboring galaxy groups, and we use our photoz results instead, or take the redshifts as variables in the mass fitting. There is also an opposite case: some neighboring galaxy clusters have the redshift records but no significant mass S/N peaks-we still include them in the fitting process, but we find their fitted masses are usually small. In general, except distinct merging clusters, the mass scales of the neighboring clusters/ groups are much smaller than those of the main mass halos in LoVoCCS (since they are selected by high X-ray luminosity), and our tests show that the neighbors have low fitted masses and have small effects on the mass fitting of the main halos (a main halo can have multiple small-scale peaks showing substructure on the mass map at small aperture R ap ). The goal of LoVoCCS is to measure the masses of those main halos, which correspond to massive clusters, rather than the masses of neighboring (e.g., radial distances  r 200 ) less-massive clusters/groups.
We then fit the masses of individual peaks simultaneously by seeking the global minimum of the total difference between the theoretical g 1,2 and the calibrated, measured g 1,2 . The κ and γ 1,2 produced by individual peaks are additive, but g 1,2 are not.
To estimate the uncertainty of the fitted mass of a cluster, we make a histogram of the fitted masses generated by resampling (randomly selecting 50% of source galaxies each time). We build a probability density function (PDF) of the cluster mass using the kernel density estimation. The measured mass is then estimated by the median of the PDF, and its uncertainty is estimated by 34.1% (assuming a Gaussian-like distribution) of the cumulative distribution function on each side of the median and then scaled by a factor of 1 2 (McCleary et al. 2020; Figure 13).
quality_check-We perform consistent tests on each cluster field, e.g., whether the coadd catalog has achieved our aimed depth, to examine the quality of our data, and to determine whether there is a potential problem. The functionality of this step will be described in Section 4.2, and some parts of it have been mentioned in previous steps (e.g., the photometric comparisons with DES and DECaLS).

Data Products and Results
In this section, we present the LoVoCCS data products, data quality, and early science results. As an example, we show the results of cluster A3911 and A3921, which are exposure-time complete in LoVoCCS and have not been well studied by lensing.

Data Products
The data products of each LoVoCCS target include the following: 1. Image products: coadd single-band images for general science studies, color-combined coadd images for visualization, mass maps, and overlays with optical and X-ray images, single-epoch calexp images for transient studies; 2. Catalog products: coadd single-band catalogs, merged multiband coadd catalogs (with derived quantities such as Figure 10. Counterparts in A3921 (z = 0.094). The peak at the lower left can be spurious since it covers only ∼1/4 of the number of galaxies compared with the central peak, and we did not find distinct optical counterparts in that region. Note, at the central peak, R r 0.1 1200 pix 316 0.55 Mpc 2 s ap~~¢ ¢~~, and the discrepancy may come from the halo shape (more details in Section 4.3). photoz), and single-epoch catalogs as counterparts of the image products.
The data products will be published on NOIRLab website one year after the completion of each target.

Statistics and Data Product Quality
We assess the quality of the data products by analyzing the statistics of various representative quantities. We consider a sample of well-measured sources that are <1°from the field center (few Mpc or the length scale of virial region) to include the data with sufficient depth, and 10 > ¢ to reduce the blending effects caused by BCG, cluster member galaxies, and ICL, unless otherwise specified.

Depth
We first verify whether the coaddition has reached the expected target depth of LoVoCCS by comparing the magnitudes and magnitude uncertainties of coadd point sources (PSF magnitude) and extended sources (CModel magnitude) in the sample. In Figures 14 and 15, we find that in both fields both types of the sources have close depths (A3911 is slightly deeper by ∼0.2 mag). CModel gives slightly larger uncertainties because of extended measurement areas and parameter fittings. The magnitude scatter at a fixed uncertainty is caused by CCD gaps and dither pattern. Though the median depth is ∼0.4 mag brighter than the target depth, the maximal depth is very close to or even slightly fainter than the target depth, and the difference between the maximal and target depth is likely caused by filtered visits/CCDs, weather conditions, and the archival data.
Another approach to assess depth is to build a histogram of magnitudes for a small interval around a S/N value (Abbott et al. 2018). We find that A3911 is ∼0.2 mag deeper than A3921. The median PSF magnitude is ∼0.3 mag deeper than CModel and is close to or even slightly fainter than the target depth (within 0.2 mag). The maximal CModel magnitude reaches the target depth (within <0.1 mag). (See Figures 16  and 17.)

Completeness
A property related to depth is completeness. We study the completeness of source detection by the number counts of magnitudes, and determine the limit magnitude of 50% completeness with an interpolated trend of number counts for PSF (or CModel) magnitudes. We consider both PSF and CModel magnitudes of all sources because at the faint end Figure 11. Example of an optical r-band coadd image (gray scale inverted for clarity) of A3911 overlaid with lensing mass S/N map (color coded), and hot gas density (the square root of the X-ray flux observed by Chandra; pink) contours. The distances between their peaks are small. The background oversubtraction near bright and large objects is amplified by arcsinh image scaling for display (Section 5). The objects affected by saturation are filtered out in the final catalog. galaxies tend to be unresolved and can be incorrectly classified as stars. Also, because the cluster fields are selected to be not obscured by the MW, the number density of stars is about 5 times smaller than that of galaxies at the twenty-third magnitude. We find the limits are ∼25.5 in A3911 and ∼25.4 in A3921. Similar to the results of depth, here generally A3911 is ∼0.1 mag deeper than A3921, and the limit in PSF magnitude is ∼0.1 mag deeper than CModel magnitude. The slope of the histogram at the bright end (23 < m < 24) in base-10 logarithmic frame is ∼0.3 in r, i, z bands, which is consistent to the values in Benítez et al. (2004) and LSST Science Collaboration et al. (2009). (See Figures 18 and 19.)

Shape and PSF Measurement
The number density of background source galaxies is essential for the quality of cluster-lensing mass measurements. The original number density of galaxies (extended sources) in the sample is ∼50 per arcmin 2 in both A3911 and A3921. The LSST gold sample cut (i < 25.3; close to our 50% completeness limit because we only reach LSST 1-2 yr depth) leads to ∼40 per arcmin 2 , and the i < 23.8 cut (S/N > 20, full completeness) gives ∼18 per arcmin 2 ; then adding the shape measurement cuts results in ∼14 per arcmin 2 (removes ∼22% sources from the result of the i < 23.8 cut). These are consistent with the LSST Science Book (LSST Science Collaboration et al. 2009). Using the cuts mentioned in mass_map (except the area cut), we find that the density drops to ∼8 per arcmin 2 (the photoz-quality cut alone reduces ∼80% sources of the original sample).
We check the star shape residual after PSF-correction by comparing the second moment shapes of all stars that are used for PSF-modeling and the PSF model at their positions. We test both SDSS shapes (base_SdssShape_psf and base_ SdssShape) and HSM shapes (ext_shapeHSM_HsmPsf-Moments and ext_shapeHSM_HsmSourceMoments). The result shows ∼90% of the stars have shape residual <0.01 and ∼97% are <0.02. (See Figures 20 and 21.) We also compute the correlation between the PSF-corrected shapes of galaxies and the original shapes of stars (PSF) to check the PSF-correction systematic errors (similar to McCleary et al. 2015McCleary et al. , 2020. We use the cuts from mass_map (except the area cut) to select a sample of galaxies. We use the same cuts except the BPZ and HSM cuts, use the cuts on PSF magnitude instead of CModel, and flag calib_ps-f_used==1 to build a sample of stars. We then use the HSM shapes of galaxies and SDSS second moment shapes of stars to calculate the correlations in Equations (11) and (12), where R¢ is the plane position of an object, R is the unit vector of R denoting a direction, and the average is for all pairs separated by a distance R ( R R = ). In principle, the shapes of Figure 12. Counterparts in A3921. The BCG is surrounded by X-ray contours. We will discuss the large discrepancy between the lensing and X-ray /BCG peaks in Section 4.3.
galaxies are irrelevant to PSF after the PSF-correction, and thus star-galaxy shape correlations vanish. In practice, the shape noise can cause fluctuations in the correlations. We find that all correlations are small (∼10 −4 , i.e., percent-level of reduced shear; similar to Mandelbaum et al. 2018a

Astrometry
We examine the LoVoCCS astrometry by comparing its stellar astrometry with Gaia DR2, which use the International Celestial Reference System at epoch J2015.5; most surveys use J2000, and the proper motions will cause coordinate shifts. We find that in both cluster fields ∼97% sources have difference <0″.1, ∼ 80% sources have difference <0″.04 (∼the standard deviation of the Gaia-matched sample) , and∼30% sources have difference <0″.01. The outliers are mainly caused by blending, star-galaxy separation, or shredding in deblending. The asymmetry of the difference distributions possibly results from proper-motion correction, the asymmetry of the observation area on the celestial coordinate grid, and the distance to the Galactic plane (i.e., star density; Abbott et al. 2018). (See Figures 22 and 23.)

Photometry
We test the cluster red-sequence versus the spec-z and photoz of galaxies on the magnitude-color diagram. We use all Figure 13. Histograms of 1000 resamples of global mass fitting for A3911 and A3921 respectively (top: (a), (b)); A85 and its neighboring clusters A87 and A89 (bottom: (c), (d)). The solid curves represent the kernel density estimates of the probability density functions of the measured masses. The fitted mass of A85 is close to the results in Herbonnet et al. (2020) (<1σ) and McCleary et al. (2020) (<3σ), but the instruments and/or analysis pipelines are different. In addition, we find that a further augmented multipeak fit including cluster substructure and a nearby galaxy group makes no significant difference (<1σ) to the total mass of A85. The kernel density estimates for A87 and A89 are skipped because of their small lensing masses; instead, we estimate their upper limits using the median values respectively: 1.09, 0.25 × 10 14 M e . The coordinates and redshifts of A85, A87, and A89 are from the NED database. Using the lensing peak of A85 reported in Table 1 to perform a single-peak fit gives M 10.76 10 2.17 2.60 14  -+ . well-measured galaxies in the catalog to obtain enough cluster member galaxies for showing the red-sequence in this test. In Figures 24 and 25, spec-z shows a group of bright-red cluster galaxies that are concentrated along the red-sequence, and some fainter-blue cluster galaxies are off the red-sequence. Though photoz (BPZ) slightly biases red cluster galaxies toward higher redshifts (0.1 < z b < 0.15), the spread of photoz on the red-sequence is still tight, and thus we set a cut at z b > 0.15 to select background galaxies for cluster-lensing mass measurement; though some cluster galaxies could have higher photozs.
Overall, these test results are consistent with our expectation.

Science Results
Our science results in this work are as follows: (1) We summarize our global mass fitting results of A3911, A3921, and A85, and their lensing peaks in Table 1. The lensing mass S/N of A3911 is ∼4, ∼5 for A85, and ∼2 for A3921. The X-ray mass (Piffaretti et al. 2011) andSZ mass (Planck Collaboration et al. 2016) of A3911 are within 1 error bar of its lensing mass, and within 2 error bars for A3921 and A85. The significant mass differences at A3921 and A85 could come from their dynamical states; details are given below. A more in-depth comparison of the lensing/X-ray/SZ masses and peaks of an ensemble of clusters will be presented in a followup paper.
(2) A3911 shows an agreement between the lensing, BCG, and X-ray peaks with small offsets ( 1.5 160 kpc; ¢F igure 11); the SZ peak is close to the BCG ( 18 ¢¢). While A3921 (Figure 12) shows correspondence between the BCG and X-ray peaks (offset 8 ; ¢¢ and the SZ peak is only 1.2 ¢ away from them), a notable difference between baryon and mass peaks exists ( 6 0.6 ¢M pc). We find no distinct bright large member galaxies near the mass peak. We do notice a concentration of member galaxies (through red-sequence and spec-z) along the NE-SW overdensity region (more details below). For A85, we see the existence of the mass substructure (since it is at about half of the distance from A3911/A3921 to us) near the cluster center when using small Schirmer apertures (similar to McCleary 2017) and mass clumps corresponding to nearby clusters and galaxy groups; though some mass clumps lack optical counterparts or published redshifts. When making the mass S/N maps of A85, we find the S/N of the central peak keeps increasing with the aperture size and maximizes at ∼7 when R ap = 32,000 pix. While the BCG and X-ray peaks are close ( 7 ¢¢), the WL, BCG, and SZ peaks are 3 ¢ (∼190 kpc) away from each other, which is likely due to its ongoing relaxing process (McCleary 2017, and references therein). Additionally, because of the binning, Figure 14. Magnitude vs. uncertainty for individual point sources and extended sources in u, g, r, i, z bands in the selected region of A3911. The blue points are the PSF magnitudes of point sources, and the orange points are the CModel magnitudes of extended sources. The solid, dashed, and dotted lines denote S/N = 5, 10, 20 respectively. The median depth is u ∼ 24.5, g ∼ 25.0, z ∼ 23.3 at S/N = 10 and r ∼ 24.3, i ∼ 23.8 at S/N = 20 (the zero-point calibration has been applied; same below). Here we present the detection limit only, and therefore the uncertainty from the correction is not included. The clumps of points isolated from the curves mainly result from blending, CCD gaps, and shredded sources deblended from bright and large objects. the resolution of the mass map is 26 , and there could be a small bias 1  ¢ due to shape noise (McCleary et al. 2020), or possibly cluster triaxiality.
(3) Previous studies show that A3921 is a merging cluster between SE (main) and NW (secondary) components, and the subcomponent is moving along SW/NE direction rather than "face to face" (Ferrari et al. 2005;Pranger et al. 2013). The XMM-Newton observation with a wider FOV reveals that the spherical and hydrostatic assumptions cannot provide a good fit for the gas distribution on a large scale (Belsole et al. 2005). Therefore, the true mass distribution is likely to have a distorted shape, and the correlation between lensing mass map overdensity and member galaxy distribution suggests that our result is not necessarily unreliable. However, for the goal of accurate cluster-lensing mass measurement, a more complicated model has to be used, because all masses mentioned in (1) could be biased because of their assumptions. We will study this in the future, e.g., using simulations to build more accurate mass profiles.

Error Budget
In this section, we discuss the sources of the uncertainties other than the shape noise. They are not fully accounted for in the mass error bars and are independent of each other. However, we expect that they have small effects on the mass measurement; we give the details below.
1. Astrometry. It has small effect on the mass measurement because we use binning when building the mass map, and the reduced shear is a slow function of the radial distance at weak-lensing regions. Also, the quality check (Figures 22, 23) shows that it has very small differences from Gaia. 2. Photometry and photoz. As shown earlier, our photometry has small scatter and zero-point offsets from DES, but we have no data set to directly compare u-band photometry without color-terms. These have small effects on our magnitude cuts and depth, but can potentially affect photoz. Because photoz mainly depends on colors, the bias is generally random; it also has a very small effect on the BPZ magnitude prior. Our photoz measurement has ∼5% errors on average for sources at i  21, and even smaller errors when member galaxies are excluded. The error can be larger for fainter sources, but fortunately, for low-z clusters, the lensing distance ratio changes slowly for the sources at middle and high redshifts, and thus the photoz errors cause small variance on the measured mass. For instance, we tested and found that a 5% variance on the redshifts of all source galaxies would cause a 4% change on the mass results, and the resampling (i.e., shape noise) could also contribute to that difference.
In order to estimate the cluster member (residual) contamination, we consider (a) using color-color (CC) cuts (Medezinski et al. 2010), or (b) using a higher limit on the photoz, to select a new sample of background sources and rerun the mass fitting. For (a), we find u − i versus g − z gives the best separation between low-and high-redshift galaxies after comparing with galaxy SED models, and the selected ones are almost in a subset of the previous sample (selected by r magnitude and photoz)the unselected ones even include some low-z background galaxies and show lensing signal. In (a) and (b), on average the fitted mass increased by 10% with higher error bars, which could come from the cluster member contamination, if we assume it is less frequent in the selected color space or higher photoz, but could also be caused by the reduced shear sensitivity at lower redshift. In Figure 7, a large fraction of the outliers are caused by blending, but the parameter "blendedness" is not able to filter out all blended objects. Blending causes the deblender to redistribute the photons from neighboring sources, and we note that the outliers do not reside in a special region in the color space. If we assume the rate of the faint objects that have excessive photozs is close to the value in the bright end, then it is ∼15% at low redshift (Figure 7). Given that the ratio between the numbers of the background and foreground/cluster galaxies is ∼3 (Figure 8), we estimate the dilution of measured, reduced shear is ∼15/(3 + 0.15) ∼5%. It is possible that photozs for faint objects are more biased, but for faint objects, the ratio of foreground/cluster galaxies to background galaxies will be smaller, leading to an effect on mass estimates of a few percent, which is close to the result of the CC cuts. These results are similar to those of Medezinski et al. (2018) but using different instruments and pipelines. In a subsequent paper, we will investigate more fully the uncertainties in photoz and the use of CC space cuts using a larger sample of clusters. However, to thoroughly solve this issue, large numbers of image simulations are required (e.g., Herbonnet et al. 2020). The study of the (cluster) blending effect on background source photoz and shear measurements is beyond the scope of this paper, and we expect it is at percent-level on average. 3. Shape measurement. We have shown that the PSFmodeling and PSF-correction have no strong biases, and the error from shear calibration is likely at percent-level. In fact, the shape noise of source galaxies causes much larger errors (∼30%) 56 on the lensing mass measurements of individual clusters. 4. Uncorrelated large-scale structure (LSS), correlated LSS (2-halo term), and halo triaxiality. In our case, the uncorrelated LSS mainly consists of background clusters and large-scale matter distribution (cosmic noise) along the line of sight. Background clusters can be detected by red-sequence or redshift (with a spatial concentration of member galaxies), or mass S/N maps (especially after a photoz cut). We see no clear sign of background clusters that would strongly affect the mass measurements in the three cluster fields. Cosmic noise could contribute to 56 This can be estimated from the current mass error bars produced by resampling. If the source galaxies are circular before lensing (i.e., no shape noise), the variance of the fitted mass from resampling should be much smaller.
uncertainties as large as statistical errors for NFW profiles (Hoekstra 2003), but it can be reduced by using widefield data and the photozs of individual sources, and fitting the masses of the clusters along the line of sight (Hoekstra et al. 2011), and therefore we expect it would be smaller here. The 2-halo term can cause a bias at ∼5%-10% (Becker & Kravtsov 2011; Hamana 2011) on the mass measurement. In our case, A3911 and A3921 could be more affected than A85 because of their redshifts, and we found that adding a cut at ∼2r 200 gave a change of ∼6% on average. In the future, we will attempt to include the 2-halo term (e.g., Oguri & Hamana 2011;McClintock et al. 2019) and test the other profile models in the mass fitting. Finally, the projection effect of the halo triaxiality can lead to a scatter of ∼20% (Becker & Kravtsov 2011;Herbonnet et al. 2019).

Known Issues in Data Processing and Analysis Pipeline
There are a few known issues in our current pipeline; though they have small effects (depends on the science goal) on data products. New algorithms are being developed and will be implemented in future data processing and data release.
1. Local sky oversubtraction and global sky subtraction.
The default local sky background detection method causes background oversubtraction near bright and large objects and biases their photometry. But it has small effects on faint and small objects. We are experimenting with large-scale sky corrections on scales larger than a chip; the algorithm has worked successfully on HSC (Aihara et al. 2019), but its implementation on DECam is on-going. Also, the global sky subtraction can bias the shape measurement and photometry of small faint objects (Aihara et al. 2022;Li et al. 2022). For the purpose of the cluster-lensing measurement, the bias from the default method is not strong because most source galaxies are small, and in the future, we plan to produce different data products with and without the global sky subtraction to fulfill corresponding science goals. 2. Satellite trails. We find that the "persistent" sky/source (Aihara et al. 2019) determination for artifact rejection is not perfect especially near CCD gaps; however the streaks are usually filtered out in catalogs because they are flagged by "too many peaks." This issue can also be solved by viewing the images for coaddition and manually cleaning out visits that are strongly affected by satellite trails. 3. Joint calibration (jointcal). Transients and variable stars could potentially affect the photometric calibration; we expect these cases are at very low frequency. The number of visits is usually smaller at the edge of the target than at the center, and the imaging quality of the telescope is usually lower near the edge of FOV. These problems can be avoided by quality cuts on the coadd catalog. We also notice a small variance of the bestpolynomial-fitting order when we process the data of other clusters; though the variance has small effects on photoz. Currently the best order is determined by Figure 17. Counterparts for A3921. The median depth is u ∼ 24.5(24.2), g ∼ 24.9(24.6), z ∼ 23.5(23.2) at S/N = 10 and r ∼ 24.6(24.2), i ∼ 24.0(23.7) at S/N = 20 for point (or extended) sources.
comparing the jointcal-calibrated photometry with DES before coaddition; in the future, we plan to test the scatter width of stellar loci. 4. Galactic extinction. An extinction correction based on the SEDs of individual sources is expected to produce better results; because most sources lack spectroscopic observation, the inference of the source SED becomes essential (Galametz et al. 2017;Sevilla-Noarbe et al. 2021). For stars, the extinction can be overestimated for the ones that are inside or in front of the dusts (i.e., closer to the observer), and thus our photometric correction can be affected because those stars are overcorrected for the extinction and are actually less blue. But the MW dust is generally within 1 kpc (Schlafly & Finkbeiner 2011)-we expect this bias is small. Furthermore, most LoVoCCS clusters are selected to have low extinction. 5. Band response. The atmosphere and weather conditions can cause variance on the band response especially at blue bands (percent-level). Some of the effect can be averaged out during the calibration and coaddition; we did not observe huge differences when we compare our photometry with DES. 6. The selection of cosmological parameters. We tested rerunning the mass fitting with commonly used values h = 0.7, Ω m = 0.3, Ω Λ = 0.7, Ω b = 0.05, and we found the difference was 3%, which is much smaller than the error bar and was partially caused by the resampling. 7. As we mentioned earlier (Section 3.1), the correction to the B/F effect has not been implemented. Though the effect on shear measurement is at percent-level and hence small for the mass measurements of individual clusters, it will be significant in cluster ensemble studies. HSC and DES have both implemented the B/F correction (Bosch et al. 2018;Morganson et al. 2018), but there is no actual case of using the correction in obs_decam as a reference. We will seek to apply the correction to our pipeline in the future.

Summary
In this work, we first introduce the LoVoCCS survey by presenting its background, motivation, and design. Then, we give the target sample and observation strategy. Next, we show the details of our pipeline, which consists of the public LSST Science Pipelines and our own analysis steps, from the processing of raw DECam images to the analysis of lensing signal. We use two clusters, A3911 and A3921, which have not been well studied by lensing, and one cluster A85, which has previous lensing studies, to demonstrate our pipeline step by step. A3911 and A3921 are calibrated using SkyMapper for photometry, while A85 is calibrated by PS1 and SDSS; all three cluster fields are calibrated by Gaia for astrometry. After that, we show our data products and run the various tests to check the quality of our data; in general, the results are as expected, and we did not find significant biases. We then present our science results: A3911 shows agreement of masses and peak coordinates between lensing, X-ray, and SZ; A3921 Figure 18. Histograms of all source magnitudes in u, g, r, i, z in the A3911 field. The blue and orange curves describe PSF and CModel magnitudes respectively. For reference, the dashed and dotted lines ∝10 0.3m and ∝10 0.34m respectively, where m is a magnitude value. The 50% completeness limits are u ∼ 25.5(25.5), g ∼ 25.7 (25.6), r ∼ 25.9(25.7), i ∼ 25.6(25.4), and z ∼ 25.1(25.0) in PSF (or CModel) magnitudes (fitting 23 < m < 24). has discrepancies between lensing and baryon results because of its merging state; we obtain similar results of A85 compared to previous studies using different pipelines; all three X-ray luminous LoVoCCS clusters are at low redshift (z < 0.1) and have masses at the scale of 10 14 − 10 15 M e . Finally, we discuss the potential issues in our analysis.
A subsequent LoVoCCS Collaboration publication(s) will contain the ensemble results for ∼30 clusters, including offsets between lensing/baryon distributions, masses, substructure proportion and distribution, multiple-peak mass fitting, the gas fraction in the total mass, tomography (reduced shear as a function of redshift), and a first attempt of constructing scaling relations.
Forthcoming LoVoCCS Collaboration publications will present the analysis of specific science topics including jellyfish galaxies (their motion in the clusters and star formation); dwarf galaxies and LSB features, mergers and tidal effects, and their statistics versus dark matter distribution; star formation in the  Left: shape residual of the stars used for PSF-modeling in the whole A3911 field. The purple curve is a Gaussian with μ = 0 and σ = 0.003 for reference. i = 1, 2. Right: star-galaxy shape correlations in the A3911 field. The blue points, orange triangles, and green stars are C 1,2,3 correlations respectively. The value at distance 0 (i.e., self) is filtered out. The error bars are estimated by N s in each bin, where σ is the standard deviation, and N is the number count in each bin. We also added a cut of r CModel < 23.5 for galaxies to improve computation efficiency and reduce possible noise.  -Δδ plane, where δ cl is the cluster decl. The red contours show the density of the points using kernel density estimation. Right: histograms of the differences. The center peaks approximate to Gaussian. The green reference curve is a Gaussian with μ = 0 and σ = 0.02 rather than a fitting result. low-z cluster environment; multiwavelength studies to compare the optical results with X-ray, SZ, and radio observations; cosmological tests using the cluster ensemble data; synergies with other instruments or surveys, etc.
We will continue observing the remaining LoVoCCS clusters and processing their data in the next ∼2 yr; new processing and calibration algorithms and analysis tools will be developed and applied to future LoVoCCS data.
We thank the referees and Prof. Gregory Tucker and Prof. Jonathan Pober at Brown University for their comments on this paper. We thank the LSST Data Management team for teaching us how to use the software. We thank the members of the Observational Cosmology, Gravitational Lensing and Astrophysics Research Group at Brown University for contributing to the tests of the run_steps pipeline on different galaxy clusters. We thank the CTIO observation support team as well.
I.D. and D.C. are thankful for support from NSF (No. AST-2108287; Collaborative Research; LoVoCCS: The Local Volume Complete Cluster Survey). G.W. acknowledges support from HST program numbers GO-15294 and GO-16300, and grant number 80NSSC17K0019 issued through the NASA Astrophysics Data Analysis Program (ADAP). Support for program numbers GO-15294 and GO-16300 was provided Figure 24. Comparison between cluster red-sequence and the redshifts of cluster galaxies and field galaxies in the A3911 field (the cluster redshift is 0.0965). Left: the orange points stand for all galaxies. The color-coded circles are galaxies that have spec-zs. Right: galaxies that are color-coded by photoz. We only select z b < 0.3 for clarity. Figure 25. Counterparts in the A3921 field (the cluster redshift is 0.094). Note. We convert WL M 200c to M 500c assuming NFW profiles and the c-M relation in Child et al. (2018;fit: individual, all). by NASA through a grant from the Space Telescope Science Note. The clusters that have been processed are marked by " * ." The coordinates are from MCXC.

C.2. Spatial Photometric Comparison with DES
We show the spatial distributions of the differences between LoVoCCS and DES DR1 photometry without galactic extinction correction in Figures 32 and 33. The peak and median offsets are generally smaller than 0.01 mag, but we notice an offset 0.02 mag at Y band that needs further investigation (its depth can be the main cause).  . Stellar photometric comparison in the A3911 field. A spatial comparison is given in Appendix C.2. We select bright (g, r, i < 20; z < 19) unsaturated stars (>16 derived from base_PixelFlags_flag_saturatedCenter in the LSP and matches with SkyMapper DR2) in each band in both catalogs. The central peaks (fraction ∼93%) approximate to a Gaussian with μ = 0 and σ = 0.01; the scatter is due to the precision of dome flats and the depth of the reference catalog. The long tails possibly result from exposure depth, blending, and star-galaxy separation. Left: LoVoCCS (coadd) vs. DES DR1 (weighted average) stellar (PSF) photometry. Right: the photometric comparison with extinction removed. The r-band extinction is ∼0.03 mag near the cluster center.  figure) and with (the right figure) extinction removed. The results are similar to the ones of A3911. The r-band extinction is ∼0.06 mag near the cluster center. Figure 31. Stellar photometric comparison in the A85 field. We use the best-polynomial order for jointcal from A3911 and A3921 because A85 is not covered by DES. The figure shows the comparison of the stellar photometry between LoVoCCS and DECaLS DR8 (only includes g, r, z bands) with 16 < g, r < 20, and 16 < z < 19 in the A85 field after removing the median zero-point difference in each band (<0.03 mag). Figure 32. Spatial distribution of the difference between LoVoCCS and DES DR1 stellar photometry in the A3911 field. In the left panel of each image, the empty hollows are caused by the masking of bright saturated stars. Compared with the center of the field, the larger bias (0.04 mag) near the edge is caused by the relatively smaller number of visits for jointcal to calibrate; we expect its effect on photoz is negligible, and we only select the central patches for lensing analysis. In the right panel of each image, the blue histogram shows the distribution of the magnitude differences (summarized in Figure 29), and the orange curve is a Gaussian with μ = 0 and σ = 0.01 for reference rather than a fitting result. Figure 35. The A85 optical image overlaid with the E-mode mass S/N map contours. The top BCG is surrounded by X-ray contours. The offset between the BCG/Xray peak and the lensing peak is 3 ¢.  . We find similar patterns in our maps but using different processing and analysis pipelines.