Legacy ExtraGalactic UV Survey with The Hubble Space Telescope: Stellar Cluster Catalogs and First Insights Into Cluster Formation and Evolution in NGC 628^∗

We refer the reader to C15 for a description of the standard data reduction of the LEGUS imaging data sets, which are currently released at the Web site https://archive.stsci.edu/prepds/legus/. Each LEGUS target has a standard and homogeneous filter coverage provided by archival and newly acquired imaging with either the Wide Field Camera 3 (WFC3) or the Advanced Camera for Surveys (ACS). All of the galaxies have WFC3 imaging in the F275W and F336W filters. Three other bands cover the blue optical (ACS/F435W or WFC3/F438W), visual (ACS or WFC3 F555W or F606W), and NIR (ACS or WFC3 F814W). Hereafter, the conventional Johnson passband naming the UV, U, B, V, and I band will be used for the five filters, although the filter throughput is not converted to the standard Johnson system. We refer hereafter to the V band as the reference frame. The photometry provided by the LEGUS analysis is in the Vega magnitude system. Reduced science frames have been drizzled to a common scale resolution, which corresponds to the native WFC3 pixel size (004/px). The frames have all been aligned and rotated north up.

A description of the LEGUS imaging available for NGC 628 is given in Table 1. The data set consists of a mixture of archival and newly acquired data that have been reduced according to the standard LEGUS approach (see C15 for details).

2.2.1. Automated Steps

A custom pipeline, legus_clusters_extraction.py (version 4.0), has been developed to produce initial cluster candidate catalogs that contain, for each source, an ID number, the position in pixel coordinates, the final photometry including errors (in the Vega system), the concentration index (CI), and the number of filters in which the source was detected with a photometric error $\leqslant 0.3$ mag. A readme file includes all of the key information about the galaxy, the parameters used to build the catalog, and the content of the columns.

The pipeline comprises six consecutive steps and produces several diagnostics that help to fix key parameters to produce the cluster candidate catalogs for each galaxy. We present here a detailed description of each step in the pipeline and its application to NGC 628. This galaxy has been observed in two different pointings, NGC 628c (inner region) and NGC 628e (outer region). The two data sets have been analyzed separately, because of different combinations of cameras and filters (see Table 1).

• 1.  
  Step 1—Source extraction. The pipeline uses SExtractor (Bertin & Arnouts 1996) to extract source positions from the white-light image produced with the five standard LEGUS bands (see C15 for the method used to produce white-light images). In cases where the white-light image is not optimal (e.g., it shows numerous artifacts at the edges of the single frames, or the I band in low stellar density fields is dominated by red giants), we replace it with the reference V-band frame. The SExtractor input parameters are optimized to extract objects with at least a 3σ detection in a minimum of five contiguous pixels. These numbers change when the white-light or the V-band filter is used for detection. We do not apply any filtering to the image but we use small background cells to make sure that the strong gradient in the background between the arms and inter-arm regions does not affect our capability of detecting sources. These initial catalogs are visually screened to understand whether potential clusters have been missed and any possible improvement in the SExtractor configuration file can be made. The configuration file used for the source extraction will be released together with the final LEGUS cluster catalogs for each galaxy. At the end of this first step, we detect 6272 (NGC 628c) and 4539 (NGC 628e) candidate clusters.
• 2.  
  Step 2—Determination of CI and aperture radius parameters. In this step, the user selects a training sample of objects that is clearly stars (i.e., isolated, bright, and with steep luminosity profiles) and another sample that is clearly clusters (i.e., isolated, relatively bright, and with shallow luminosity profiles) using the reference frame (V band). The pipeline performs aperture photometry on these sources using radii of increasing size (from 1 px to 20 px with a step of 1 px). The local sky background annulus is set at 21 px and is 1 px wide. From this photometry, the pipeline calculates the CI (the magnitude difference between apertures of radius 1 pix and 3 pix) and the curves of growth for the input lists of stars and clusters. In Figure 2, we show an example of the plots produced by the pipeline in Steps 2 and 3. The top-left panel shows the CI distribution recovered for the star and cluster control sample. From this plot, the user selects a CI value that separates the visually selected stars and clusters as cleanly as possible. In the case of NGC 628c, the value we apply is 1.4. This value is an initial guess that can be iteratively adjusted in Step 3. As we will discuss and further investigate in Step 3 and in the completeness analysis, a selection criterion based on the CI cut allows us to remove a significant amount of stars from the catalog, but it will also remove compact clusters that appear unresolved on the frame (see Section 2.3 for a detailed discussion). The middle panel of Figure 2 shows the growth curve analysis on the stellar and cluster reference sample. The median stellar growth curve (thick red line in the plot) shows how the flux increases for an averaged unresolved source on this specific reference frame. This curve is linked to the stellar PSF produced by a given combination of detector and filter, and will change from galaxy to galaxy. The median cluster growth curve shows (thick blue line in the plot) how the flux is distributed in a resolved cluster spread function (CSF).⁴³ The increase in flux is slower in the resolved cluster than in the averaged stellar one, explaining why the CI is larger for resolved (clusters) than unresolved (stellar-like) sources (the CI method was first introduced by Holtzman et al. 1992; Whitmore et al. 1993). Inspection of the growth curves allows the user to fix the size of the aperture radius used to perform photometry. For each galaxy, we fix the aperture radius to do cluster photometry to the smallest integer number of pixels containing at least 50% of the flux within the aperture (aperture corrections for the missing flux are discussed in Step 4). We select the smallest aperture to reduce as much as possible the risk of contaminations from neighboring sources. We will refer to this value as the science aperture, to distinguish it from other apertures used in the analysis.
• 3.  
  Step 3—CI selection and multiband photometry. The CI distribution of all the extracted sources in Step 1 is plotted together with the CI cut value selected in Step 2. The top-right plot of Figure 2 shows the distribution for NGC 628c. As already discussed in the literature (e.g., Chandar et al. 2010; Adamo et al. 2015), the CI distributions of point sources that are unresolved (stars) will show a narrow peaked distribution (e.g., peak in the right panel of Figure 2). On the other hand, clusters are resolved and they show a broad range of sizes (e.g., Ryon et al. 2015), thus their CI distribution appears significantly broader than the stellar one (this can be seen in the distribution as the prominent wing extended to large CI values). The goodness of the chosen CI value for the cut is then directly tested on the reference frame (V band) via visual inspection of the objects in the catalog that satisfy this condition. We check whether a smaller CI introduces a large contamination of stars. The chosen best CI value is, eventually, indicated as a selection criterion. All sources with a CI smaller than the reference value (1.4 for the ACS V band of NGC 628c and 1.3 for the WFC3 V band of NGC 628e) are removed from the initial catalog. We then perform multiband aperture photometry on the sources that pass this selection step, using the science aperture fixed in Step 2 and a local sky annulus located at 7 px with a width of 1 px. The same science aperture and local sky annulus is used in all the five filters.
