PHANGS-HST: New Methods for Star Cluster Identification in Nearby Galaxies

PHANGS-HST has adopted the DOLPHOT photometry package (v2.0, Dolphin, 2000) as the principal source detection code in our pipeline. Specifically, we use the WFC3 and ACS modules, and main distribution tarball downloaded⁵⁵5http://americano.dolphinsim.com/dolphot/ on 4 Dec 2019. This choice was motivated by a preference to have a single, unified catalogue for selection of both candidate clusters and stellar sources (which are used to define multi-scale associations; see Lee et al., 2021; Larson et al., 2021). PSF-fitting source detection also takes full advantage of our high-resolution $HST$ imaging, de-blending closely neighbouring sources in an optimal manner.

Source detection using DOLPHOT is governed by an extensive set of parameters. A detailed list of the parameters used is given in Appendix B (Table 8). Briefly, we detect sources to the $3.5\sigma$ level using simultaneous PSF-fitting on the individual exposure (*flc.fits) images in all filters. The drizzled $F555W$ ($V$) band image is used as the positional reference for each target, meaning that exposure level detections in $F275W$, $F336W$, $F438W$, $F814W$ are joined into a unique source list with sky positions based on the drizzled $F555W$ world coordinate solution (WCS). Each detected source is photometered in all the bands, with non-detections flagged. We adopt the Tiny Tim PSF library (default within DOLPHOT) since the DOLPHOT-specific implementation of Anderson PSFs (Anderson, 2016, 2018) leads to systematic differences across the field (A. Dolphin priv. comm.). DOLPHOT makes refinements to the PSF and implements aperture corrections based on a subset of well-detected, isolated stars in each chip. Magnitudes are calibrated onto the Vega system using STScI-supplied zeropoints in the image headers. The photometric uncertainties reported by DOLPHOT tend to be underestimated, as they do not account for the noise contributed by neighbouring sources in crowded fields (see Williams et al., 2014, where it is shown that artificial star testing can be used to constrain the true uncertainties).

DOLPHOT is run using a set of python wrappers originally developed for the LEGUS project (by PHANGS-HST member L. Ubeda), and later adapted for PHANGS-HST (by D. Thilker). These wrappers organize the data prior to photometry, generate a complete parameter file specific to the data set, and manage execution of the DOLPHOT package codes.

We adopt the $V$-band DOLPHOT source list as our principal inventory of astrophysical sources (clusters and stars) in each observed target. Testing during pipeline development showed that occasional high-confidence clusters were not included in the DOLPHOT catalogue. In order to assess the degree to which DOLPHOT is effective at recovering clusters for our PHANGS-HST distance range, we match the DOLPHOT catalogue against the LEGUS Class 1 and 2 cluster population for NGC 628 (Adamo et al., 2017). We find that $\sim 98\%$ of LEGUS clusters have a DOLPHOT source within 2 pixels, and for $\sim 80\%$ the closest match is within 1 pixel. Given the methodological difference in detection algorithms (LEGUS used SExtractor; PHANGS-HST uses DOLPHOT), it is unsurprising that small positional shifts in the source positions are encountered at this level.

Because clusters with large angular sizes can be missed by DOLPHOT, we also run the DAOStarFinder code (in astropy/photutils), which is better suited to detecting extended objects. DAOStarFinder employs a convolution-based source identification method, whereas DOLPHOT requires explicit local maxima. We adopt a DAOStarFinder kernel with a FWHM of 2.5 pixels ($0.1\arcsec$, e.g. slightly broader than the PSF, so as to emphasise extended objects), and only add distinct new sources, plus those where the summed DOLPHOT catalogue flux within 2 pixels is $>2.5\times$ fainter than measured by DAOStarFinder. Typically only $\lesssim 1$% of all eventual candidate clusters originate from DAOStarFinder. The DOLPHOT catalogue augmented by DAOStarFinder sources is referred to as the PHANGS-HST “all-source” detection list.

DOLPHOT is used only for source detection in our cluster pipeline⁶⁶6The DOLPHOT PSF-fitted magnitudes are used in a parallel pipeline devoted to point sources.. Once the all-source detection list (DOLPHOT+DAOStarFinder) is ready, circular aperture photometry is performed using the photutils python package.

Because most clusters appear round in $HST$ images, we measure photometry in a series of circular apertures with radii of 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0 pixels. The background level is measured in a circular annulus 7–8 pixels away from the centre of each source, then subtracted off. This is a fairly standard background annulus when working in crowded fields imaged with $HST$ (e.g., Sabbi et al. 2018), and balances the need to measure the background level close to each source while being sufficiently distant ($\approx 14-29$ pc for typical galaxies in our sample) to minimise light from the source itself within the background annulus. We estimate the sky value as a robust median, using the sigma_clipped_stats function to mitigate the influence of neighbours. The sigma-clipped mean and standard deviation are also recorded by our pipeline, to allow for investigation of inevitable systematic and random errors in the adopted background. Calculations of the photometric uncertainty follow the prescription of F. Masci⁷⁷7Eq. 13 in http://wise2.ipac.caltech.edu/staff/fmasci/ApPhotUncert_corr.pdf with specific allowance for spatially-correlated noise in the pixels of the source aperture and sky annulus.

For sources selected as candidate star clusters (see Section 4), the total magnitude in each filter is determined by measuring the flux in a 4 pixel radius ($\approx 12$ pc at the median PHANGS-HST galaxy distance) aperture, which captures $\sim$50% of the total flux, and then applying an aperture correction (e.g., Adamo et al., 2017; Cook et al., 2019).

A detailed description of how PHANGS-HST aperture corrections are determined is given in Deger et al. (2021). Briefly, human inspection-based selection of an aperture correction sample (consisting of highly confident clusters, compact associations, and stars) in each target enables measurement of the correction for comparatively isolated sources. One of us (Bradley C. Whitmore, hereafter BCW) inspects the drizzled image products before our pipeline run and manually identifies 10 sources of each cluster Class (1, 2, 3). For each source, we compute aperture corrections to convert the $r=4$ pixel ($7{-}8$ pixel background annulus) $V$-band photometry to an estimated total (sky-subtracted at $r=20{-}21$ pixel) magnitude. Rather than letting the relative aperture correction between bands float on a galaxy-by-galaxy basis, which would add noise, we normalise our average aperture corrections as offsets with respect to the measured $F555W$ aperture correction in each target. These normalised corrections versus band are then combined into a mean survey-wide aperture correction offset table, after outlier rejection. For each galaxy we then adopt the final aperture correction as the galaxy-specific $F555W$ value plus the survey-wide offsets for other bands. The offsets are $-0.19$, $-0.12$, $-0.03$, $0.00$, $-0.12$ mag in $F275W$, $F336W$, $F438W$, $F555W$, and $F814W$, respectively. For the galaxies processed at this time, the $F555W$ aperture corrections range from $0.6$ to $0.85$ mag.

We refer to the complete photometric database resulting from these measurements as the ’all-source’ catalogue, to which morphological information (e.g. MCI) is later added (Sec. 4.2.1).

Synthetic cluster sources are added to the $V$-band image for each target galaxy, since these images are the ones used for source detection and cluster selection. The procedure is generalized sufficiently that we can easily add artificial sources to any PHANGS-HST band, although this is not currently implemented in the pipeline. However, we envision future work along these lines for the purpose of machine learning cluster classification training (e.g. Wei et al. 2020, Whitmore et al. 2021).

Random positions which are uniformly distributed are generated for the synthetic clusters. We considered using positions weighted to the observed source locations, but decided against it so that this step can be run independently from the rest of the pipeline. Note that if positionally-weighted properties are desired for the synthetic clusters, these can be selected post-facto. The location of a cluster influences detectability and photometric scatter/bias (with objects in complex and/or high background regions of a galaxy being harder to identify and measure than those in uncrowded and/or faint environments). Thus for future catalogue completeness analysis (beyond the scope of the current paper) we expect to impose such post-facto positional weighting.