• 4.  
  Step 4—Averaged and CI-based aperture correction. We apply two different aperture correction methods both widely used in the literature.In the first case, we produce averaged aperture corrections using the cluster control sample produced in Step 2. The correction is estimated as the difference between the magnitude of the source recovered at 20 px (sky annulus at 21 px, 1 px wide) minus the magnitude of the source obtained using the science aperture (sky annulus at 7 px, 1 px wide). Not all of the sources in the control samples are detected in all the filters. To avoid clusters not detected in some filters from skewing the averaged values of the aperture corrections, we reject sources that have corrections outside a certain range. We show for example the aperture correction distributions recovered for the cluster control sample of NGC 628c in the five LEGUS band in the bottom panels of Figure 2. The vertical dashed lines show the limits within which we use the values to produce the average correction. These limits can be adjusted by the user for a given galaxy. In the case of NGC 628c, some of the sources are very faint or not detected in the bluest filters; therefore, their corrections become extreme. In Table 1, we list the final averaged aperture corrections and the errors (standard deviations) for each filter of each pointing in NGC 628. These values are added to the photometry produced in Step 3. The standard deviation of the aperture correction recovered from the sample is added in quadrature to the photometric error to take into account the uncertainties introduced by the averaged correction. With this method, the differences in the sizes of the clusters are not taken into account. However, as one can see from the recovered values in Table 1, the differences in the corrections as a function of waveband are quite small. Therefore, averaged aperture corrections do not change the shape of the spectral energy distribution (SED) of the source—they only change the normalization. This will slightly affect the mass estimates (not the derived ages and extinction values), but it will be well within the average 0.1 dex uncertainties that the SED fitting methods produce.In the second approach, we produce aperture corrections based on the CI of the source in each band. The method to derive the CI versus aperture correction relation will be described in a forthcoming paper (D. O. Cook et al. 2017, in preparation). This method has initially been developed to produce an aperture correction for very extended clusters (partially resolved in their stellar components) detected in the LEGUS dwarf galaxies and afterwards extended to the rest of the LEGUS targets. The CI-based aperture correction for each source is calculated by measuring the CI in each band and applying the CI–aperture correction relation determined for the appropriate instrument and science aperture size. While this method takes into account the size (indirectly measured via the CI) of the source as a function of waveband, it can change not only the normalization of the SED but also the shape (so it will affect the mass, age, and extinction). This effect will be particularly large in faint sources and in crowded regions where the estimate of the CI becomes more uncertain.At the end of Step 4, we have two sets of aperture corrections calculated. Within the LEGUS naming convention, the two final photometric catalogs produced in Step 5 with these two correction methods and derived analyses are referred to as AV_APCOR and CI_BASED.
• 5.  
  Step 5—Final products. The quality of the photometry is checked in this final step. We remove all of the sources that have not been detected in at least two filters (the reference V band and either the B or I band) with a photometric error ${\sigma }_{\lambda }\leqslant 0.3$ mag. In each band, we correct the photometry of these sources by the foreground Galactic extinction (Schlafly & Finkbeiner 2011). Two automatic final photometric catalogs are produced, e.g., the automatic_catalog_ngc628c_avgapcor.tab and the automatic_catalog_ngc628c_cibased.tab. Both catalogs contain the ID of the sources, positions in pixels and in R.A. and decl., multiband photometry and errors, CI, and a flag indicating the number of filters in which each source was detected. In the case of NGC 628c and e, the automatic catalogs produced at the end of this step contain 3086 and 593 cluster candidates, respectively.These sources, however, are likely not only star clusters. Among them there are still interlopers, i.e., pairs and multiple stars in crowded regions, background galaxies, bright stars in the galaxy that have CI slightly larger than the limit used here, foreground stars, and objects and artifacts at the edge of the science frames. To minimize the contaminations of our final catalogs, we visually inspect a fraction of the sources found in the automatic catalogs that satisfy some extra selection criteria. We select for visual inspection all the sources that are detected in at least four bands (in the automatic_catalog_ngc628c_avgapcor.tab) with a photometric error below 0.3 mag and have an absolute V-band magnitude brighter than −6 mag. The V-band magnitude cut is applied to the V-band photometry obtained with the CI-based aperture correction. This allows us to include slightly fainter, very diffuse clusters since they have large values of the CI and therefore get a larger aperture correction.
• 6.  
  Step 6—Missed clusters. During the visual classification of the sources it is also possible to add sources that have not been included in the visual inspection catalog. If such extra sources are not found in the automatic catalogs then we use Step 6 to produce their final photometry in the same homogeneous way as done for the bulk of the other sources. The number of added sources is typically below 1% but it can change from galaxy to galaxy.

Zoom In Zoom Out Reset image size

2.2.2. Classification of Cluster Candidates

At this stage of the reduction, three or more independent classifiers from within the LEGUS team visually inspect the cluster candidates. A tool has been made available within the LEGUS collaboration to perform this task. This tool visualizes the object to be classified in two frames, in the reference V band and in a combined three-color image. Using an interactive window, the user can also visualize the contour, the radial profile, and the surface plot of each source in the reference band (see Figure 3) using the IRAF task IMEXAMINE. Two circles on the images show the aperture used to perform photometry (inner) and the location of the local sky (outer circle). Each source is inspected in the single-band image to check whether the light is extended (cluster) or has a stellar PSF (star). The difference in the light distribution of a cluster and a star at the distance of NGC 628 can be easily discerned in the top row panels of Figure 3. The source at the center of the inset is a cluster and it has an FWHM = 3.71 px. A star is visible in the same cutout and its light is much less extended (FWHM = 2.4 px) than the one of the cluster. The contour and surface plots (middle and right panels) confirm the different behavior in the two objects.

Zoom In Zoom Out Reset image size

Based on this inspection, each object in the catalog gets assigned one of four defined classes. In Figure 3 we show an example object for each class taken from the central pointing of NGC 628. Class 1 clusters are compact and centrally concentrated objects, with an FWHM more extended than the stellar one. They show a homogeneous self-consistent color. With respect to class 1, class 2 systems are clusters with slightly elongated density profiles and less symmetric light distribution. Class 3 can potentially be less compact and homogeneous clusters, more likely compact associations⁴⁴ , and they show asymmetric profiles and multiple peaks on top of diffuse underlying wings, which suggest the presence of a possible concentration of low mass stars. Finally, class 4 objects are single stars or artifacts, or any other interlopers (e.g., background galaxies, foreground stars) that can affect the cluster analysis. The bottom row of Figure 3 shows an example of an object detected with extended CI, morphology, and CSF profile similar to a class 3 object. However, the stars belonging to this system have clearly different colors. We consider these objects as a chance superposition along the line of sight of two stars; thus, they are excluded from the cluster catalog.

Within the LEGUS collaboration, we are also testing a machine-learning (ML) recognition of the sources. The visually inspected catalogs are fed to the algorithm as training sets for the subsequent classification of sources. Testing is in progress and the results will be published in a forthcoming paper (Grasha et al. 2017, in preparation). The cluster catalogs of a third of the LEGUS galaxies have currently been inspected by humans. We are providing visual classification of a fraction of the catalogs of another third of galaxies, which are used as training sets by the ML algorithm. The goal is to perform the visual classification of the remaining galaxies using a purely ML-based approach.

To investigate the impact of the assumptions and the selection criteria adopted to produce the LEGUS cluster catalogs, we have performed several completeness tests. The tests are run with a custom-made pipeline, legus_cct.py (LEGUS Cluster Completeness Tool), available within the LEGUS collaboration. The pipeline runs in five consecutive steps, i.e., creation of synthetic sources, source extraction, photometry, recovery fractions, and final outputs. One of the input files is the same used to run the standard legus_clusters_extraction.py. The second input file specifies quantities necessary to create the synthetic frames containing the simulated clusters, e.g., input parameters for the BAOlab software (Larsen 1999). A stellar PSF generated with the IRAF task PSF using selected stars in the science frames is also provided. This PSF is then convolved with MOFFAT15 profiles of specified effective radii, R_eff. The resulting extended sources are then randomly distributed in a blank frame of the same dimension as the science one. A magnitude range is provided as well. To overcome crowding problems, we only insert a thousand clusters per loop. The frame containing synthetic clusters is then added to the real science frame. In the next step, all the sources are then extracted using the same procedure as in legus_clusters_extraction.py. For all the sources extracted, it estimates the CI and the photometry (including averaged aperture correction) in the same way as for the real clusters. Using the known position of the simulated clusters, the software estimates a source recovery fraction before and after the CI cut is applied.

As a first result of these simulations, we have investigated the relation between the CI and R_eff of increasingly extended sources, with a particular emphasis on the impact that the CI cut has on removing compact clusters as a function of distance of the galaxy. In Figure 4, we show the recovered relation between CI and R_eff for the ACS and WFC3 F555W science frames, used as reference frames for the CI selection. As an example, we refer to the study case NGC 628. The applied CI cuts of 1.4 mag and 1.3 mag for the ACS and WFC3 F555W correspond to a cluster ${R}_{\mathrm{eff}}\sim 1$ pc. This does not imply that all clusters of 1 pc size are systematically removed from the catalogs. We investigate the recovery fraction of sources in the V-band frame before and after the CI cut is applied. In Figure 5, we illustrate the results for NGC 628c. Before the CI cut is applied, we have 100% recovery for sources of 1 pc as well as more extended ones. However, the recovery of very compact sources goes below 50% after the CI cut is applied. Only half of the sources with ${R}_{\mathrm{eff}}\sim 1$ pc are extended enough to make it into the selection. A smaller CI will include too many stellar objects, so we conclude that at the distance of NGC 628, the CI cut applied removes a fraction of unresolved clusters with sizes of 1 pc and below. This choice does not introduce biases in our cluster analysis because observationally, no clear relation between the size and mass or luminosity of star clusters has been found (see Ryon et al. 2015 and references therein). Moreover, an ongoing analysis of cluster sizes in NGC 628 shows a clear log-normal distribution which is well above the detection limits of 1 pc (Ryon et al. 2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 6, we illustrate the final completeness limits of the two pointings in NGC 628. In Table 1, we list the magnitude limits corresponding to 90% recovery of sources with a detection threshold above 3σ (column 8 Table 1) within a minimum of five contiguous pixels. For the same band, the differences in the exposure times (longer exposures) have a larger impact on the recovery fractions than the differences in crowding between the inner and outer region (i.e., compare the recovery magnitude for the F336W and F435W for the inner and outer pointing). The completeness software has been run independently for each band. Therefore, these values must be considered as indicative of the detection limits intrinsic to the science frames. However, we are not taking into account that we apply to our catalogs additional selection constraints (detection in at least two filters, or in the case of visually inspected sources in four filters), visual inspection, and that real clusters have different luminosities at each wavelength depending on their age, mass, and extinction.

Zoom In Zoom Out Reset image size

Figure 7 shows a comparison between the LEGUS final catalog (green circles are class 1 and 2, blue circles class 3) and a catalog from Whitmore et al. (2014, shown as red circles) for a small region in NGC 628. The Whitmore et al. catalog was produced automatically based on Hubble Legacy Archive (HLA) observations, including measurements of the concentration index and algorithms that remove close pairs and multiple hits in very crowded regions. A relatively bright limit is used (${M}_{V}=-8$) for the HLA-based catalog since no manual inspections are made.

Zoom In Zoom Out Reset image size

The correspondence is fairly good, with 77% of the HLA-based sources being included in the LEGUS catalog. Making a similar magnitude cut in the LEGUS catalog (which goes roughly 2 mag deeper) results in a recovery rate of 73% of the HLA-based catalog being found in the bright end of the LEGUS catalog. There are two primary differences: (1) the LEGUS catalog includes approximately five more (out of about 100 in the region of overlap used for the comparison) bright class 1 objects that are clearly good cluster candidates and were missed in the HLA (probably due to using too conservative a limit for the concentration index), and the LEGUS catalog includes about 15 more compact associations—i.e., type 3—than found in the HLA-based catalog. Although selection effects are an important consideration in any study, we note that Chandar et al. (2014) compared completely manual, hybrid (like with LEGUS), and completely automatic catalogs made by two different studies and found that the selection did not result in major changes to the primary results of the studies (i.e., the CLF, CMF, and age distributions).

In this section we compare the photometry and analysis performed on the CI_BASED catalog to our reference catalog. For convenience, we show only the analysis performed on the NGC 628c pointing, but the outcome is very similar for the outer field.

The advantage of using average aperture corrections is that this method does not change the shape of the SEDs or the intrinsic colors of the clusters. As we see in Table 1, the differences between the applied aperture corrections change only slightly between different bandpasses (i.e., it reflects the change in the PSF as a function of filter/camera). Therefore, such a correction will change the normalization of the cluster SED, but not the shape. This method, however, does not take into account that, at the distance ranges of the LEGUS targets, clusters are partially resolved, that is, their CSF changes as a function of the cluster size. However, since the intrinsic shape of the SED is preserved and only the normalization is affected, this method may overestimate the mass of very compact clusters or underestimate the mass of the very extended sources, but the uncertainties will be within the 0.1–0.2 dex in logarithmic age and mass usually produced by the fitting method.

The CI-based aperture correction takes into account the relation between the size of the cluster and the required aperture correction. The relation has been derived in each band using simulations of YSCs of varying sizes (see D. O. Cook et al. 2017, in preparation). This method has the great advantage of taking into account the cluster sizes. The limitation, however, resides in the uncertainties produced by the CI estimates as a function of wavelength. This means that the CI-based method not only changes the normalization of the SED but also the shape. Therefore, the uncertainties propagate in the estimated mass, age, and extinction of the source.

In Figure 10, we show the same UV and optical color–color diagrams but derived using the photometry of the CI_BASED catalog. We notice that while the overall location of the clusters and associations are similar in the two photometric catalogs, the contours of the CI_BASED photometry are more extended on the left side of the tracks, a spread that is larger than the photometric uncertainties. In general, one expects sources identified as clusters to diffuse on the right side of the evolutionary tracks because of reddening, while the spread on the left side is mainly produced by photometric errors and uncertainties in the calibration.

Zoom In Zoom Out Reset image size

In Figure 11, the analysis of the residuals (i.e., the difference between the observed and the best model-integrated fluxes) of class 1, 2, and 3 clusters in each band does not show significant differences in the residuals of the F336W filter in both the AV_APCOR and CI_BASED catalogs. The residuals of the F275W filters in the CI_BASED catalog are more concentrated around zero (the best match) than the AV_APCOR ones. The opposite trend is observed in the AV_APCOR residuals of the BVI filters, where they show a narrower distribution around the best match than the CI_BASED ones. The latter are less peaked and show a tail toward positive values (i.e., the observed magnitude is brighter than the one predicted by the best model) in the B and V bands and negative values in the I band. Overall, we see that the differences between the reduced ${\chi }^{2}$ obtained from the fit performed to the photometry of the AV_APCOR and CI_BASED catalogs show a more significant negative tail, suggesting that the fit to the CI_BASED photometry is slightly worse.

Zoom In Zoom Out Reset image size

A direct comparison of the recovered ages, masses, and extinctions of class 1, 2, and 3 objects of the AV_APCOR and CI_BASED catalogs is shown in Figure 12. Although on average we see correspondence along the one-to-one line, some deviations are also significant. The most important difference in the CI_BASED catalog is the appearance of the pronounced peak at around 5 Myr, visible in the top-right distribution of Figure 12. According to the top-left panel, these systems have been assigned in the AV_APCOR catalog both younger, similar, or older ages. Since the mass-to-light ratio is smaller at younger ages, the objects that become younger in the CI_BASED analysis will also have smaller stellar mass; therefore, we can explain why the AV_APCOR and CI_BASED mass histograms (central right panel of Figure 12) differ in the mass bins between 1000 and 10⁴ ${M}_{\odot }$. In the CI_BASED catalog, these objects have been assigned masses $\leqslant {10}^{3}\,{M}_{\odot }$. In total, we estimate that about 20% of the class 1, 2, and 3 systems that are in common in the two catalogs have differences in ages larger than 0.1 dex, which is the average uncertainty recovered for age estimates obtained with deterministic methods. This small fraction of deviating sources may explain why we do not see significant differences in the recovered CLFs, CMFs, and disruption rates (Section 5) of the cluster population in NGC 628 when AV_APCOR and CI_BASED catalogs are used. Further investigation of these two approaches is presented in D. O. Cook et al. (2017, in preparation). In the next section, we will report the results of the analysis performed on the YSC population of NGC 628 using only the AV_APCOR reference catalog.

Zoom In Zoom Out Reset image size

The luminosity function of YSCs is typically described as a power-law function in the luminosity space, ${dN}/{dL}\propto {L}^{-\beta }$. However, we fit the function in logarithmic space, where the luminosity is replaced by the magnitude, so that $\mathrm{log}[d(N)/{dM}]\,\propto \theta \times M$, where $\beta =2.5\times \theta +1$ (e.g., Whitmore et al. 2002; Haas et al. 2008).

We apply two different techniques, equally used in the literature, to analyze the CLF of NGC 628. In Figure 13 we report the observed luminosity functions in each band, using cumulative distributions (left panel), like in Bastian et al. (2012), and bins containing the same number of objects (right panel), following the Maíz Apellániz & Úbeda (2005) approach.

Zoom In Zoom Out Reset image size

The luminosity properties of the clusters are directly observable, not affected by age or mass determinations. However, their luminosity distributions depend on both the detection limits of our data sets and by the adopted extraction procedure combined with the selection criteria we impose to yield the final catalog. To build the CLF, we select only sources that have been visually classified as class 1, 2, or 3 (thus they are brighter than −6.0 mag in the V band and detected in four filters with photometric error smaller than 0.3 mag) and are younger than 200 Myr. In Figure 14, we show the age–mass diagnostic diagram including the recovered 90% detection limits in the four bands required for detection and the V-band cut at −6.0 mag, all converted to limiting masses as a function of age. We observe that the V-band cut applied for the visual inspection is more conservative than the detection limits of our data set. However, we notice here that the resulting flattening of the distributions at the low luminosity ends is produced by the combination of a sharp magnitude cut combined with detection limits of the science frames and the method used to produce the final position catalog of cluster candidates. The age limit of 200 Myr enables us to directly compare the CLF to the CMF (see Section 5.3). In total, we count 733 (370) class 1 and 2 (class 3 numbers are indicated inside brackets) objects in the F275W filter, 846 (404) in the U band, and 851 (408) in the BVI bands before the age cut. After the age cut is applied, we are left with 703 (369), 778 (397), and 783 (401) class 1 and 2 (class 3) objects in the UV, U, and BVI bands, respectively. To prevent incompleteness from affecting our analysis, we perform the fit of the binned and cumulative distributions from the brightest bin (object) down to the bin (object) with magnitude brighter than 22.12, 22.26, 23.50, 23.50, and 23.0 mag in UV, U, B, V, and I, respectively. These limiting values have been chosen to avoid the shallower regions of the distributions and are more conservative than the 90% completeness magnitudes reported in Table 1.

Zoom In Zoom Out Reset image size

In the case of an equal number of object bins, the slopes and the associated uncertainties are produced by the IDL package LINFIT, which takes into account the weighted error associated with each bin. In the case of the cumulative functions, the error analysis has been performed with bootstrapping techniques. To take into account how the photometric uncertainty associated with each point affects the final recovered slope, we perform Monte Carlo realizations of 1000 cumulative distributions. With each observed magnitude, we associate an uncertainty extracted from a Gaussian distribution with standard deviation equal to the maximum tolerated error of 0.3 mag. Each cumulative realization is thus fitted. The final error associated with the observed slope is the standard deviation of all the recovered indexes.

In general, both methods produce slopes that are close to an index of −2. We notice that the recovered slopes for the binned data are shallower than in the case of cumulative distributions. This is mainly the result of the differences between the two techniques (see Section 5.3), hence slopes determined using binned and cumulative distributions cannot be directly compared. Some important features are, however, observable in both analyses. We find a clear steepening as a function of increasing wavelength, i.e., the recovered slopes become significantly steeper than −2 in the BVI filters. This is true for both class 1 and 2 clusters and class 3 associations (see insets in Figure 13). The distributions of the two bluest filters show an extended flat peak that in the cumulative distributions appears as a significant curvature. Moreover, the cumulative distributions at all wavelengths show a clear steepening at the bright magnitude ends. Such steepening is not observed in the binned distributions because while a variable-size binning technique mitigates biases introduced by equally spaced binning approaches (Maíz Apellániz & Úbeda 2005), it also tends to wash out small-scale variations (see Section 5.3). Because of the small number of very luminous objects, the brightest bin of the CLF has a width of about 1.5–2 mag encompassing the range where the steepening in the cumulative function is observed. Our simulations, in Section 5.3, show that the departures from a single power-law function happen mainly at the bright (massive) end, thus the cumulative analysis is better suited to investigate the shape of the luminosity and mass functions.

The luminosity function is a direct observable of the underlying mass function integrated over time. The trends we observe in the CLF of NGC 628 suggest a dearth of luminous (massive) clusters/associations and thus simultaneously analyzing the CLF and CMF can tell us something important about how clusters form and evolve in this galaxy.

To derive masses, we need to perform an extra step where the ages and internal reddening of the sources are extracted. As discussed in K15, the stellar physical properties derived with our deterministic method are severely biased at very low masses due to stochastic variations from the small number of stars. Figure 14 of K15 shows the one-to-one comparison between Yggdrasil deterministic and SLUG stochastically derived cluster properties, suggesting that important deviations occur at cluster masses below 5000 ${M}_{\odot }$. This mass limit was, in recent years, widely adopted in the YSC analysis based on deterministic approaches. We assume the same mass limit in our current analysis. Figure 14 shows the age–mass diagram of the clusters/associations in our two HST pointings of NGC 628. A mass cut of 5000 ${M}_{\odot }$ gives us complete detection up to a stellar age of 200 Myr. Sources falling within the shadowed areas are, thus, not included in the analysis of the CMF and disruption rates presented hereafter. In total, we count 320 class 1 and 2 clusters and 42 class 3 objects that pass this mass and age selection. We notice that our mass cut is very close to the completeness limits of our catalogs at the last age bin at 200 Myr. We tested whether using a higher mass cut at 8000 ${M}_{\odot }$ produces different outcomes. We do not see any change in the recovered CMF properties, only higher errors because of the smaller number of clusters available for the analysis.

The observed mass function of class 1 and 2 (orange triangles) and class 3 (blue dots) systems are shown as cumulative distributions in Figure 15. The observed cumulative distributions are fitted using the IDL maximum-likelihood fitting package MSPECFIT (Rosolowsky 2005). We perform two different fits to the cumulative distributions, a single power-law function, in the integral form $N(M^{\prime} \gt M)\propto {M}^{\alpha }$, and a power-law function with a truncation at the upper mass end, i.e., $N{(M^{\prime} \gt M)\propto {N}_{0}[(M/{M}_{\star })}^{\alpha }-1]$ (see Rosolowsky 2005 and references therein for a complete discussion of the formalism). The resulting fitted parameters for the two functions (single power law in the top panel, truncated power law in the bottom panel) are included in the insets of Figure 15. In Table 2 we list the recovered values for the class 1 and 2 population. M_⋆, the index ${\alpha }_{\mathrm{SF}}$, and N₀, which is the number of objects more massive than ${2}^{1/\alpha }{M}_{\star }$, are determined for a truncated function, while the pure power-law fit provides the index ${\alpha }_{\mathrm{PLF}}$. As described in Rosolowsky (2005), if the resulting N₀ is significantly larger than 1, then a truncated CMF form is preferred to the more traditional single power-law function. When ${N}_{0}\lt 1$, the truncation mass is unconstrained and thus a single power-law fit is sufficient. In Table 2, we also include errors. The errors associated with the observed maximum cluster mass, M_max, and fifth most massive cluster mass, ${M}_{\max }^{5\mathrm{th}}$, have been computed during the SED fitting procedure and described in Section 3. The errors associated with the best-fitting parameters have been computed using deviations from 1000 iterations of bootstrap trials.

Zoom In Zoom Out Reset image size

Table 2.  Parameters Describing the Cluster Mass Function in NGC 628^a

Region	Number	Radius	Mmax	${M}_{\max }^{5\mathrm{th}}$	N0b	${\alpha }_{\mathrm{SF}}$ b,c	${M}_{\star }$ b	${\alpha }_{\mathrm{PLF}}$ c
		(kpc)	(105 ${M}_{\odot }$)	(105 ${M}_{\odot }$)			(105 ${M}_{\odot }$)
bin 1	107	0.46–3.19	${1.62}_{-0.26}^{+0.13}$	${1.11}_{-0.30}^{+0.07}$	2.29 ± 3.40	−1.84 ± 0.13	4.85 ± 2.02	−1.90 ± 0.10
bin 2	107	3.22–4.53	${2.16}_{-0.19}^{+0.18}$	${0.45}_{-0.10}^{+0.11}$	3.77 ± 2.80	−2.15 ± 0.12	0.98 ± 0.42	−2.25 ± 0.12
bin 3	106	4.56–10.15	${4.36}_{-0.27}^{+0.14}$	${0.43}_{-0.11}^{+0.04}$	5.67 ± 4.17	−2.00 ± 0.13	1.04 ± 0.59	−2.13 ± 0.09
all	320	0.46–10.15	${4.36}_{-0.27}^{+0.14}$	${1.40}_{-0.08}^{+0.09}$	7.58 ± 4.20	−2.03 ± 0.07	2.03 ± 0.81	−2.09 ± 0.06

Notes.

^aThe table only includes observed and fitted values of class 1 and 2 cluster population with ages smaller than 200 Myr and masses above 5000 ${M}_{\odot }$. ^bMass function parameter fits, computed via the maximum-likelihood method of Rosolowsky (2005). If ${N}_{0}\gg 1$, a truncated CMF form is appropriate, while ${N}_{0}\lesssim 1$ indicates a single power law is more appropriate. ^c ${\alpha }_{\mathrm{SF}}$ is the slope derived assuming a truncated mass function, and ${\alpha }_{\mathrm{PLF}}$ has been derived assuming a pure power-law function. See the text for details.

Download table as:  ASCIITypeset image

In general, we observe that both a very steep single power-law fit and a truncated function fit with a slightly flatter index can reproduce the observed mass distribution for class 3 objects. However, the number of associations is very small (42), and thus it is not possible to impose any further constraint.

On the other hand, the analysis of the mass distributions of class 1 and 2 systems yields, for both a single power law and a truncated function type fits, slopes very close to −2. However, as already noticed during the analysis of the CLF, the approximation of the CMF by a single power-law function (see the top panel of Figure 15) overestimates the expected number of clusters at the upper mass end of the distribution. A fit to the observed CMF of class 1 and 2 with a truncated power-law function (bottom panel of Figure 15) yields a similar slope, but it mitigates the differences at the high mass end of the CMF distribution. The resulting N₀ (see value listed in Table 2) is larger than 1, suggesting that the latter function provides a better fit to the observed CMF. Thus, a truncated function with slope ${\alpha }_{\mathrm{SF}}=-2.03$ and ${M}_{\star }\sim 2.0\times {10}^{5}\,{M}_{\odot }$ is the statistically favored description of the observed CMF of NGC 628. However, it is important to notice that the number of clusters more massive than $5\times {10}^{4}\,{M}_{\odot }$ is about 22 and only half of those clusters are more massive than 10⁵ ${M}_{\odot }$ so the constraint on ${M}_{\star }$ is weak and the uncertainties on N₀ large.

As an exercise, we try to estimate the expected number of clusters more massive than ${M}_{\star }$. Using the combination of far-UV and 24 μm fluxes of the area covered by the LEGUS pointings of NGC 628, we estimate an SFR of about 0.59 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. Assuming that the SFR was constant for the last 200 Myr, we estimate that a total stellar mass of $1.18\times {10}^{8}\,{M}_{\odot }$ has been formed in the region. Using the cluster formation efficiency definition given in Adamo et al. (2015) and clusters in the age range between 1 and 100 Myr (same as the age range to which the estimated SFR is sensitive to), we derive for this region of NGC 628 a cluster formation efficiency of 12%. This means that 12% of the total stellar mass of $1.18\times {10}^{8}\,{M}_{\odot }$ is in bound clusters, i.e., $1.42\times {10}^{7}\,{M}_{\odot }$. Using the latter amount as the total stellar mass in clusters, we can estimate the number of clusters more massive than ${M}_{\star }$. Observationally, we find two clusters more massive than ${M}_{\star }$. Assuming a pure power-law mass function of slope −2.09 (with upper mass $1.\,\times \,{10}^{7}\,{M}_{\odot }$), we estimate that five clusters more massive than ${M}_{\star }$ should have formed in the last 200 Myr. A Schechter-type function, as described by Equation (3), results in one cluster more massive than ${M}_{\star }$. The estimated total stellar mass in clusters results in cluster numbers that are consistent with the observed ones but does not produce any definitive proof that can help to discern the real shape of the upper mass function. Therefore, the solution provided by a pure power-law function with slope −2.09 cannot be discarded.

We also test whether there is any variation in the CMF properties as a function of galactocentric distance by performing a radial analysis of the CMF for the class 1 and 2 clusters. In total, our sample contains 320 objects more massive than 5000 ${M}_{\odot }$ and younger than 200 Myr. We divide the sample in three radial bins of increasing distance from the center of the galaxy, each containing the same number of objects so that we remove the size-of-sample effect. We then determine the ${M}_{\star }$, N₀, ${\alpha }_{\mathrm{SF}}$, and ${\alpha }_{\mathrm{PLF}}$ of the mass function of each bin using both a truncated and a more traditional power-law function as described above. The recovered values are listed in Table 2. In bin 1, the recovered ${M}_{\star }$ is larger than the most massive cluster observed in the bin, and the resulting N₀, including the uncertainties, is very close to 1. This means that the shape of the upper mass end of the CMF in the inner bin remains unconstrained and the solution with a single power-law fit is as likely. In the second and third bins, between 3 and 10 kpc, ${M}_{\star }$ is consistently 1 × 10⁵ ${M}_{\odot }$ and does not decline significantly. The recovered N₀ is larger than 1 but the uncertainties are large, too. In Figure 16, we plot the derived ${M}_{\star }$ (cyan dots and blue bands, including the uncertainties), the masses and uncertainties of the most massive cluster, M_max (red triangles and orange bands), and the fifth most massive cluster, ${M}_{\max }^{5\mathrm{th}}$ (green triangles and bands), observed in each bin. The derived ${M}_{\star }$ decreases significantly between the innermost and the other two bins. In the two outermost bins, the observed masses of the most massive clusters are significantly more massive than the constrained ${M}_{\star }$, but the numbers of clusters close to the determined ${M}_{\star }$ are small, which makes the statistics (e.g., N₀) quite uncertain.

Zoom In Zoom Out Reset image size

Overall, we observe that if the recovered N₀ is larger than 1, the analyzed cluster population has formed at least a few clusters with masses close to and larger than the truncation value, ${M}_{\star }$. This behavior suggests that ${M}_{\star }$ is not a sharp truncation and that the mass function is likely stochastically sampled. In the literature, it has been suggested that the CMF of YSC populations in local galaxies can be described with a Schechter function of slope close to −2 and a rapid exponential decline above a certain truncation mass (Gieles et al. 2006; Larsen 2009; Bastian et al. 2012). The Schechter function differs from a pure power-law function only at the upper mass end of the distributions and could in principle explain the disagreement between the low numbers of observed massive clusters with respect to the expected one from the extrapolation of a pure power-law function. The probability to form clusters more massive than ${M}_{\star }$ declines exponentially but is not null.

We notice that the ${M}_{\star }$ of class 1 and 2 systems recovered for NGC 628 is very similar to the one retrieved for M83 (Adamo et al. 2015), in agreement with the evidence presented in Larsen (2009), who suggested that ${M}_{\star }$ in local spiral galaxies is about a few times 10⁵ ${M}_{\odot }$.

In the attempt to understand the uncertainties imposed by the low number statistics and, at the same time, link the observed CMF to the CLF we use simulated cluster populations. We stochastically sample a pure power-law mass function (we will refer to this population as run A) and a truncated one in the form of a Schechter function (run B), with the same slopes and ${M}_{\star }$ derived for the cluster (class 1 and 2) population of NGC 628. In both runs, we sample the mass function from 2 × 10² to 100 × M_max ${M}_{\odot }$. We use $1.5\times {10}^{4}$ objects so that the resulting cluster population will approximately have the same number of clusters more massive than 5000 ${M}_{\odot }$ as observed in NGC 628 (i.e., 320 class 1 and 2 objects). For each cluster, we stochastically assign an age between 1 and 200 Myr, i.e., we assume that star formation has been constant during this time range. Based on the evidence produced by the analysis of cluster dissolution timescales in Section 5.4, we also include mass-dependent disruption using the formula:

$$\begin{eqnarray}&&{M}_{i}({t}_{i})={({M}_{i}^{\gamma }-\gamma {t}_{i}/{t}_{0})}^{-\gamma },\end{eqnarray} \tag{ 1 }$$

where we assume $\gamma =0.65$ and t₀ = 4.91 × 10⁵ years (values obtained from a maximum-likelihood fit to the data; see Section 5.4 for details). We then create 1000 realizations of each population.

In Figure 17, we include the observed CMF of NGC 628 (orange dots) and we overplot the median (red solid line) mass cumulative distribution and the distributions containing 50% (dashed red line) and 90% (dotted red line) of the 1000 Monte Carlo realizations sampled from an underlying pure power-law mass function (run A, top panel) and Schechter truncated function (run B, bottom panel). To build the cumulative distributions of the simulated populations we apply the same mass cut as in the observations at 5000 ${M}_{\odot }$. When comparing the loci occupied by the simulations with respect to the observed CMF we see that 90% of realizations of run A overestimate the number of clusters more massive than 10⁵ ${M}_{\odot }$ (the locations of the upper mass end of the observed CMF is off the 90% limits of realizations). We test the null hypothesis that the upper part of the median distribution of run A and run B realizations and of the observed CMF are drawn from the same parent population. Since the differences between the pure power-law function and the truncated function is at the upper mass end of the distributions, we run both a Kolmogorov–Smirnov (KS) and an Anderson–Darling (AD) test using clusters more massive than 10⁴ ${M}_{\odot }$. We recover similar probabilities from both statistics with p(AD) and p(KS) of ∼0.3 and ∼0.9 for run A and run B, respectively. This result does not discard any of the two functions, but it yields a marginal preference for a truncated mass function.

Zoom In Zoom Out Reset image size

Next, we analyze the resulting luminosity functions produced from the two run A and run B simulated cluster populations with the goal of understanding if the underlying CMF has a truncation. We use the median cluster population of the 1000 Monte Carlo realizations to yield the CLF. We do not apply any mass cut to the simulated populations, but we use the age and mass of each mock object together with Yggdrasil models to estimate their fluxes in the UBVI LEGUS bands. We then select only sources with a luminosity brighter than ${M}_{V}\leqslant -6$ mag and brighter than the 90% completeness values in the other two bands for which we require detection for a source to enter the catalog (namely B and I). Extinction is not taken into account.

The resulting cumulative and binned CLFs as a function of waveband of the two simulated cluster populations are illustrated in Figure 18. As done in Section 5.1 for the observed CLFs, we fit a single power-law function to the simulated CLFs. To recreate the limiting detection depth of each waveband, we use the same low luminosity limits as for the observed distributions.

Zoom In Zoom Out Reset image size

Since the underlying CMF of the simulated cluster populations is known, we can directly verify the effect of using either binning or cumulative techniques on the resulting distributions. We notice that the recovered slopes ${\beta }_{\mathrm{PLF}}^{\lambda }$ and ${\beta }_{\mathrm{SF}}^{\lambda }$ listed in the left panel of Figure 18 are about 0.2 steeper than the ones listed in the right plot of the same figure. As discussed in Section 5.2, the way the equal-number object binning and cumulative distributions are built make them sensitive to different properties of the distributions. Since we are interested in understanding whether the CMF has a truncation at high mass, we prefer to analyze cumulative distributions.

The cumulative distributions in the left panel of Figure 18 should be compared to the observed cumulative distributions of class 1 and 2 objects (cumulative functions illustrated with triangles in the left plot in Figure 13). First, the resulting indexes of run A (${\beta }_{\mathrm{PLF}}^{\lambda }$) and run B (${\beta }_{\mathrm{SF}}^{\lambda }$) show that the effect of fading when comparing the slope of the U band to the ones in the redder filters is between 0.1 and 0.15. The differences in the recovered slopes for the observed CLFs in the ultraviolet and optical are larger (about 0.2–0.3). However, we note that our simulations do not include extinction, which preferentially affects the bluer bands and could explain the larger flattening of the observed ultraviolet slopes.

Second, we notice that the recovered CLF slopes for a population of clusters sampled out of a pure power-law CMF (run A) are very close to the initial slope of the corresponding CMF (i.e., ${\alpha }_{\mathrm{PLF}}=2.09$). The slopes obtained fitting a CLF sampled out of a CMF with a power-law shape (slope −2) and a truncation at the high mass end (run B) are steeper than −2. We also notice a fundamental difference at the bright luminosity end between the cumulative distributions of the two runs (pure power law versus power law with a truncation). The presence of a truncation at the high mass end of the CMF (run B) results in significant deviations at the bright end (i.e., solid triangles) of the CLF, with respect to the population sampled out from a pure power-law CMF (run A, open circles). The drop at the bright luminosity observed for the run B CLF (filled circles, left panel of Figure 18) is very similar to the ones observed in the real CLF (class 1 and 2, left plot in Figure 13). This trend reinforces the evidence of a presence of a truncation in the underlying CMF of NGC 628.

In this section, we probe the evolution of the two classes of objects visually identified: class 1 and 2 likely being clusters, and class 3 being a heterogeneous sample dominated by compact stellar associations.

We apply the same mass cut at 5000 ${M}_{\odot }$, which ensures that we will be complete in their recovery rate up to 200 Myr, but to look now at the number densities of objects per unit time (dN/dt) as a function of increasing age. A decreasing rate of objects as a function of time ${dN}/{dt}\propto {t}^{\delta }$ is historically interpreted as evidence for YSC dissolution within the galaxy (see Lamers 2009 for a short review). In Figure 19, top panel, we plot the change in number densities of class 1 and 2 (orange dots), class 3 (blue diamonds), and the whole sample (green triangles) using bins of the same width (left panel). The shadowed areas are the regions of the diagram where incompleteness mimics a rapid disappearance of objects. Therefore, these regions are excluded (e.g., see Lamers 2009). Overall, we observe that the disruption rate of class 3 systems is certainly more significant than class 1 and 2. When the two classes are analyzed together, the resulting disruption rate is coincident with the one of class 1 and 2, probably because this latter class is much more numerous.

Zoom In Zoom Out Reset image size

5.4.1. Cluster and Compact Association Evolution During the First 10 Myr

5.4.2. The Change in the Number Density of Clusters and Associations as a Function of Time, Mass, and Galactocentric Distance

Between 10 and 200 Myr, the left plot of Figure 19 shows that the change in the number density of class 1 and 2 clusters is consistent with mild disruption. We do not observe any drastic dissolution of this class of objects up to 200 Myr. On the other hand, the number density of class 3 systems keeps declining by a factor of 3 in the same age range, suggesting, as already observed in the color–color diagrams analysis, their rapid disappearance within 100 Myr (the slope we recover is $\delta \sim -0.7$). These short timescales are in agreement with the clustering analysis performed by Grasha et al. (2015) on the class 3 population. The lack of clustering after 40 Myr may be caused by the quick dissolution of these systems. On the other hand, since class 1 and 2 survive longer, their lack of spatial clustering is possibly the result of the randomization of their positions because they move away from their birthplaces. The clustering analysis of stellar structures performed by Gouliermis et al. (2015) in another LEGUS galaxy, NGC 6503, finds that hierarchical clustered stellar structure disappear and distribute into the stellar field within 60 Myr. This suggests that class 3 objects are well nested within the hierarchical properties of star formation, while stellar clusters, even though have formed within the same hierarchically structured ISM, as gravitationally bound systems, may follow a different fate.

Recent theoretical and numerical works (e.g., Elmegreen & Hunter 2010; Kruijssen et al. 2011) point out that tidal forces exerted by GMC encounters in a hierarchical ISM can reproduce the steady decrease in cluster numbers over time as a result of a decreasing ISM density in each cluster's environment as it drifts away from its birth place. Higher gas densities toward the centers of spiral galaxies or in starburst systems can increase the overall dissolution rate of YSCs.

In Figure 8, the color–color diagrams of the inner and outer cluster population of NGC 628 show interesting features, suggesting that the outer cluster population is preferentially older. We attempt here to investigate whether we observe any change in the recovered dissolution time as a function of galactocentric distances. We split our sample of class 1 and 2 objects into a central and outer bin containing the same number of objects. In the central panel of Figure 19, we show the recovered disruption rates for the inner and outer bins. We observe that the number density of clusters in the outer bins is consistent with being constant between 10 and 200 Myr ($\delta \sim 0.0$), while in the inner bin we observe higher dissolution rates $\delta \sim 0.3$. The observed trends are consistent with theoretical expectations, suggesting a higher mass-independent disruption in denser galactic ISM, thus closer to the center of the spiral galaxy.

In the right panel of Figure 19, we probe the disruption rate of low (magenta dots) and high (green diamonds) mass cluster candidates. We observe that low mass clusters show evidence of mild disruption ($\delta \sim -0.3$) while the number densities of clusters with mass similar to or larger than 10⁴ ${M}_{\odot }$ are consistent with being constant ($\delta \sim 0.0$). This result suggests the presence of a mass dependency in the disruption rate of clusters.

Another way to investigate cluster evolution as a function of the cluster mass is to encode in a bivariate distribution, $g(M,t)$ the mass function and time dependence for formation and evolution of clusters (see Fall et al. 2009; Gieles 2009). The function $g(M,t)$ is thus the number of clusters observed as a function of time and mass expressed as

$$\begin{eqnarray}&&g(M,t)=\displaystyle \frac{{d}^{2}N}{{dMdt}}.\end{eqnarray} \tag{ 2 }$$

Integrating $g(M,t)$ over the mass provides the dN/dt distributions, while integrating $g(M,t)$ over a time range provides the observed CMF.

In Figure 20 we show the recovered mass distributions as a function of cluster age normalized by the corresponding age interval (i.e., dN/dMdt diagnostic). The distribution is built to contain the same number of objects in each bin. The dN/dMdt diagnostic is sensitive to the evolution of the mass function with time. If the number of clusters rapidly declines independently of the cluster mass then the shape of the mass function will be unchanged but the distribution at each age interval will be shifted because the number of clusters is diminishing. On the other hand, if the cluster disruption timescale depends on cluster mass, then the aging CMFs should overlap at the mass ranges untouched by disruption and deviate where disruption of low mass clusters is significant. The plot in Figure 20 shows a clear offset of the CMF for clusters younger than 10 Myr. As already seen in the dN/dt distribution, the number of clusters at this age range is significantly higher. The dN/dMdt analysis provides further insights showing that the number of clusters younger than 10 Myr is higher at any mass range in units of time. The CMFs of the other two age ranges (10–100 and 100–200 Myr) overlap at masses larger than 10⁴ ${M}_{\odot }$, while deviations become significant at lower masses. These deviations suggest that the number of clusters with masses below 10⁴ ${M}_{\odot }$ are decreasing, flattening the CMF. This trend is consistent with both the effect of detection limits causing the loss of an increasing number of low mass clusters at older ages and higher disruption rates of clusters below 10⁴ ${M}_{\odot }$ (consistent with the trends observed in the right plot of Figure 19). As pointed out in Section 5.2, our mass cut of 5000 ${M}_{\odot }$ is very close to the completeness limits at 200 Myr (see Figure 14). Incompleteness can mimic a mass-dependent cluster dissolution, and we cannot exclude that the observed flattening is in part caused by incompleteness.

Zoom In Zoom Out Reset image size

To understand whether the trends observed in the bivariate distribution of the CMF as a function of time are compatible with mass-dependent disruption we perform a maximum-likelihood fit to our data. We refer to Gieles (2009) for the formalism behind the fit and to Bastian et al. (2012) for an application to the YSC population in the M83 galaxy. The fit is performed assuming a Schechter function that describes the probability to form a cluster of mass M_i in the time interval t and $t+{dt}$ as

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}N}{{{dM}}_{i}{dt}}=A\,{M}_{i}^{-2}\,\exp \left(-\displaystyle \frac{{M}_{i}}{{M}_{\star }}\right),\end{eqnarray} \tag{ 3 }$$