Our fundamental goal in this portion of the pipeline is to generate a database of synthetic sources which adequately covers the plausible range of parameters expected for clusters in the PHANGS-HST imaging. This implicitly means we need to sample the age–mass plane, range of $A_{V}$ values and all plausible morphologies, and add enough realisations of clusters with a given set of parameters to fully represent the galactic environment and quantify the measured scatter in morphological metrics. To the degree that the entire cluster parameter space is represented in a target’s source population, the locus traced by morphological metrics most frequently seen in our synthetic database should resemble the observed distribution of confirmed clusters. Even so, we stress that our aim in the current paper is not to forward-model predict the expected distribution of metrics for a hypothetical cluster population, but rather to delineate the bounds of the distribution containing objects of interest.

We assign physical characteristics (age, mass, $A_{V}$, and morphology) to each synthetic cluster, to enable completeness studies (as a function of cluster type) and build a foundation for the forward-modelling framework alluded to above. We permit ages log($t$) = [6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0] yr and extinction values of $A_{V}$ = [0.0, 0.25, 0.5] mag (in addition to Milky Way foreground), plus a wide range of Moffat profiles which sample intrinsic size and radial shape; see Sec. 3.2 for details. The chosen limits of our model parameter space cover the vast majority of clusters found in nearby galaxies, with the particular exception of the very youngest, still embedded population. These particularly elusive, dust-enshrouded clusters (e.g. Lada & Lada, 2003) will be one focus of our approved JWST program (GO-2107) targeting half of our PHANGS-HST sample. Their expected numbers are uncertain and highly dependent on environment (Romita et al., 2016). For each combination of age, $A_{V}$, and morphology, we generate 1000 cluster realisations from a flat mass distribution between $10^{3}$ to $10^{5}~{}M_{\odot}$, sampling every $10^{2}~{}M_{\odot}$. Fixed step sizes in age, mass, $A_{V}$ are used for convenience. If we eventually further explore forward-modelling applications (as mentioned above), we would then likely switch to randomized sampling to eliminate any possibility of unwanted biases. Population synthesis models within CIGALE (Boquien et al., 2019) are used to predict the $V$-band magnitude for each synthetic cluster, based on its assigned age, mass, and $A_{V}$. In this way, we establish consistency between our synthetic clusters and the SED-fitting (Turner et al., 2021) of the PHANGS-HST pipeline. In total, we inserted some 4.5+ million synthetic clusters in each PHANGS-HST target, 200 objects at a time, simulating diverse cluster morphologies accounting for the PSF and incorporating Poisson noise on a pixel-by-pixel basis as described in Sec. 3.2 and 3.3. The same procedure used to determine photometry of real sources is applied to the synthetic ones.

We adopt Moffat profiles to represent our synthetic clusters, since these give a particularly good fit to the measured profiles of young clusters. This choice is well-supported by observational studies (e.g. Elson et al., 1987; Larsen, 1999; McLaughlin & van der Marel, 2005), although there are other reasonable options like King models (King, 1962, 1966). To reduce computational burden, our Moffat models are circularly symmetric, though in future work we may include asymmetric synthetic cluster models as they would be a superior representation of Class 2 clusters found observationally. Our model based selection (Sec. 4.2.2) is ‘loose’ enough that we do not expect to be biasing against Class 2 objects in the present work. However, if we later attempt forward-modelling of the cluster population, this complexity should be included.

The Moffat surface brightness profile is given by:

where $\mu_{0}$ is the central surface brightness, $a$ is a characteristic radius, and $\eta$ is the power-law exponent of the profile wings. Note that $\eta$ is equal to $\gamma/2$ in equation (1) of Elson et al. (1987). As $\eta$ approaches infinity the Moffat profile becomes a Gaussian. The characteristic radius can be expressed as a function of full-width-half-max (FWHM) and $\eta$:

Studies of clusters in nearby galaxies typically measure the effective radius $r_{\text{eff}}$, which is defined to be the radius of a circle containing half the integrated light of the model. For a Moffat profile, $r_{\text{eff}}$ is given by:

which is only valid for $\eta>1.0$. Although the effective radius is undefined for $\eta\leq 1.0$, and the integral of the profile infinite, the Moffat parameterisation can still be evaluated (and matching cluster morphologies observed) for this regime. In the real-world case it simply means that the cluster profile is truncated or lost to the background at some large radius.

The parameter space of circular Moffat clusters spans core size (via $a$ or FWHM) and halo power-law slope (via $\eta$). We adopted 8 values of $\eta$ (= 0.75, 0.875, 1.0, 1.25, 1.5, 2.0, 3.0, 4.0) and pre-PSF-convolution FWHM ranging from 0.5 to 7.0 pix, with step sizes of 0.25 pix. In total, taking all permutations across this parameter space, we allow for 216 (8$\times$27) distinct cluster morphologies. Models with high $\eta$ values (e.g. 3, 4) have minimal halo structure, and are dominated by their progressively Gaussian core. Conversely, low $\eta$ values correspond to very ‘fluffy’ cluster morphologies. Guiding the choice of minimum cluster core FWHM for our models, parametric fitting studies of marginally resolved stellar clusters have shown that reliable structural parameters can be obtained for effective radii approaching down to $10{-}20$% of the PSF FWHM, dependent on signal-to-noise and cluster morphology (e.g. Larsen, 1999; Ryon et al., 2017; Brown & Gnedin, 2021). We expect to have a morphological limit of $r_{\text{eff}}\sim 1$ pc at 16 Mpc (median PHANGS-HST distance, Anand et al. (2021)) with WFC3/UVIS.

Regardless of our particular choice for parameterising the structure of model clusters, we require a method to generate and realistically insert them into the observations in a way that mirrors the observational strategy (e.g. multiple exposures, dithering) and instrumental characteristics (e.g. resolution, noise). It is particularly important that the intrinsic cluster profile be convolved with the PSF. The assumed PSF is adopted from the empirically determined PSF products of J. Anderson, using a PSF representation appropriate to the centre of the WFC3/UVIS camera field of view⁸⁸8https://www.stsci.edu/~jayander/STDPSFs/ ⁹⁹9https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf. This choice allows the synthetic cluster database to be generated in advance of cluster detection/selection pipeline runs.¹⁰¹⁰10In later versions of our cluster catalogue we may adopt the PSF determined by DOLPHOT instead of generic representations.

Although it was developed to study galaxy morphology via parametric model fitting, IMFIT (Erwin, 2015) is well suited for our purpose of generating PSF-convolved models, especially as the package includes Moffat model capabilities. It is adopted for the initial steps of our cluster model realisation process (the stages of generating the intrinsic structure and convolving with the PSF), then a custom python script is used to do the remaining tasks, including adding noise, reprojecting exposure-level models to the drizzled grid, coadding, and finally adding the synthetic clusters to the PHANGS-HST images.

Images of real clusters may or may not be centred at the centre of a pixel. To account for this effect, which can become important for compact clusters, for each ‘batch’ of clusters that is inserted into an image, we adopt one of five different centres offset from the pixel centre by (0.0, 0.0), (0.25, 0.25), (0.5, 0.5), (0.25, 0.0), and (0.5, 0.0), to allow the centres to range from the centre to a pixel edge. Provided our step size is sufficiently fine, these offsets are sufficient to represent any actual cluster centre, as all other offset choices are equivalent via reflection and rotation.

For each batch of 200 synthetic clusters, of fixed age, $A_{V}$, FWHM, $\eta$, and subpixel location, we use IMFIT to generate a noise-free Moffat model using factor of 4 oversampling. The Anderson PSF appropriate to the filter and camera of the synthetic run (in our case $WFC3/UVIS/F555W$ almost always) is convolved with the model to produce a noise-free template synthetic cluster (still oversampled by a factor of 4) for each of our dither positions (generally three). Next we use these to generate native-scale ($0.03962\arcsec/\mathrm{pix}$) models (again one per dither position / $HST$ exposure) which were subsequently scaled to the integrated intensity corresponding to the adopted mass [magnitude, see beginning of this Section], independently per cluster. At this point we compute a Poisson variate for every significant pixel in the penultimate flc-exposure-level-appropriate model realisations, in a manner that accounted for the local area specifically around each cluster position as recorded in the unmodified image. Finally, we co-add the noisy models (per cluster) and add the result to the drizzled $V$-band image.