where M_i depends on time. If only cluster mass-dependent disruption is taken into account, the disruption time can be described as ${t}_{\mathrm{dis}}={t}_{0}{M}^{\gamma }$, with t₀ and γ depending on the galactic environment (see e.g., Lamers et al. 2005), and the mass evolution is described by $M={({M}_{i}^{\gamma }-\gamma t/{t}_{0})}^{1/\gamma }$ (Lamers et al. 2005; Gieles 2009). In the fit, γ is fixed to the average value found in local spirals ($\gamma =0.65$). The fitting algorithm thus finds the values of ${M}_{\star }$, t₀, and t₄ that maximize the likelihood. The t₄ is the timescale necessary to dissolve a cluster of ${M}_{4}={10}^{4}\,{M}_{\odot }$, and it depends on t₀ as ${t}_{4}=({t}_{0}/\gamma ){M}_{4}^{\gamma }$. The results of maximum-likelihood fits performed on the YSC population of NGC 628 are shown in Figure 21. The fit is performed on class 1 and 2 clusters more massive than 5000 M_⊙, brighter than the magnitude cut ${M}_{V}=-6$ (V = 23.98 mag) mag, and younger than 300 Myr. The age–mass diagram on the left side panels shows the YSC population taken into account in the fit. The mass cut removes low mass objects up to an age of 200 Myr. The magnitude cut becomes important at older ages because it removes clusters more massive than 5000 ${M}_{\odot }$ at ages older than 200 Myr. The fit is done for the entire cluster population (top), and the inner (center) and outer (bottom) regions of the galaxy corresponding to two bins containing the same number of clusters. The maximum-likelihood fit performed on the entire cluster population finds ${M}_{\star }=2.2\,\times {10}^{5}\,{M}_{\odot }$ and a timescale for dissolution of a 10⁴ ${M}_{\odot }$ cluster of ${t}_{4}=190\,\mathrm{Myr}$. These values are in agreement with the ${M}_{\star }$ found in Section 5.2 and the evidence of slow disruption in the dN/dt analysis performed above. With the same technique, we derive ${M}_{\star }$ and t₄ for clusters of the inner and outer regions containing the same number of objects, as defined above. The cluster population located within the inner bin has the same ${M}_{\star }$ value found for the entire galaxy (${M}_{\star }=2.2\times {10}^{5}\,{M}_{\odot }$) and a timescale for dissolution of a 10⁴ ${M}_{\odot }$ cluster of ${t}_{4}=130\,\mathrm{Myr}$. The clusters in the outer region have a factor of 2 longer disruption timescales with ${t}_{4}=270\,\mathrm{Myr}$ and ${M}_{\star }=2.0\times {10}^{5}\,{M}_{\odot }$. The maximum-likelihood fitting analysis confirms the results and trends obtained from the analyses of the CMF and the dN/dt distributions. The average timescale for cluster disruption in the galaxy are long enough to produce a shallow dN/dt distribution. However, the age interval to which our analysis is sensitive too is long enough that we should start to see the effect of disruption on low mass clusters. Indeed, in this respect the dN/dMdt analysis is a more sensitive diagnostic than the simple dN/dt one. The longer timescales for the disruption of clusters in the outer part of the galaxy are also in agreement with the expectation from theoretical works (e.g., Elmegreen & Hunter 2010; Kruijssen et al. 2011) and the results obtained for the M83 spiral galaxy (Bastian et al. 2012).