We begin with the set of sources detected by DOLPHOT, and impose a limit of $\mathrm{S/N}\geq 5$ and crowding $\leq 0.667$ mag (rationale on this in Sec. 4.4) for a source to enter cluster candidate processing. DAOStarFinder then adds in a very limited number of objects missed by DOLPHOT, meeting the conditions for augmentation described at the end of Section 2.1. Imposition of the DOLPHOT $\mathrm{S/N}\geq 5$ cut (even though our cluster candidate sample selection later demands $\mathrm{S/N}\geq 10$) excludes potentially unreliable compact sources from the determination of whether a given DAOStarFinder detection should be added to the all-source detection list, or if it was already preferentially well split into two or more DOLPHOT detections of reasonable quality. Next, aperture photometry and MCI computation is conducted for each source. From the resulting combined “all-source catalogue” (= bright+uncrowded DOLPHOT plus DAOStarFinder), we cull out poorly measured sources by requiring detections in $\geq 3$ bands with photometric error $\leq 0.3$ mag. We further require that V is one of these bands. These general criteria reduce the all-source catalogue to only $\sim 4{-}9\%$ of its initial entries. As described below, we further down-select this subset by considering two populations: (1) typical sources which meet a full set of morphological criteria (see Sec. 4.2) and (2) those that are so bright they are surely clusters or artefacts (brighter than the Humphreys–Davidson (H-D) limit, Humphreys & Davidson, 1979), for which we impose a more lenient set of criteria (see Sec. 4.3).

DOLPHOT’s goodness of (PSF-)fit metric ($\chi^{2}$) and sharpness parameter are both limited to the inner portion of sources ($\leq 3$ pixels for $\chi^{2}$, and the central pixel alone for sharpness), which is not ideal for selecting sources broader than the PSF. For this reason we do not use these quantities for candidate cluster selection.

One of the primary aims of our study is to improve on the cluster selection methodology used in previous studies, which has mostly relied on a single concentration index (CI, e.g. Holtzman et al., 1992; Whitmore et al., 1993; Adamo et al., 2017; Cook et al., 2019). High-resolution $HST$ images of clusters in galaxies at distances between $\sim 5{-}25$ Mpc allow for constructing maximally-informative (though still azimuthmally-averaged) constraints on cluster morphology, enabling more robust discrimination of point-like or cluster-like appearance than has been done to date.

We compute Multiple Concentration Index (MCI) metrics to probe radial profile shape over a wide range of distance from a source. These metrics allow us to compare the shape to that of a fiducial cluster morphology (via the specific choice of normalisation constants). In the subsections below, we formally define the MCI_in and MCI_out parameters used in our cluster candidate selection method, give examples of how cluster morphology varies with location in the MCI plane, and provide motivation for the area(s) of the MCI plane in which candidates are selected.

As noted in Sec. 4.1, we treat the brightest sources in our augmented DOLPHOT+DAOStarFinder catalogue slightly differently than the vast majority of (fainter) possible cluster candidates. The luminosity distribution for stars is strongly limited at the bright end. The Humphreys–Davidson limit (Humphreys & Davidson, 1979, a.k.a. H-D limit) effectively delineates a demarcation in the Hertzsprung–Russell diagram above which stars can only briefly exist, in an unstable state. This limit can be cast into an equivalent boundary in the ($V{-}I$, $I$) CMD. There is colour-dependence, but for our purposes we adopt $M_{V}=-10$ as the H-D limit. Clusters exist both fainter and brighter than this magnitude, but given the rarity of stars brighter than the H-D limit, it is reasonable to assume that nearly all sources brighter than the H-D limit are clusters. For this reason, we accept as a cluster candidate any source surviving the cuts outlined as general criteria in Sec. 4.1 and having $M_{V}\leq-10$ if it also meets the very lax condition of having $-0.55\leq\mathrm{MCI}_{\mathrm{in}}\leq 0.45$ (since they can sometimes be [nearly-]saturated or have close neighbours within 5 pixels, leading to corrupted MCI values). We call these objects ‘super H-D candidates’. They are marked as red symbols in Fig. 9. In practice, super H-D candidates are a mix of luminous clusters and artefacts. See the caption of Fig. 9 for more information, but note that contaminating foreground stars (being misinterpreted at the distance of the target galaxy) can also sometimes make it into this subset, and are later discarded to the best of our ability during classification. Objects of this sort should appear within or near the stellar exclusion region, however genuine luminous clusters having intrinsically compact morphology (e.g. Smith et al., 2020) would also manifest as super H-D objects found in the exclusion region – particularly for the more distant galaxies in the PHANGS-HST sample.

Point-like super H-D cluster candidates present a particular challenge during classification. We can utilize Gaia Bayesian distance estimates (Bailer-Jones et al., 2018) to weed out foreground stars, but such parallax-based ancillary information is not always conclusive since the method depends on sky location with respect to the Galactic model and because parallax uncertainty is excessive for some super H-D sources (e.g. our assumed H-D limit of $M_{V}=-10$ implies $m_{V}=20$ for a galaxy at 10 Mpc). The appearance of an Airy ring or sharp diffraction spikes can also help to weed out foreground stars. When neither Gaia nor PSF-wing source morphology helps, we rely on contextual hints from where in the galaxy image a point-like super H-D source is found. For instance, if it lies in a complex of active star formation or in a comparatively luminous region, and the angular surface density of other point-like super H-D sources nearby in the celestial sphere is low enough that chance projection on such areas is unlikely, then we generally accept the cluster candidate as being a bona fide cluster. A dedicated, follow-up study of super-HD clusters will explore this issue in more detail. Finally, one might wonder how complete our catalogues are for super H-D clusters. Whitmore et al. (2021) looks at this issue for NGC 628, concluding that our catalogues approach 90% completeness in this bright regime, despite the complexities of presented by point-like super H-D candidates.

As it is originally intended to be a point-source finding and PSF-fitting photometry routine, DOLPHOT sometimes chooses to represent an extended source as two or more very closely spaced point sources, rather than passing through only one dominant component. Although well motivated for resolved stellar photometry of crowded fields, and good for close pairs of stellar sources (as in the context of our PHANGS-HST individual star photometry and association analysis), this is a drawback in our case of cluster detection. We included a general condition on DOLPHOT crowding in our candidate selection to mitigate this situation from the start of our cluster-specific work. See Sec. 4.1 for specific implementation, but we only allow DOLPHOT sources with crowding $\leq 0.667$ mag. This crowding cut eliminates a majority of cases of duplicative detection, but at the very end of our cluster candidate selection we also explicitly impose a ‘doubles radius’ cut. After the preliminary candidate list is established, we sort candidates by count rate and work through the list in order of decreasing count rate, disqualifying any fainter DOLPHOT artefact (= ‘duplicate’) neighbouring candidate(s) at separations $\leq 2.5$ pixels (hence quasi-unresolved by $HST$). We cannot make the doubles radius any larger without starting to lose bona fide double objects (not artefacts) of which the fainter source could conceivably be a cluster.

DOLPHOT deblending artefacts that do squeak through this cut are later eliminated by human inspection, or, in the specific case of machine learning classified populations, via a secondary proximity cut to remove redundant links to the same physical object. The additional need for this secondary cut is best illustrated by Class 3 sources, which are operationally defined as groups of 4 or more point-like sources within 5 pixels radius (further discussion in Whitmore et al. 2021). Our detection code will generally return each point source in such groupings as a candidate, and subsequently most of these candidates (except perhaps those on the edge of a group) will be confirmed as Class 3 objects by machine learning. Even so, there is only one actual grouping. Therefore, after ML classification is completed (Sec. 5.2), we sort Class 3 by decreasing count rate and eliminate ‘redundant’ objects within a rejection radius of 5 pixels. Interactive inspection of the results showed that the method works as expected, keeping the brightest peak within each Class 3 grouping. The same operation is conducted for Class 1+2 classified candidates jointly – as a purely conservative step at this time – even though any pair of true clusters at separations $<5$ pixels would lose its fainter member. We expect to re-evaluate the Class 1+2 redundancy cut as the PHANGS-HST analysis progresses further. We do wish to clarify that Class 3 and Class 1+2 objects are not allowed to eliminate cross-class. That is, a Class 1 or 2 cluster can be allowed within the $2.5{-}5.0$ pixel range from a Class 3 compact association.