Zoom In Zoom Out Reset image size

To verify the effect of incompleteness at the low mass end as a function of age, we repeat the maximum-likelihood fit using the same mass cut at 5000 ${M}_{\odot }$ together with a more conservative luminous cut, i.e., M_v = 22.93 mag, i.e., one magnitude brighter than previously done. The mass and magnitude cuts result in a selection of clusters more massive than 5000 ${M}_{\odot }$ up to 100 Myr and objects brighter than 22.93 mag at older ages. We obtain values of ${M}_{\star }$ and t⁴ within a factor of two from the analysis performed above with a less conservative magnitude cut. With a limiting brightness of ${M}_{v}=-22.93$ mag, we retrieve ${M}_{\star }\,=1.7\times {10}^{5}\,{M}_{\odot }$ and ${t}^{4}=170\,\mathrm{Myr}$ for the entire sample, and ${M}_{\star }=2.2\times {10}^{5}$ and $1.4\times {10}^{5}\,{M}_{\odot }$, and ${t}^{4}=100$ and 530 Myr for the central and out region of NGC 628, respectively. This further test reinforces the evidence, found above, of a mass-dependent component necessary to describe cluster disruption in NGC 628.

CLF and CMF are powerful tools to investigate the formation of YSCs.

In the spiral galaxy NGC 628, using the LEGUS data set, we recover the properties of the CLF from the UV to the NIR. The luminosity distributions of the cluster (class 1 and 2) population is close to the slope −2, while the compact stellar associations (class 3) appear to have steeper slopes. The analysis of increasing size of stellar aggregates (Elmegreen et al. 2006) also finds slopes close to −2. Hence, star formation is consistent with a hierarchical process, driven by turbulence, from the largest scales, i.e., star-forming complexes, down to the densest and smallest physical scales, i.e., star clusters. Our analysis also reveals some interesting variations in the CLF for both clusters and associations. At each waveband, we observe a steepening at the brightest end of the cumulative distributions. Extinction cannot explain the steepening at the brightest end of the CLF. If more luminous clusters are more extinguished than faint ones, it would imply preferential extinction as a function of luminosity, a trend not observed in studies of local galaxies.