The simplest distillation of our cluster candidate selection method is to say that we evaluate each entry in our all-source catalogue (Sec. 2.1), first requiring that photometric and morphological quality assurance conditions are met (Sec. 4.1), then allow high quality sources to qualify as clusters either owing to their position in the MCI plane (Secs. 4.2.2 for the maximally inclusive, model-guided ML sample and 4.2.3 for the smaller, semi-empirically selected human sample), or because they are more luminous than the H-D limit (Sec. 4.3, included in both samples). Final checks to prevent any double counting are also implemented as just described. To aid the reader in understanding this complex selection method, we refer them again to the overview flowchart at the start of this Section (Fig. 1). The path in the flowchart from source detection to science-ready cluster catalogues ends with cluster candidate classification, which is the topic of the following Section.

Via interactive inspection of the $V$-band ($F555W$) image and also a colour composite made from $B$, $V$, $I$ bands, one of us (BCW) classified objects in the human sample. Inspection consisted of radial profile analysis and comparison of morphology with known stars during variation of the image intensity transfer function. A detailed description of the procedure is given in Whitmore et al. (2021). The large number of candidates in some targets prevented classification of all objects in this sample, and candidates were inspected down to a magnitude limit that included $\sim 1000{-}1500$ objects per target. The limit for the four galaxies in this paper ranged from $m(V)$ = 23.0 to 24.1 (see notes on Table 6 for details).

We applied the deep transfer learning ResNet-18 (18 layer residual, He et al. 2015, hereafter ResNet) and VGG-19_BN (Visual Geometry Group 19 layer with batch normalisation, Simonyan & Zisserman 2014, hereafter VGG) convolutional neural network models of Wei et al. (2020) to all of our candidate clusters, even those not included in the human sample. Specifically, we adopted the Wei et al. models trained using LEGUS-BCW human classifications for ten galaxies. As ML for cluster classification is an actively emerging field, future successful development will benefit from testing various approaches. We direct the interested reader to Grasha et al. (2019) and Pérez et al. (2021). We note that the method of Grasha et al. (2019) was also used for cluster classification required by the analysis of Messa et al. (2018a).

Our classification processing was accomplished on Amazon Web Services (AWS) cloud computing hardware, using GPU instances launched on demand. Running on a single GPU (NVIDIA Telsa V100-SXM2 16 Gb) we were able to attain a classification rate of $\sim$0.7s per candidate.¹⁴¹⁴14A step-by-step tutorial of our ML procedure is given in Whitmore et al. (2021), and linked at http://www.stsci.edu/hlsp/phangs-hst. For each of ResNet and VGG we evaluated ten independent models. From this set of results we obtained the mean, median, mode classification, plus standard deviation, for each network. The mode was adopted as the final single network classification. In Sec. 6.3, we experiment with various ways to attain a joint classification based on ResNet and VGG together. It is worth stressing that we do expect our ML classification accuracy to improve in the future, as we are starting further training experiments based on a representative set of BCW human classifications for PHANGS-HST candidates and from synthetic cluster populations, rather than relying on LEGUS classifications for typically more nearby galaxies than in our sample. As such, we do not consider the current ML classifications to be finalized.

In Table 5, for each analysed galaxy, we summarise information about the objects which are detected by our method and pass our selection criteria. As described in Table 4 we process NGC 628 as two separate fields without mosaicking them together due to camera/filter differences. Furthermore, we present the statistics for NGC 3351 broken into two columns because additional observations were obtained as part of PHANGS-HST, increasing spatial coverage to better match the ALMA observations, thereby increasing the exposure time in the region of overlap between LEGUS and PHANGS-HST imaging. Specifically, we analyse: (i) a LEGUS-only data set (NGC 3351-S in Tab. 4) and (ii) a mosaic constructed from all available images for NGC 3351.

Table 5 principally gives counts pertaining to our method, but we also include the total number of sources in the LEGUS ‘automatic’ catalogue and the subset which were classified by the LEGUS team. This allows an initial assessment regarding the efficacy of our new techniques.

Row (1), giving the combined number of DOLPHOT and DAOStarFinder sources, shows that the initial set of all-source detections coming from DOLPHOT is overwhelmingly large compared to the actual number of reliably identifiable stellar clusters. At least for the distance range we probe, this all-source census is dominated by stars, including some of marginal significance (as is customary for PSF-fitting resolved stellar photometry prior to selection of ‘good stars’ via quality assurance metrics). Row (2), listing counts after application of general photometric selection criteria (first paragraph of Sec. 4.1 and also Sec. 4.4), demonstrates the drastic reduction in the tally of potential candidates due to our limits on photometric error (in multiple bands). Row  (3) and (4) show the number of objects from Row (1) meeting our separated conditions of morphological appearance and luminosity, described in Secs. 4.2 and 4.3, respectively. In this table, the general conditions of Row (2) are not enforced for tabulating Rows (3) and (4), indicating the small fraction of the initial all-source catalogue having suitably accurate MCI measurements, and the minimal population of objects brighter than the H-D limit. Row (5) gives the number of sources meeting our general selection criteria and being selected as a candidate only due to being very luminous.

Our eventual selection, summarised in Fig. 1, is the result of demanding the set of general criteria and either the morphological or luminosity based criteria. We report the final number of cluster candidates so selected in Row 7 (ML sample) whereas Row 6 (human sample) gives the subset distinguished as occupying the semi-empirical selection region or being in the super H-D tail. Though not broken down in the table, most ($\sim 50{-}70\%$) candidates originate either within the semi-empirical selection region or the highest resolution model cluster region (0.01 binning), or their common zone. Candidates added due to their MCI plane location in the 0.02 or 0.04 bin-size model cluster regions are typically comparable in number to each other, each additional layer of area (at lower resolution) adding between $20{-}30$% to the candidate sample. Our most distant target, NGC 1566, is different, with the candidates originating more evenly across the MCI plane.

Rows (8) and (9) give the LEGUS candidate counts in total, and those with LEGUS classifications, respectively. Comparison of our ML sample size (Row 7) with that of the LEGUS automatic catalogue (Row 8) shows that despite cuts (e.g. Rows 2, 3, 4) our tally of candidate clusters is about a factor 2$\times$ larger than LEGUS. This is probably due to a combination of our allowance of more compact cluster morphologies as candidates ($\sim$4% with CI < LEGUS limit) and the substantial difference in source detection methods, since DOLPHOT deblends neighbouring sources more effectively than SExtractor. We find that candidates allowed into the sample only by virtue of their absolute magnitude (being too bright for a single star) are rare (Row 5), with generally between several to $\sim$30 in each target. NGC 1566 is an exception, with 130 super H-D candidates not otherwise qualified based on morphological metrics, presumably due to a combination of being much further away than the other galaxies (18 Mpc versus $\lessapprox$ 10 Mpc) and having nearly triple the HST-footprint-integrated SFR of any other target.