The cause of the steepening of the CLF could be connected to the nature of the cluster formation process. The CLF is a direct observable tracer of the underlying CMF integrated over time. If the CMF is a pure power-law function of slope −2 at any age and cluster dissolution is independent of the cluster mass, the resulting CLF should consistently be a function with the same slope. The steepening that we are observing at the bright end of the distributions suggests that we find fewer luminous (massive) clusters than expected for a pure power-law CMF.

Increasing evidence in the literature suggests that the YSC mass function is better described by a truncated power law. Johnson et al. (2017), studying the YSC population of M31, find that the CMF in this galaxy has a Schechter-type form, with ${M}_{\star }$ ∼1 × 10⁴ ${M}_{\odot }$. Cosmological simulations of the Milky Way type of galaxies, where YSCs are implemented as star formation units, show that the resulting YSC populations have CMF that are better described by a Schechter function and ${M}_{\star }$ scales as a function of SFR (Li et al. 2017).

In the case of NGC 628, both a steeper power law ($\alpha \sim 2.1$) or a Schechter function ($\alpha \sim 2.0$ and ${M}_{\star }=2.0\times {10}^{5}\,{M}_{\odot }$) can equally reproduce the observed CMF. We investigate whether the steepening observed in the mass function, which suggests a dearth of very massive systems (if a pure power-law function is assumed), is a sign of an underlying soft truncation. The difficulty in constraining the upper mass distribution of the CMF in local spirals has already been pointed out by Larsen (2006). The combination of cluster formation being a stochastic process (e.g., Adamo & Bastian 2015) and relatively low SFR combined with detection limits and rapidly fading luminosities above a few hundreds of Myr make it very challenging to put strong constraints on the shape of the cluster upper mass distribution, because of low number statistics.

Another fundamental aspect that gives important constraints on the formation and evolution of YSCs is understanding how they disrupt. Star-forming regions and stellar complexes are hierarchical in space and time (Efremov & Elmegreen 1998), thus their crossing times are comparable to their ages and appear to dissolve on timescales below 60 Myr (e.g., Pellerin et al. 2012; Crocker et al. 2015; Gouliermis et al. 2015). On the other hand, the fate of YSCs, forming in the densest peaks of these very regions, is not yet well understood.

In our analysis, we look at the number density of stellar systems between 1 and ∼200 Myr. We are not able to derive the dynamical status of the stars within our objects, but since the average size of our class 1 and 2 systems peaks at 3 pc (Ryon et al. 2017), their crossing times, for those older than 10 Myr, are much shorter than the stellar ages so we can consider them likely gravitationally bound. We see a clear decline between the number of objects during the first 10 Myr and the next age bin, suggesting a significant loss of both class 1 and 2 and class 3 systems. The recovered number densities of YSCs versus associations in the age range 10–200 Myr are significantly different. Compact stellar associations (class 3) rapidly decline in number and disappear on timescales (∼50 Myr) comparable to those of hierarchically structured star-forming regions. Objects that we have identified as potentially bound clusters (class 1 and 2) show close to constant number densities ($\delta \sim -0.2$) in the dN/dt analysis.

Small disruption rates within the first 100–200 Myr have also been reported in the literature for other local galaxies, e.g., in M31 (Fouesneau et al. 2014), LMC (Baumgardt et al. 2013), and M83 (Silva-Villa et al. 2014). However, we notice here that these results are not ubiquitous. A compilation of previous works (Adamo & Bastian 2015) clearly shows that disruption rates of YSCs may change significantly from galaxy to galaxy and become very high for YSCs in hostile environments like the Antennae system (Whitmore et al. 2010).

Our results suggest a steeper disruption slope for YSCs located in the inner portion of the galaxy ($\delta \sim -0.3$). This finding is observationally supported by the difference in the number densities of YSCs in regions of the color–color diagrams (Figure 8) corresponding to young ages. In the inner pointing, the number of YSCs is larger than in the outer pointing at young ages, while similar number densities are observed at older ages between the populations of the two pointings. If star formation has been constant in the last few hundreds of Myr, then these differences suggest a longer survival time in the outer regions. A larger coverage of the galaxy would certainly improve the results, as recently shown for M83 (e.g., Chandar et al. 2010, 2014; Bastian et al. 2012; Silva-Villa et al. 2014).