Drawing cluster candidates from the ‘human sample’ (see Table 5), one co-author (BCW) assessed source morphology and assigned classes following the method in Sec. 5.1. The depth of the human classifications varied with the specific target, going deeper in cases with few candidates and shallower in highly populous targets. We are able to use these classification data to assess the effectiveness of our selection criteria, and (in the future, with human classifications aggregated across all PHANGS-HST targets) to improve the training of ML networks. Because human classifications are confined to the portion of the MCI plane we call the semi-empirical selection region, with the exception of rare super H-D sources, we assess compact cluster selection purity ($\equiv$ # Class 1+2 / # Candidates) and Class 3 contamination ($\equiv$ # Class 3 / # Candidates) only within this area, leaving the outer model cluster regions beyond this for a later study.

Figure 13 shows the Class 1+2 selection purity attained in all targets, as a function of $V$-band magnitude. We present the information both as a computed purity with Poisson error (bottom row, red markings), and as the associated counts versus magnitude (top row, green and grey histograms). Purity generally increases from $\sim 35{-}50$% at 23 mag to $\sim 70{-}85$% brighter than 21 mag. In NGC 1566 the data suggest a flatter relationship. Note that our ability to accurately constrain the purity worsens as the population becomes brighter due to reduced counting statistics. Nevertheless, the overall trend across the sample is clear and the typical purity within the semi-empirical selection region is in excess of the comparable statistic for LEGUS (44%, 45%, 37%, 44%, 32% for NGC 628-C, NGC 628-E, NGC 1433, NGC 1566, and NGC 3351 (LEGUS-only), respectively). In Fig. 13, we also present a purity assessment for Class 1+2+3 (black markings in bottom row), finding that this metric including compact associations alongside compact clusters is typically $10{-}15$% higher.

Our method does not recover a majority of the Class 3 objects identified by LEGUS, but this is by design since we are treating these multi-peak compact associations with the dedicated method of Larson et al. (2021). Specifically, our selection areas in the MCI plane were determined on the basis of only Class 1 and 2 single-peak objects, both synthetic and observed. Class 3 objects generally lie below (more negative MCI_out at fixed MCI_in) Class 1 and 2 objects. Some compact associations (Class 3) do survive our selection as the selection areas for clusters (Class 1 and 2) had to be made broad enough to retain outlying Class 1 and 2 objects. Such surviving Class 3 objects should be considered contaminants, and we only keep them for completeness and comparison to prior work (LEGUS). We strongly advocate using the multi-scale associations of Larson et al. (2021) instead. This is particularly important if the emphasis is on studying young star formation products, as our Class 3 objects would be highly incomplete in such a context.

We achieved our goal of minimising Class 3 contamination (# Class 3 / # Candidates) amongst the objects considered as possible clusters. For the depth attained by PHANGS-HST human classification, such contamination ranged from 8 to 16% across the galaxies in our study. We do not compare to similar metrics for LEGUS, as they aimed to recover Class 3 compact associations in addition to Class 1 and 2 clusters (incidentally, this is a likely reason for their usage of SExtractor as a detection algorithm).

As currently implemented in our pipeline, we conduct ML classification using two different neural network architectures (ResNet and VGG), each pre-trained as described in Wei et al. (2020). For each network we evaluate the ensemble of class predictions ten times, allowing us to determine a mode cluster classification for each source. Just as with human classifiers, these modes do not always agree between networks or between networks and the BCW human classification. We are then faced with deciding how to best use the ML results, in the context of our goal to provide cluster classifications with the highest completeness (fraction of true Class 1 or 2 clusters classified as such) and lowest contamination (population of non-clusters classified as clusters). We have several options:

• (1)

  Identify the network which agrees with human classification of Class 1 or 2 most frequently. That is, select a ‘preferred network’ (ResNet or VGG) and adopt its output as the final classification.

• (2)

  Combine the modes of each network (per candidate) into a ‘mode consensus classification’, by accepting a candidate as a compact cluster (either Class 1 or 2, but unspecified) if either network predicts that and the other network does not predict Class 4 (star or artefact).

• (3)

  Evaluate the mode of all 20 (ten ResNet and ten VGG) predictions in a unified sense, to yield a ‘combined network classification’.

We now examine the outcomes associated with these options, and then use the difference between resulting compact cluster counts to: quantify the uncertainty of our ML classification methods, and to bracket the actual number of detectable Class 1 and 2 clusters in our data set.

Option 1 (‘preferred network’) is the most straightforward way to proceed, though it should be expected to provide lower accuracy than other choices – since the confusion matrices of Wei et al. (2020) suggest the accuracy of each network is approximately on par with the consistency achieved by a single human if classifications are blindly redone after a period of years. There is no redundancy against misclassification associated with specific downfalls of a network. However, it has the benefit of simplicity and retains distinction between Class 1 and 2 (the ‘mode consensus classification’ would not). In order to decide which network is preferred we compared the classifications for PHANGS-HST human Class 1 and 2 objects (independently per target) to ResNet and VGG classifications. We find that in 3 of 5 targets the number of identically classified sources is highest for the VGG network. For NGC 1433, ResNet is more frequently in agreement with human classification. The two networks are about equivalent versus Human classification for NGC 3351. Relaxing the set of PHANGS-HST Human classified sources to also include Class 3 (so 1, 2, or 3), we find that VGG classifications are more closely aligned with the Human assessment in 4 of 5 targets, and VGG and ResNet are about equivalent for NGC 1433. Therefore, we take VGG as our ‘preferred network’. Whitmore et al. (2021) independently come to the same conclusion, using different methods. Comparison of ResNet and VGG mode predictions for human Class 4 sources, shows that ResNet is very slightly more likely to agree with human classification that a contaminant is not a cluster. Table 6 lists the number of compact clusters (Classes 1 and 2) and compact associations (Class 3) identified in each of our targets (plus NGC 3351 LEGUS-only) via the preferred network, VGG. For completeness we also give the tallies according to ResNet.

Option 2 (‘mode consensus classification’) is an attempt to synthesise the classification output of both trained networks, and, in doing so, provide mitigation against misclassification by VGG or ResNet operating alone. As defined above, this option can be thought of as demanding that a source is agreed upon to not be a star or artefact (by both network architectures), and has been classified as 1 or 2 by at least one of them. As such, it is a balance between completeness (inclusivity) and contamination. Mode consensus classification results for all five targets are given in Table 6. The compact cluster counts are lower than for the preferred network option, but presumably benefit from reduced contamination.

Option 3 (‘combined network’), in which the mode is evaluated directly over set of 10+10 (ResNet+VGG) models, is also a means of synthesising all available runs of ML. It has the advantage over option 2 of retaining the discrimination between classifications 1, 2, 3, and 4. However, it is unclear if the possible drawback of adding information resulting from a non-preferred network (here ResNet) is overcome by the increase in statistical significance due to more ML classification votes. We include counts from the ‘combined network’ option to inform future work. In general, combined network cluster counts are between those of VGG and ResNet, or slightly lower than either. We note this ‘combined network’ option may have more meaning in a situation where networks (possibly even the same architecture) are trained according to different strategies, such as training with observed clusters (e.g. Wei et al., 2020) versus training with synthetic clusters and artificial stars. Discussion of other varied training strategies is given in Whitmore et al. (2021).

The choice of which option to use will ultimately be driven by the specific science use-case. If high completeness or retention of information regarding Class 1 versus 2 is important then we advise adopting the straight VGG mode classification. When forming distribution functions of cluster properties from such VGG classified objects, one could opt to weight by the fraction of VGG votes agreeing with the VGG mode, effectively placing emphasis on the most certain classifications. For those studies that require low contamination and can accept the inability to distinguish Class 1 from 2, we suggest use of the ‘mode consensus classification’. We emphasise that the census of clusters obtained in this manner will still be more complete than what is provided via PHANGS-HST human classification. Presently, we discourage use of the ‘combined network’ classification option until such time as diverse training strategies have been implemented.

Table 6 indicates that the ratio of Class 1 to Class 2 for ML (here VGG) is typically a factor of 2 to 4, but this same ratio is around 1 to 1.5 for human based tallies. The main reason for this is the brighter magnitude limit for human classification and does not imply systematic differences in how candidates are classified (between ML and human). For instance, human Class 1 clusters are generally brighter than ML Class 1 because of the limits imposed by inspection.

What can we say regarding the uncertainty of Class 1 and 2 cluster counts in our analysis? For each galaxy/field we interpret the ‘mode consensus classification’ tally as a lower limit on the number of Class 1+2 sources, and the maximum of VGG and ResNet Class 1+2 cluster counts (VGG and ReSNet individually) as an upper limit. Adopting the mean of these three tallies as a representative count, we can express the bracketed range as a percent uncertainty with respect to the mean. We find $\pm$13%, $\pm$9%, $\pm$22%, $\pm$18%, $\pm$16% respectively for NGC 628-C, NGC 628-E, NGC 1433, NGC 1566, and NGC 3351. Unfortunately, analogous percentage uncertainties from LEGUS analysis are not available for comparison.

It is vital to develop an understanding for potential differences between the Class 1 and 2 cluster populations identified by our new method and the techniques used by others. As noted earlier in the text, we chose to compare against the LEGUS survey, as several galaxies are in common between the projects and, even for targets that differ, equivalent five band $NUV$-$U$-$B$-$V$-$I$ $HST$ data are employed. We anticipated both random and systematic differences in the resulting catalogues, with the former associated primarily to classification uncertainty and the later linked to changes in source detection / candidate selection procedures. We did not undertake a head-to-head comparison of the Class 3 objects (compact associations) resulting from our cluster identification pipeline versus LEGUS, as our methods were tuned to preferentially identify single-peaked cluster-like objects.

We first generate a set of sources to guarantee a fair comparison between the surveys. We account for differences in photometry and coverage. LEGUS classified cluster candidates down to $M_{V}=-6$, effectively setting the faint end for comparison. Not surprisingly, slight differences in aperture corrections, photometric scatter (from minimal shifting of candidate centres), and occasional changes in aperture size suggest that in some cases we should back off from the $-6$ limit to ensure that objects brighter than $M_{V}=-6$ in PHANGS are also above this limit in LEGUS. The magnitude limits were allowed to vary by target and for machine learning classification-based comparisons versus human classification. The adopted values are given in Table 7. For our ML-based comparison the limit is generally up to $0.2$ mag brighter than $M_{V}=-6$, though LEGUS only classified sources down to $M_{V}=-7$ in NGC 1566. For our human classification-based analysis, the limit is even brighter for NGC 628-C, NGC 628-E, and NGC 1566, as these galaxies have very rich cluster populations and it was not practical for our expert human classifier (BCW) to inspect the entire population. Finally, of importance when setting the scope of comparison, for NGC 628-C and NGC 628-E the observational data were a mix of WFC3/UVIS and ACS/WFC imaging. These cameras have a different size on the sky and were not forced to have identical orientation. LEGUS only searched for clusters in the intersection of WFC3/UVIS and ACS/WFC coverage, whereas we allow cluster candidates anywhere with coverage in at least three bands, including $V$ ($F555W$). Accordingly we cut back the footprint of the survey comparison to the region with LEGUS classifications. Our eventual ‘comparison sets’ are defined as the union of Class 1 or 2 clusters from either PHANGS-HST or LEGUS meeting the magnitude and sky area criteria described above. Table 7 gives the number of sources in the comparison sets as the first row in each section of the table.

Within this comparison set for each target we tallied the number of Class 1 or 2 clusters: (a) found by PHANGS-HST, (b) found by LEGUS, (c) found by both surveys, (d) found by PHANGS-HST but not LEGUS, and (e) found by LEGUS but not PHANGS-HST. These cluster counts are listed in Table 7, along with their associated percentage of the entire comparison set. The percentages are also presented graphically in Fig. 14. In a panel devoted to each target and the type of PHANGS-HST classification (ML or human), the inner ring of each doughnut plot shows basic results of the survey comparison. Relative to the entire comparison set (think of this as anything that could possibly be a Class 1 or 2 cluster), the surveys agree at a significant percentage (case (c) above), approximately 50% for ML (median 47%) and slightly higher for human classifications (median 57%). Agreement is weakest for NGC 1433, with 38% for ML and 40% for human classification. Survey agreement also falls for NGC 1566, which we believe is a consequence of its substantially larger distance. This also makes sense in terms of NGC 1433, which Anand et al. (2021) have shown is actually at a distance of 18.6 Mpc rather than the 8.3 Mpc we assumed (in this paper only) for consistency with LEGUS. The percentages discussed above should not be interpreted in terms of completeness, since some objects will be identified as clusters by one survey but not by the other, and because the source list itself includes at least some non-clusters from both surveys. Rather, this $\sim 50{-}60$% agreement reflects that both LEGUS and PHANGS-HST are likely returning Class 1 and 2 cluster samples with at least this level of accuracy.

At face value, the higher and generally rather similar percentages of the comparison set reported for PHANGS-HST (case (a) – blue+green in Fig. 14) and LEGUS (case (b) – red+green) suggest that recovery of the true Class 1 or 2 population is indeed higher than the intersection statistics above require as a minimum. The median percentage recovered by the surveys for ML is 73% and for human classification is 77%, calculated by integrating over both surveys without preference. Of course, both surveys are likely subject to contamination and incompleteness. Contamination in particular will make the case (a) and (b) percentages skewed misleadingly high for the survey with the issue and low for the other survey. Without a more detailed look at the objects seemingly missed by either survey it is difficult to interpret any further than to say that both surveys most likely recover on the order of 75% of true Class 1 or 2 clusters meeting the magnitude/footprint constraints, and if systematic identification/misclassification exists then this figure will drop or rise accordingly depending on which survey has the problem.

Blue and red shaded portions of the outer rings of the doughnut plots in Fig. 14 aim to better understand the types of sources for which the surveys find different classifications. We begin with objects identified as Class 1 or 2 by PHANGS-HST but not by LEGUS (shown in blue). The wedges labelled L3 and L4 were classified as Class 3 and Class 4, respectively, by LEGUS (but Class 1 or 2 by PHANGS-HST). Approximately one quarter to one third of the PHANGS-HST-only Class 1 or 2 population falls into this category. We note that this happens more often in NGC 3351 than for other targets (42% for ML, 58% for human). In the case of human classification, this may be due to the greater exposure depth (which was not used in ML), in the overlapping LEGUS and PHANGS-HST fields. An even larger contribution ($\sim 40{-}70$%) to this misclassification category comes from sources which were somehow not even considered by LEGUS (labelled ‘no ID’). These result from intrinsic differences in the types of sources detected by SExtractor versus DOLPHOT, or due to automated cuts made by LEGUS after the detection stage (CI, number of bands with good photometry). A minority of sources, belonging to doughnut wedges labelled with asterisks, were at least considered as LEGUS candidates and seem to have met the conditions for LEGUS classification given by Adamo et al. (2017) but did not receive an eventual class for some reason we have not been able to determine. Next, we examine the sources deemed Class 1 or 2 by LEGUS but not by PHANGS-HST (red region in the outer doughnut). A first category is comprised of sources we flag as ‘doubles’ (see Sec. 4.4), where two objects are detected within 2.5 pixels of one another. Extended sources at such small separations, even if real and not a redundant DOLPHOT detection of the wings of the brighter extended object, would be challenging to classify accurately. We find that $10{-}20$% of the LEGUS Class 1 or 2 objects not included by PHANGS-HST fall into this category. Whitmore et al. (2021) also comment on the preferential inclusion of ‘doubles’ by LEGUS, and find they are typically assigned to Class 2. Differences in classification between surveys can also be seen as many PHANGS-HST Class 1 or 2’s are called Class 3 or 4 by LEGUS. In terms of Class 3 (see ML3, H3 doughnut wedges), they amount to $\sim$15%. One comment applicable to both survey classification efforts is the difficulty in classifying objects in crowded environments – this is particularly relevant to the distinction between Class 3 and Class 1+2. For example, imagine taking a confident Class 1 or 2 cluster and placing it in a region with many point sources or even other clusters. Based on the relative location of objects, it is possible that such a Class 1 or 2 is frequently judged by the classifier (either ML or human) as an apparent Class 3 compact association because the only evidence remaining of cluster-nature is the very inner profile of the source, which can be overridden or overlooked due to environment. We have no reason to believe this happens more often in either survey, and the ML3/H3 fractions are indeed similar to H3 in Fig. 14. A far more significant (and systematic) classification trend between LEGUS and PHANGS-HST is represented by the ML4/H4 wedges which dominate the LEGUS-only Class 1 or 2 population. In most targets more than 50% (and up to $>70$%) of LEGUS Class 1 or 2 sources not called 1 or 2 by PHANGS-HST are instead classified as contaminants by PHANGS-HST. We return to this point below. Finally, a small percentage of LEGUS Class 1 or 2 clusters are omitted from the PHANGS-HST Class 1 or 2 census for reasons we did not track down (red outer wedge labelled with asterisk), but possibly these sources failed cuts on MCI error.

In the case of the human comparison set, we have more detailed information on the appearance of sources LEGUS called Class 1 or 2 and PHANGS-HST classified as contaminants (that is, Class 4 in the original LEGUS system, and now subdivided by BCW into Classes 4.1, 4.2, …, 4.12 in order to specify the variety of contaminant cluster candidate, e.g. single star, stellar pair, stellar triple, saturated star, diffraction spike, galaxy nucleus, background galaxy, etc. – see Whitmore et al. (2021) for details). In Fig. 15, we show doughnut charts for these human Class 4 (‘H4’) sources, splitting them into up to six BCW subclasses (others of the 12 were not encountered). Specifically, the H4 sources are classified as a single star, pair or triple stars, redundant (another accepted Class 1, 2, 3 within 5 pixels), edge artefact, or too faint to classify. This categorical assessment of contributors to the H4 LEGUS-only population is most informative/reliable in cases having many sources. For this reason NGC 628-E and NGC 1433 should be viewed cautiously, as they each have $\leq 15$ H4 sources (only 8 in NGC 628-E).

The first observation to be made from Fig. 15 is that sources appearing to BCW as single stars are a minority in this population, amounting to less than 25% in all cases. The bulk of H4 sources are thought to be pairs or triples of stars according to the PHANGS-HST human classification. In general pairs constitute $\sim 45\%$, and triples $\sim 33\%$, of the H4 LEGUS-only population. Combined with the already flagged ‘doubles’ from the analysis above, we conclude that the majority of objects called Class 1 or 2 by LEGUS and not by PHANGS-HST are actually multiple stellar sources not meeting the source surface density needed to qualify as Class 3 associations. Whitmore et al. (2021) describes how this distinction is set by counting peaks within a 5 pixel radius. This effect is probably accentuated by situations in which SExtractor does not deblend such multiples and places a source between the actual peaks. This is a particular strength of using the PSF-fitting source detection algorithm (DOLPHOT) for our work, though it comes with the price of spurious detections in the wings of bright objects (largely eliminated by our doubles check, Sec. 4.4) and the need to explicitly eliminate redundant counting of associations due to proper deblending of subclumps (last paragraph of Sec. 4.4). The subclasses not yet described (redundant, edge, too faint) remain in the minority across the entire H4 population, with faint objects only cropping up in NGC 1566 and edge artefacts in NGC 1566 and NGC 628-E. The slightly more numerous ‘redundant’ population can be attributed to differences in source detection between surveys.

Lastly, we emphasise one final point from Fig. 14 concerning exact classification agreement. In the outer ring of the doughnut plots, in two shades of green, we indicate the fraction of Class 1 or 2 clusters in common between surveys that were given the same Class (e.g. 1 and 1; 2 and 2) with the darker wedge, and those that disagreed (e.g. 1 and 2; 2 and 1) with the lighter colour. Quite commonly the survey classification agreed on specific Class for this subset of the the overall comparison set. We direct the reader to Whitmore et al. (2021) for a detailed, extensive discussion of confusion between classes, including a discussion of confusion matrices for cross-survey (PHANGS-HST versus LEGUS) and internal (VGG versus ResNet, human versus VGG, etc.) comparisons.

Complementing the limited PHANGS-HST-ML to LEGUS-human classification results (top row of Fig. 14) of this Section, we direct the reader to the more extensive analysis of similar issues by Whitmore et al. (2021). In particular, Whitmore et al. (2021) introduce use of the colour-colour diagram as a ’figure of merit’ to further test the relative performances of the different methods.

This paper is not intended to provide astrophysical analysis of the cluster populations recovered by our methods, but rather to focus on the adopted procedures. Nevertheless, without presentation of some scientifically relevant plots, it would be impossible to fully appreciate the improvement we have attained relative to previous works, except in the sense of the comparison described in Sec. 6.4. In Fig. 16, we show the PHANGS-HST Class 1+2 clusters identified by human and ML classification for NGC 628-C. Fig. 17 shows the spatial distribution of clusters for the same field (trimmed to the WFC3/UVIS+ACS/WFC footprint) colour-coded by age in both contexts. These figures displaying our results for one target are both internally instructive, emphasising the more expansive cluster population recovered using ML classification relative to human-based efforts, and also can be compared versus the results of Adamo et al. (2017) who used the same data. Appendix A contains equivalent figures for all the fields analysed in our paper.

A limiting step in most studies of extragalactic star clusters is the manual effort needed to remove artefacts from candidate lists. It is difficult to automatically identify candidate star clusters with a high degree of certainty due to environmental complexity (effects such as crowding and high background), in addition to resolution limits. The MCI approach we introduced was originally designed to reduce the number of candidate star clusters that would need to be manually examined to remove artefacts such as pairs of stars, background galaxies, edge effects, etc. The basic strategy was to leverage photometry from several apertures rather than just the two used to determine the concentration index (CI), which has been the standard way to filter out artefacts. In this way, selection criteria could be enforced that insured that the profile followed a typical profile of a cluster rather than a star, pair of star, or other more complicated morphologies. Original trials using this approach were largely successful, reducing the size of the candidate list by about 50%.

During implementation, several other considerations reduced this percentage somewhat, while enabling other improvements to the overall strategy. One example is that the selection routine for both stars and cluster candidates has been unified to use a single primary method (DOLPHOT), rather than using two completely different methods employed by many other projects. This eliminates the necessity for subsequently cross matching the stellar and cluster catalogues after totally independent selection/classification decisions have been made, which would otherwise be a significant source of noise and complexity at a later stage of a project. However, a side-effect of the use of DOLPHOT is that it introduces a much larger initial candidate list since DOLPHOT’s selection routine uses an iterative method to optimally deblend neighbouring sources and find fainter, inconspicuous sources in the wings of bright objects. This is the main source of the growth of our candidate cluster lists. Another similar, but less important effect from a growth standpoint is the development of a more sophisticated method of removing point-like objects, allowing us to include proportionately more (of particularly compact) cluster candidates than possible using only the concentration index, hence improving the completeness levels.

The MCI approach also enables other potential improvements, many of which will be more evident in later stages of the PHANGS-HST project. For example, the development of cluster models needed to define selection regions for the MCI values will be useful in the analysis of various physical characteristics of the clusters, such as whether the objects are likely to be bound or unbound, whether certain kinds of potential profiles (e.g. very diffuse) are not actually found in nature, and how the profiles of clusters change as a function of age.

In summary, while some of the initial goals of the MCI method have not been fully realised (i.e. the reduction in the candidate list), the overall improvements to the project is evident and is expected to have long range benefits relative to the PHANGS-HST science goals and cluster studies in general.

Our fully integrated pipeline for cluster detection, identification, and ML classification produces reliable catalogues of clusters with minimal human effort. These catalogues have clear advantages over prior efforts, even using the same observational data and neglecting the relative ease with which they can be generated. This subsection is devoted to highlighting such scientific benefits in broad terms.

The advantage of reproducibility, and presumed replicability, should not be understated. Even within LEGUS, the premier cluster study prior to PHANGS-HST, attaining reproducibility was challenging (Pérez et al., 2021). The burden of human classification relative to the limited pool of trained team member classifiers ($\sim 12$) led to compromises in the uniformity of the adopted procedure. For instance, in general for each galaxy only a small subset (3 people) of classifiers examined the candidates (Adamo et al., 2017). This could impart systematic differences in the effective classification scheme between targets, if individual volunteers gave more weight to different aspects of the classification procedure. In cases where the number of candidates per target were very high, the classification task was divided across the galaxy, assigning stripes of coverage to different subsets of the classification group. It should be said that membership in each subset was shuffled, and the mode was adopted as final judgement amongst the classification votes per candidate, in hopes of minimising systematic effects. The adopted LEGUS method represents the best that can be achieved with a small group of human volunteers. ¹⁵¹⁵15Experimentation with citizen science based classification was undertaken by LEGUS as well, to potentially overcome the limits of person power, but the lack of interactive inspection tools (image contrast, radial profiles) on current platforms led to the conclusion this was not yet feasible. We note that in the more resolved case of the Magellanic Clouds, and M31 / M33 (with $HST$), citizen science cluster identification has proven rather successful (e.g. Johnson et al., 2012).

For PHANGS-HST, at no stage in the ML classification branch of our pipeline does human judgement enter the procedure in real time, so in one sense our analysis is reproducible. However, the training of the ML network embodies the non-reproducible¹⁶¹⁶16Whitmore et al. (2021) highlight the fact that even an expert reviewer, conducting classifications of the same objects separated by a period of several years, will assign a small fraction of clusters to differing classes. expert knowledge of cluster classification of co-author BCW, as well as the collective efforts of LEGUS classifiers (used for training by Wei et al. 2020 alongside BCW-only classifications). Because there are humans involved in the construction of the training sets, the outcome of the procedure may change if the training set is constructed by a different human[s], though Wei et al. (2020) conclude that LEGUS and BCW training sets seem to produce comparable accuracies. We expect ML classification accuracy to improve with further training method experimentation and as the database of human classified objects grows and becomes fully representative, both factors described by Whitmore et al. (2021). Systematics can still exist in our ML classification, dependent upon the overall suitability of the training sets to the actual objects being classified, and possibly on the degree to which the training sets are balanced among classes (clusters versus stars versus artefacts) or even among cluster morphologies within a class (physically large clusters versus compact). In regard to overall suitability, the majority of galaxies in the full PHANGS-HST sample are more distant than the LEGUS targets that were used by Wei et al. (2020). This is one of the main reasons that we still undertake human classification in our project, with the goal of accumulating a training sample that probes a more fully representative range of distance spanning the Local Volume (e.g. LEGUS) out to objects beyond the Virgo cluster (some of the PHANGS-HST targets). In essence, the construction of a training set should be seen as an integral part of the ML procedure, not as something that precedes it. To summarise, we hope that full adoption of ML classification by the star cluster community (see also important work by Pérez et al. 2021 and Bialopetravičius et al. 2019) will lead to replicability amongst groups that often have somewhat different data sets, and such replicability can build consensus regarding our understanding of cluster formation / evolution / destruction.

Secondly, and in a more practical sense, another methodological advantage is that by relying on ML classification we enable complete evaluation of all candidates down to our detection (selection) limit. Naturally, classifications of fainter objects can be expected to individually be less confident (distinguishing between Class 1 and Class 2 for example) than for bright analogues in the same environment (see Whitmore et al. 2021). However, integrated over an ensemble of many objects the cluster / not-a-cluster ML evaluation (that is, Class 1 or 2, or not) is still reliable near the faint limit of our candidate selection ($V$-band $\mathrm{S/N}\sim 10$), based on spot checks of deeper human classification versus the ML classification. We emphasise further that by classifying our complete list of candidates we are able to include the often substantial, easy to classify population of faint candidates occupying isolated, low surface brightness, low confusion environments. These typically get missed by human classification efforts, unless specifically remedied post facto (e.g. over 100 faint diffuse sources were added manually to the NGC 4449 sample in Whitmore et al. 2020).

Unfortunately, improved recovery of faint clusters enabled by ML is counteracted by losses in our ability to determine accurate photometric ages for young star clusters with mass $\lessapprox 10^{4}$ M_⊙ due to stochasticity (e.g. Fouesneau et al., 2014; Krumholz et al., 2015; Hannon et al., 2019; Whitmore et al., 2020). In some science use cases requiring single object characterisation, this may limit the usefulness of clusters at faint apparent magnitude, particularly as galaxy distance decreases.

Even so, ML classification for complete candidate samples dramatically increases the number of recovered clusters, translating into statistical improvement in the age–mass diagram, reaching significantly lower mass at fixed age, or alternatively greater age at fixed mass. In Fig. 16, we show the PHANGS-HST Class 1+2 clusters identified by human (open circles) and ML classification (dots). Class 2 clusters are colour-coded green. Notice how the ML population reaches masses about 0.5 dex ($\sim 3\times$) lower than the human classification. In the alternative perspective, at fixed mass (at least for intermediate age clusters), the ML catalogue probes up to a factor ten older in age. Appendix A includes age–mass diagrams for all the fields analysed in this paper. The degree of human-to-ML improvement varies within the five fields, owing to differences in depth of the human classifications.

There are numerous trickle-down scientific gains to expect coming from expansion of the identified cluster population. For example, improvements in the age–mass diagram will impact studies of the cluster mass function by allowing the bins in age that are typically used to each extend to lower mass (e.g. R. Chandar et al. in prep, note that the intermediate age clusters are often examined separately as they have survived removal from their natal environment). This will generally boost the statistical significance of mass function fit results through increased numbers, and may perhaps allow sensitivity to detect broken power law mass functions. Far better populated age–mass diagrams will yield improvement in modelling of cluster disruption (mass dependent versus mass independent) as it provides additional statistics for low mass clusters in an evolved, intermediate age state. We anticipate $\Gamma$ (the fraction of star formation occurring in bound clusters, e.g. Bastian, 2008; Kruijssen, 2012) determinations will increase, given that some of the newly recovered clusters are of Class 1 and have ages significantly larger than their crossing time (thus are likely bound at high confidence).

Another trickle-down scientific gain from the increased number of clusters is improvement of studies that use clusters as ‘clocks’ to estimate timescales for physical processes via cross-correlation analyses (such as cloud disruption, via comparison to the ALMA CO catalogues unique to PHANGS-HST). This is due to the increased surface density of the recovered cluster population in the ML census versus human classification outcome, as best appreciated in Fig. 17. Note the improved clarity of patterns in age in the ML cluster distribution. Appendix A presents the full set of ML classified catalogue spatial distributions. In the timescale ‘clock’ application, having a more complete cluster catalogue means that the power of apparent physical association between cluster and cloud will contribute to the cross-correlation function at the ‘correct’ length scale more frequently than otherwise in the limited human catalogue. Viewed in another way, thinking in terms of nearest neighbour analysis, incompleteness in a cluster catalogue might lead to the ‘wrong’ neighbour being selected. This being said, measurements (e.g. Cook et al., 2019) of $\Gamma$ in main sequence galaxies like those in the PHANGS-HST sample are generally $\lessapprox 30$%, and will almost certainly remain less than 50% even with expanded cluster catalogues. Therefore, in order to capture the majority of the star formation activity, we stress the importance of using other localized tracers such as the most compact multi-scale associations (e.g. Larson et al. 2021) as well for timescale ’clock’ applications.

Lastly, although discussed at length in the previous Section, we end the presentation of scientific benefits associated to our methods by reiterating the improvement we achieved with respect to excluding (often random) pairs and triples of point-like sources from the Class 1 and 2 census. This reduction in contamination is partially due to our specific attention to this issue but also due to our adoption of DOLPHOT rather than SExtractor as source detection code. Our choice to use DOLPHOT as a starting point also strategically facilitates a self-consistent inventory of point sources and clusters, and enables determination of multi-scale stellar associations (e.g. Larson et al., 2021). Fig. 18 illustrates the combined census of star formation products attained by PHANGS-HST by displaying the associations and clusters of NGC 628-C together. It is abundantly clear that a complete view of recent (few 10⁷ yr) star formation requires methods that capture structures spanning from very likely bound ($\Pi>$ many), possibly but indeterminately bound ($\Pi\lessapprox$ few), to unbound and already dissolved/disrupted. Another illustration of the innovative PHANGS-HST approach to combined recovery of clusters and associations can be found in Fig. 9 of Lee et al. (2021), highlighting a region of NGC 3351 with all components of our multi-observatory PHANGS program displayed. With the union of PHANGS cluster, association, CO cloud (from ALMA), and H ii region (from MUSE and traditional narrowband imaging) catalogues we will be able to place unprecedented constraints on the cyclic process of star formation in nearby galaxies.