The Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER). III. The Mass Function of Young Stellar Clusters in M33

For each cluster in our sample, we extract an optical CMD in the F475W and F814W filters. These passbands yield the deepest CMDs available from the PHATTER data. Specifically, we use CMDs composed of stars that lie within the photometric aperture (R_ap) derived in Johnson et al. (2022), which corresponds to approximately three times the cluster half-light radius. We assume all stars within R_ap are cluster members, and all members are within R_ap. We make no correction for mass that lies outside R_ap. Based on experiments with synthetic clusters, we expect <30% of light to fall outside our photometric aperture. The median aperture correction for clusters in our sample is −0.04 mags, with the largest correction being −0.69. This suggests we may be losing ∼4% of the light from our clusters in typical cases, and thus underestimating clusters logM by 0.02 dex. This number is much smaller than the uncertainties in our M_c values derived below, and we don't correct our masses for this aperture correction. We characterize the surrounding field star population using an annulus that spans 1.2–3.4 R_ap. This annulus has 10× the area of the cluster aperture. The background is fit along with the cluster models; the scaling of the background is a free parameter in the fit.

We use the MATCH software package to perform maximum-likelihood CMD fits and derive constraints on cluster properties following techniques described in Dolphin (2002). Unlike its typical use in determining time-resolved star formation histories (SFHs), we use MATCH in a more limited simple stellar population mode. Here, the model CMD is composed of a population drawn from a single time bin, rather than a linear combination of populations from multiple time bins.

The MATCH code uses theoretical isochrones to populate synthetic CMDs according to input parameters of age, dust extinction (A_V ), distance, metallicity, stellar initial mass function (IMF), and binary fraction. Synthetic populations are created from unique combinations of input parameters, convolved with a model of observational completeness and noise derived from artificial star tests (ASTs), to produce a simulated CMD. This simulated CMD is combined with a background model created from the CMD of stars lying in an annulus surrounding the cluster, then both the model and background components are scaled to best reproduce the observed cluster CMD. The fit quality is calculated according to a Poisson likelihood, and the code iterates through a grid of input parameter value combinations to map the distribution of probability.

We adopt the M33 distance modulus of 24.67 (de Grijs & Bono 2014), and assume a metallicity [M/H] of −0.15 ± 0.25 based on the range of present-day gas phase metallicity in M33 (e.g., U et al. 2009). For young clusters, the age is heavily weighted by the main sequence and thus has little metallicity sensitivity. For consistency with previous studies, our assumptions in stellar modeling follow the cluster fitting of Weisz et al. (2015) and Johnson et al. (2016). Briefly, we adopt a binary fraction of 0.35 with a uniform mass ratio distribution, a Kroupa (2001) stellar IMF between 0.15 and 120 M_⊙, and Padova stellar models (Marigo et al. 2008) with low-mass asymptotic giant branch tracks from Girardi et al. (2010).

For each cluster, we perform 25,000 ASTs to ensure accurate characterization of photometric completeness and noise, encompassing a wide range of CMD positions and cluster radii. Input positions for ASTs are distributed based on the measured half-light radii of the clusters, and assume a King (1962) profile with a concentration of 10.

We compute CMD fits for a grid of age (6.6 < log(Age/yr) < 9.0) and dust extinction (0 mag < A_V < 2 mag), deriving the mass at each grid point from the best-fit CMD model. We use relative likelihoods derived across the grid to obtain marginalized probability distribution functions (PDFs) for each parameter. We adopt the best-fit model for mass, age, and A_V , and assign each an uncertainty defined by the 16th and 84th percentiles of the PDF. The masses derived represent the initial masses for the clusters.

Figure 1 shows the full sample of best-fit ages and masses for the M33 cluster sample. We find the majority of the cluster age and mass PDFs to be Gaussian. In <10% of cases, the PDFs are bimodal or have a tail to lower or higher values. The median 1σ uncertainty in age is 0.11 dex, while for mass, the median is 0.04 dex. The full catalog of cluster parameters is presented in Table 1. For the remainder of the paper, we use the maximum a posteriori values for the cluster ages and masses; these are contained in the best-fit column in Table 1.

Zoom In Zoom Out Reset image size

While our results focus only on clusters <300 Myr, we note that there also exists of a relatively large number of massive ∼1 Gyr old clusters. A similar sample of clusters is not present in the M31 PHAT data (Johnson et al. 2016).

After fitting the cluster CMDs, we perform a visual inspection of each cluster's CMD fits and optical images. From this inspection, we exclude results for 33 clusters that a group of coauthors unanimously agreed were poor fits. These clusters are older clusters with few detected member stars that were poorly and erroneously fit with a young, high A_V model. These clusters are denoted by an "exclude" flag in Table 1, and are represented by the red open circles in Figure 1. We note that excluding these flagged fits does not significantly impact the CMF results (i.e., differences in parameter fits with and without these clusters are much less than 1σ). We also visually identify 13 globular cluster candidates (blue open circles in Figure 1) and exclude these objects from our sample; globular clusters are known failure cases for our CMD analysis due to limits in the age and metallicity range assumed for the fitting. Most (11/13) of these candidates appear in previous catalogs (San Roman et al. 2010), and every cluster in this subsample with an existing age estimate from Fan & Grijs (2014) was reported as >1 Gyr old.

We select a sample of clusters for mass function analysis in the age range $7.0\lt \mathrm{log}(\mathrm{Age}/\mathrm{yr})\lt 8.5$. The lower threshold is adopted because during the first 10 Myr of a cluster's life, clusters are embedded, making optical observations difficult, and our sample incomplete. In addition, young embedded groupings of stars (<10 Myr) are still forming through hierarchical merging of subclumps, making it unclear whether the resulting structure will be a long-lived, gravitationally bound cluster (Allison et al. 2010; Gieles et al. 2012; Messa et al. 2021). The upper threshold is where CMD fitting becomes less accurate due to the cluster's main-sequence turnoff dropping below the 50% detection limit for the stellar photometry.

In Figure 2, we present CMD fitting results for two example clusters that lie at the median age (∼100 Myr; LGCS-M33 110) and maximum age (∼300 Myr; LGCS-M33 156) of the cluster sample. For each cluster, we show the observed CMD, modeled CMD, residuals, and the significance of the residuals.

Zoom In Zoom Out Reset image size

We find that the uncertainty of our age estimates increases as cluster masses decrease, and that large age uncertainties also translate to less reliable mass estimates. As a result, we adopt a minimum mass of 1000 M_⊙ for our CMF fitting cluster sample. Below this mass limit, the average 1σ error in log(Age/yr) is 0.30 dex for cluster masses from $2.8\leqslant \mathrm{log}(M/{M}_{\odot })\lt 3.0$, whereas the average 1σ error improves to 0.21 dex for clusters with masses from $3.0\leqslant \mathrm{log}(M/{M}_{\odot })\lt 3.2$. We note that this selection of minimum mass has a less than 1σ variance in the inferred CMF.

Our final cluster sample for CMF fitting includes 711 clusters. The selected age and mass range is denoted by the green dashed lines in Figure 1, and the spatial distribution of the selected sample is shown in Figure 3. For this selected sample of clusters, the median 1σ age uncertainty is 0.10 dex, and the median 1σ mass uncertainty is 0.03 dex. The distribution of cluster mass errors is discussed further in Section 5.

Zoom In Zoom Out Reset image size

Once we obtained measurements of cluster ages and masses, we needed to correct the cluster catalog for completeness to properly fit the intrinsic mass distribution. The full description of the cluster sample completeness can be found in Section 4 of Johnson et al. (2022). Briefly, the completeness of the cluster sample was determined by measuring the detections of synthetic clusters placed in LGCS images (Section 2.4 of Johnson et al. 2022). The synthetic clusters are characterized by log(Age/yr) versus log(Mass/M_⊙), binned as a function of log(Age/yr). The completeness as function of log(M/M_⊙), C, is characterized using the functional form of a logistic function given by:

$$\begin{eqnarray}&&\begin{array}{l}C({M}_{50},\mathrm{log}(M/{M}_{\odot }))\\ =\,{\left(1+\exp [-k(\mathrm{log}(M/{M}_{\odot })-{M}_{50})]\right)}^{-1}\end{array}\end{eqnarray} \tag{ 1 }$$

where k sets the slope of the logistic function (fixed to 6.02 as explained in Johnson et al. 2022) and M₅₀ is the 50% completeness limit. M₅₀ is well described by an exponential function:

$$\begin{eqnarray}&&{M}_{50}(\tau )=a\times \exp ((b\times (\tau -{\tau }_{\min })))+c\end{eqnarray} \tag{ 2 }$$

where $\tau \equiv \mathrm{log}(\mathrm{Age}/\mathrm{yr})$, and ${\tau }_{\min }$ is the median cluster log(Age/yr) in the youngest bin, which is 7.09. The constants a, b, and c were fit through minimizing the χ² on the binned data set of synthetic cluster detections over the age range of $7.0\lt \mathrm{log}(\mathrm{Age}/\mathrm{yr})\lt 8.5$. The best-fit M₅₀(τ) parameters are a = 0.0303, b = 1.9899, and c = 2.9770, with a reduced χ² of 0.82.

Johnson et al. (2022) found that environment plays a significant role in completeness, largely due to the density of main-sequence stars in search fields that contain a cluster. We characterize the number of main-sequence stars per LGCS search image (∼36'' × 25'') as log(N_MS), which ranges from $2.0\lt \mathrm{log}({N}_{\mathrm{MS}})\lt 3.75$ in M33. Johnson et al. (2022) found that over this range, the 50% mass completeness is impacted by up to 0.6 dex.

We further examine the environmental completeness dependence here and incorporate it into our completeness model, by deriving a relationship between the 50% mass completeness and log(N_MS). We first exclude synthetic clusters that are greater than 2σ outliers in log(N_MS), then split the cluster into three equal bins of $\mathrm{log}({N}_{\mathrm{MS}})$, and fit the mass completeness as a function of log(Age/yr) for each bin. We find a linear trend between the average 50% mass completeness values for median values of the three $\mathrm{log}({N}_{\mathrm{MS}})$ bins, which we characterize using a linear function $m\times (\mathrm{log}({N}_{\mathrm{MS}}))+{b}_{\mathrm{nms}}$, where m = 0.7727 and b_nms = 0.6674. We adopt the value at the bin edge for $\mathrm{log}({N}_{\mathrm{MS}})$ values lower than the lowest bin and higher than the edge of the highest bin to avoid extrapolating where our number statistics are minimal. The 50% mass completeness as a function of log(N_MS) is shown in Figure 4.

Zoom In Zoom Out Reset image size

We incorporate the environmental impact on completeness from Figure 4 into the sample completeness function by having the c parameter in Equation (2) be dependent on log(N_MS), such that our completeness function becomes,

$$\begin{eqnarray}&&{M}_{50}(\tau ,{N}_{\mathrm{MS}})=a\times \exp (b(\tau -{\tau }_{\min }))+c({N}_{\mathrm{MS}})\end{eqnarray} \tag{ 3 }$$

and c(N_MS) is described by,

$$\begin{eqnarray}&&\begin{array}{l}c({N}_{\mathrm{MS}})\\ \ =\,\left\{\begin{array}{l}2.77\mathrm{log}({N}_{\mathrm{MS}})\lt 2.53\\ (m\times \mathrm{log}({N}_{\mathrm{MS}})+{b}_{\mathrm{nms}})2.53\leqslant \mathrm{log}({N}_{\mathrm{MS}})\leqslant 3.49\\ 3.47\mathrm{log}({N}_{\mathrm{MS}})\gt 3.49\end{array}\right..\end{array}\end{eqnarray} \tag{ 4 }$$

In Figure 5, we show the completeness correction with respect to the raw data. We note that this figure is a visual representation of binned data, and the completeness correction is done on a cluster-by-cluster basis.

Zoom In Zoom Out Reset image size

The likelihood function for an observed cluster with mass M is given as

$$\begin{eqnarray}&&{p}_{\mathrm{cluster}}(M| \theta ,\tau )\equiv \displaystyle \frac{1}{Z}{p}_{\mathrm{MF}}(M| \theta ){p}_{\mathrm{obs}}(M| \tau ),\end{eqnarray} \tag{ 5 }$$

where the cluster distribution is represented by p_MF(M∣θ) defined by parameter θ, and p_obs(M∣τ) is the observational completeness given as a function of cluster age τ. In order to have the likelihood integrate to 1, the normalization Z is required and given by

$$\begin{eqnarray}&&Z=\int {p}_{\mathrm{MF}}(M| \theta ){p}_{\mathrm{obs}}(M| \tau ,{N}_{\mathrm{MS}}){dM}.\end{eqnarray} \tag{ 6 }$$

We note that this normalization is calculated for each cluster, and does not refer to Bayesian evidence.

We test a pure power-law distribution where the only parameter is given by θ = (α); α is defined as the power-law index. We also adopt a Schechter (1976) functional form that has two parameters given by θ = (α, M_c ) where α is the power-law index for low-mass clusters, and M_c is the characteristic mass defining the exponential high-mass truncation. The Schechter (1976) distribution follows the form:

$$\begin{eqnarray}&&{p}_{\mathrm{MF}}{(M| \alpha ,{M}_{c})\propto (M/{M}_{c})}^{\alpha }\exp (-M/{M}_{c}),\end{eqnarray} \tag{ 7 }$$

as first used in the cluster context by Larsen (2009).

To limit the impact of clusters with low completeness, we cut out any clusters below their local 50% completeness limit, M₅₀(τ, N_MS), as expressed in Equation (3). Mathematically, this cluster selection can be described as:

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{obs}}(M| \tau ,{N}_{\mathrm{MS}})\\ \ =\,\left\{\begin{array}{l}{\left(1+\exp \left[\displaystyle \frac{-k(M-{M}_{50}(\tau ,{N}_{\mathrm{MS}}))}{{M}_{\odot }}\right]\right)}^{-1},\\ \ M\gt {M}_{50}(\tau ,{N}_{\mathrm{MS}})\\ 0,\\ \ \ \ \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 8 }$$

We find 533 clusters have masses larger than the local 50% completeness limit.

We use Bayes' theorem to derive the posterior probability distribution function for the mass function parameters given as

$$\begin{eqnarray}&&p(\theta | {M}_{i},\tau )\propto {p}_{\mathrm{cluster}}(\theta | {M}_{i},\tau )p(\theta ),\end{eqnarray} \tag{ 9 }$$

where p(θ) is the prior probability of the Schechter parameters, M_i is a set of 533 clusters, and p_cluster(θ∣M_i , τ) is the combined likelihood function for the full set of clusters defined below. We adopt a uniform tophat prior for the Schechter parameters that covers the range of published values with plenty of cushion: −3 ≤ α ≤ −1 and $3\leqslant \mathrm{log}({M}_{c}/{M}_{\odot })\leqslant 8$.

The assumption with this implementation of a probabilistic approach is that the cluster mass uncertainties are negligible. Weisz et al. (2013) demonstrated that this assumption is valid if the fractional mass uncertainties are smaller than 10%. However, when the fractional error of the masses extends to 50% and beyond, fitting results can become affected (Johnson et al. 2017). We will further discuss this assumption in Section 5.

The likelihood function for a set of clusters is defined as the product of the individual cluster mass probabilities given by,

$$\begin{eqnarray}&&p({M}_{i}| \theta ,\tau )\equiv \prod _{i=1}^{N}\displaystyle \frac{1}{{Z}_{i}}{\left({M}_{i}/{M}_{c}\right)}^{\alpha }\exp (-{M}_{i}/{M}_{c}){p}_{\mathrm{obs}}({M}_{i}| \tau ),\end{eqnarray} \tag{ 10 }$$

where the normalization becomes

$$\begin{eqnarray}&&{Z}_{i}={\int }_{{M}_{50}}^{\infty }{\left({M}_{i}/{M}_{c}\right)}^{\alpha }\exp (-{M}_{i}/{M}_{c}){p}_{\mathrm{obs}}({M}_{i}| \tau ,{N}_{\mathrm{MS}}).\end{eqnarray} \tag{ 11 }$$

Power-law Model: For comparison, we also fit for a pure power-law form of the mass function, as opposed to the Schechter form. The corresponding likelihood function for the power-law model is given by:

$$\begin{eqnarray}&&p({M}_{i}| \theta ,\tau )\equiv \prod _{i=1}^{N}\displaystyle \frac{1}{{Z}_{i}}{\left({M}_{i}\right)}^{\alpha }{p}_{\mathrm{obs}}({M}_{i}| \tau ),\end{eqnarray} \tag{ 12 }$$

where the normalization becomes

$$\begin{eqnarray}&&{Z}_{i}={\int }_{{M}_{50}}^{\infty }{({M}_{i})}^{\alpha }{p}_{\mathrm{obs}}({M}_{i}| \tau ,{N}_{\mathrm{MS}}).\end{eqnarray} \tag{ 13 }$$

Schechter function fitting results for 533 young clusters with masses greater than the local 50% mass completeness limit are shown in Figure 6. The left panel shows the completeness-corrected mass distribution for the observed sample and a Schechter function overplotted in blue that uses median α and M_c values derived from marginalized posterior PDFs. We draw 100 random (α, M_c ) samples from the two-dimensional posterior PDF and plot their corresponding Schechter forms in gray to indicate the variance in the fitted parameters. The binned histogram of observed masses is used only for visualization purposes, and not used for Schechter function fitting.

Zoom In Zoom Out Reset image size

We find the CMF to be well described by a Schechter function with $\alpha =-{2.06}_{-0.13}^{+0.14}$, and log(M_c /M_⊙) $=\,{4.24}_{-0.13}^{+0.16}$. These results are based on marginalized one-dimensional posterior PDFs shown in the right panel of Figure 6. We discuss this primary result further and place it into context with previous results in Section 5.

4.1.1. CMF Dependencies on Cluster Properties

We test for age dependence of the Schechter function fits by dividing the cluster sample into age bins of 10–100 Myr and 100–300 Myr, as shown in the top panel of Figure 7. A signature of mass-dependent cluster destruction would be a flattening of the low-mass slope of the CMF with increasing age due to the dissolution of low-mass clusters (Gieles 2009). Using clusters with masses greater than the local 50% mass completeness, we fit 254 clusters with ages between 10–100 Myr and 279 clusters with ages between 100–300 Myr. We present Schechter function parameters for the two age bins in the bottom panel of Figure 7.

Zoom In Zoom Out Reset image size

In addition to the median values presented in Figure 7, we find the maximum a posteriori best-fit values to be α = −2.04 and −2.03, and log(M_c /M_⊙) = 4.39 and 4.07 for the young and old samples, respectively. Both values are well within the 1σ range of the PDFs.

We do not observe any flattening in the slope of the CMF with increased age. There is therefore no clear evidence for mass-dependent cluster destruction being important over timescales of 300 Myr in the central regions of M33. We do notice an increase in M_c for the younger sample. In comparing the PDFs, we find an 87% probability that the younger sample has a higher M_c than the older sample. This could be evidence of a nonconstant SFH, with enhanced star formation in the last 100 Myr relative to the 100–300 Myr time period. However, we note that the constraints on the best-fit Schechter function parameters are significantly broader than for the full sample due to the decrease in number statistics in each age bin. Therefore, drawing strong conclusions about M33's SFH based on these measurements is difficult.

In addition to age, the CMF could also depend on galactic radius. Previous studies have shown a radial dependence on Schechter function truncation mass, where M_c decreases with galactic radius (Adamo et al. 2015), and in M33 there is evidence of a radial dependence on stellar age, with older ages at the center (Williams et al. 2009; Davidge & Puzia 2011). Unfortunately, this work is unable to answer this question. The region of M33 imaged by the PHATTER survey only extends to galactic radii of ∼5 kpc, and we only find 92 clusters with radii >3 kpc. Fits to radially binned samples yielded no significant trends due to these limited number statistics and the small range of galactic radii probed. Determination of any CMF radial dependence in M33 will require cluster samples that extend to larger radii.

4.1.2. Analyzing our Assumptions in Completeness

In addition to fitting a Schechter function, we fit the M33 cluster sample for a power-law distribution. Power-law fitting results are shown in the left panel of Figure 8. We find the distribution of clusters with masses greater than the local 50% mass completeness limit to be best described by α = −${2.49}_{-0.06}^{+0.06}$. This slope is much steeper than the canonical −2 power law in the literature (e.g., Krumholz et al. 2019), but still overpredicts the number of high-mass clusters by >4σ.

Zoom In Zoom Out Reset image size

Previous studies have used the number of massive clusters to test whether cluster masses follow a power law or Schechter function distribution (Johnson et al. 2017; Adamo et al. 2020). Implementing this same test here, we observe 16 M33 clusters above the fitted truncation mass of log(M_c /M_⊙) = 4.24. To compare the models to this observation, we calculate this same number for 10,000 synthetic distributions of clusters drawn from the PDFs of both the Schechter function and power-law fits presented above.

The predictions from these synthetic distributions are shown in the right two panels of Figure 8. In both panels, it is evident that the Schechter function distribution better matches the actual number of high-mass clusters (i.e., 16). The middle panel shows that the median number of synthetic clusters observed above log(M/M_⊙) of 4.24 for the Schechter function is ${18}_{-6}^{+7}$, and ${35}_{-5}^{+7}$ for the power law, where uncertainties give the 16th and 84th percentiles of the distributions. All 10,000 draws from the power-law distribution result in more than the 16 observed high-mass clusters. The Schechter function provides a much better fit to this population of high-mass clusters. Tests using other threshold masses indicate that these results are not sensitive to our exact threshold choice.

The right panel of Figure 8 shows the cumulative distribution of clusters above the mass on the x-axis. We perform a Kolmogorov–Smirnov (K-S) test on these distributions for both Schechter and power-law distributions to assess the overall goodness of fit to observed data. Due to inherent randomness in drawing distributions, we draw 100 distributions and report the median values of cluster counts above a given mass to compare cumulative distributions to observed counts, as shown in the middle panel of Figure 8. The K-S statistic for the Schechter function is 0.04 with a p-value of 0.99, indicating the observed distribution is consistent with the Schechter function. However, for the power-law distribution, the K-S statistic is 0.21 with a p-value of 0.02. Thus, the test suggests that the two samples are not drawn from the same distribution. In addition to the K-S test, we run an Anderson–Darling test, which suggests a Schechter function is consistent with our mass distribution with p > 0.25; on the other hand, the power-law distribution has a p-value of 0.008, and thus is not consistent with the data.

One final test we use to determine which model best describes the cluster sample is the Bayesian Information Criterion (BIC) test. We find a delta BIC of 6.28. We use the Kass & Raftery (1995) criterion, which classifies this delta BIC as strong evidence in favor of Schechter function parameters opposed to the power-law model. We also perform this test on the two age samples in Section 4.1.1. Due to the decreased number statistics, the delta BIC is less significant than for the full sample at 5.54 and 5.64 for the young and old samples, respectively, but each sample nonetheless provides positive evidence in favor of the Schechter function according to the Kass & Raftery (1995) guidelines.

Both the number of high-mass clusters, and mass distribution of clusters are therefore not well described by a power-law distribution, thus providing strong evidence for the existence of a truncation at higher masses.

Johnson et al. (2017) found a clear correlation between the truncation mass of the CMF and a galaxy's SFR surface density, Σ_SFR. In this section we assess whether M33 follows this trend, and compare our results to theoretical predictions of mass function truncation.

Following the methodology described in Appendix A of Johnson et al. (2017), we measure a characteristic Σ_SFR for M33. We combine GALEX far-UV and Spitzer 24 μm images following the prescription of Leroy et al. (2008) to produce a map of Σ_SFR. We use this particular SFR prescription for consistency with other galaxy measurements, but note that in the future, a CMD-based SFH estimate will be available for the PHATTER survey region (M. Lazzarini 2022, in preparation). We use an SFR-weighted average surface density, 〈Σ_SFR〉, to summarize the local, kiloparsec-scale properties of M33's disk in a way that accounts for the nonuniform distribution of star formation. We measure log(〈Σ_SFR〉/(M_⊙ yr⁻¹ kpc⁻²) of $-{2.04}_{-0.18}^{+0.16}$ for M33 within the PHATTER survey footprint.

We compare M_c and 〈Σ_SFR〉 for M33 to a compilation of measurements for other galaxies in Figure 9. We combine results from M31 (Johnson et al. 2017), M51 (Messa et al. 2018), M83 (Adamo et al. 2015), NGC628 (Adamo et al. 2017), the Antennae (Jordan et al. 2007), as well as four galaxies from the HiPEEC survey (NGC 3256, NGC 3690, NGC 4194, and NGC 6054; Adamo et al. 2020). For 〈Σ_SFR〉 measurements, we adopt the 80% Σ_SFR values tabulated in Adamo et al. (2020) for the HiPEEC galaxies, while we use measurements from Johnson et al. (2017) for M31, M51, M83, and the Antennae. While the HiPEEC 80% Σ_SFR values are not a perfect match to the SFR-weighted 〈Σ_SFR〉 measurements discussed above, this alternative form of area-weighting tends to produce a similar, focused measure of Σ_SFR.

Zoom In Zoom Out Reset image size

Johnson et al. (2017) found a nearly linear relation between M_c and 〈Σ_SFR〉 represented by,

$$\begin{eqnarray}&&\mathrm{log}({M}_{c})=(1.07\pm 0.10)\times \mathrm{log}(\langle {{\rm{\Sigma }}}_{\mathrm{SFR}}\rangle )+(6.82\pm 0.20).\end{eqnarray} \tag{ 14 }$$

We plot this relation as the black line in Figure 9. In Johnson et al. (2017), the fit was based on just four data points: M31, M51, ¹¹ M83, and the Antennae. Here, we add six new data points that show remarkable agreement with the original fitted trend from Equation (14), with a mean and maximum residual of 0.24 dex and 0.81 dex, respectively. We also use the Python package linmix (Kelly 2007) to fit this new sample of observations and constrain the relation's intrinsic scatter, accounting for uncertainties in M_c and 〈Σ_SFR〉. We derive a slope of 0.97 ± 0.13 and an intercept of 6.56 ± 0.18, which are consistent with the values fit by Johnson et al. (2017). We constrain a median intrinsic scatter of only ${0.19}_{-0.15}^{+0.39}$ dex around the newly fit relation, which is small relative to the 4 dex range of both axes.

In addition to the observational results discussed above, Reina-Campos & Kruijssen (2017) present a theoretical model that extends the cluster formation efficiency (Γ) prescription from Kruijssen (2012) and predicts a maximum cluster mass based on three environmental observables: the gas surface density, the epicylic frequency, and the Toomre Q parameter. We use rotation curve and radial surface density profiles of total gas and stars compiled in Utomo et al. (2019) to compute a prediction for the PHATTER M33 cluster sample. We adopt input parameter values ¹² that correspond to a characteristic galactic radius of 2.2 kpc, the median galactocentric radius of the PHATTER cluster sample.

The Reina-Campos & Kruijssen (2017) model predicts a feedback-limited maximum cluster mass for M33 of log(M/M_⊙) = 4.06, ∼0.2 dex below our fitted M_c value. Similarly, we calculate a model prediction for M31 using observable input values ¹³ applicable to the dominant 10 kpc "Ring-Total" region from the PHAT Γ study by Johnson et al. (2016), whose properties apply to a vast majority of that sample's clusters. The Reina-Campos & Kruijssen (2017) model predicts a feedback-limited maximum cluster mass for M31 of log(M/M_⊙) = 3.53, 0.4 dex below the fitted M_c value. While the model appears to underestimate the maximum cluster mass in both M31 and M33, the availability of an environmentally dependent model that is accurate within a factor of two to three is a fantastic step toward better understanding cluster formation.

The nearest galaxy to M33 in 〈Σ_SFR〉 is NGC628; however, there is ∼1 dex difference in values of log(M_c /M_⊙). Despite the similar 〈Σ_SFR〉, these two galaxies are quite different, and these differences may account for their different M_c values. Specifically, NGC 628 has a log(M_⋆) of 10.2 (Cook et al. 2014), nearly an order of magnitude higher mass than M33. Similarly, NGC 628's peak rotation curve velocity is ∼180 km s⁻¹ (Aniyan et al. 2018), about 60% higher than in M33 (Utomo et al. 2019). Also, NGC 628's central metallicity is [O/H] = +0.1 (Berg et al. 2015), roughly 0.2 dex higher than M33. Further, as we will discuss in Section 5, the large mass errors on the individual cluster measurements could lead to a significant upward bias in the measurement of M_c in NGC 628 that may account for some of the observed disagreement with M33.

Our CMF fitting method ignores the individual cluster mass errors that result from CMD analysis. For small errors (≲0.05 dex), this simplification has been shown to minimally impact the mass function fitting (Weisz et al. 2013). However, larger errors may result in significant biases in our inferred parameters. In particular, integrated light measurements typically have larger mass errors than the CMD-based mass estimates we have used here. In this section, we explore the potential biases in mass function parameters as a function of the cluster mass error distribution by using literature samples of both CMD and integrated light measurements. We expect large mass errors to lead to the overestimation of M_c values, due to the steeply declining mass function, which results in more clusters being scattered to higher mass values.

5.2.1. M31 Results: Impacts on the CMF from Integrated Light Masses

We first consider data from M31 where cluster masses have been derived with both CMD and integrated light methods. The CMD masses in M31 (Johnson et al. 2016) were derived using an identical technique to that described in Section 2. The cluster ages used range between 10 and 300 Myr. The integrated light ages and masses were fit by M. Fouesneau (private communication) using the method described in Fouesneau et al. (2014). The authors utilize a large number of simulated clusters that are sampled using MCMC to develop a PDF for each cluster's age, mass and extinction. We show the direct comparison of CMD masses to the integrated light masses in Figure 10. Because we use the same sample of clusters that have CMD ages, for our mass function fitting, we use the same completeness function of Johnson et al. (2017) for both cluster samples.

Zoom In Zoom Out Reset image size

We quantify the distributions of reported uncertainties for the two methods in Figure 11. The median 1σ errors on the cluster masses from the CMD method are just 0.03 dex. However, for the integrated light measurements, the median 1σ errors are 0.14 dex. Further, the median ratio of integrated light uncertainty to CMD uncertainty is 3.3. While the difference in uncertainties between the two methods is large, in Figure 10, we do not see a significant bias in the one-to-one comparison of cluster masses, with the median difference of −0.06 dex. However, there is a fair amount of scatter around this median with a standard deviation of 0.59 dex.

Zoom In Zoom Out Reset image size

The Schechter function fits to the CMD-based masses are presented in Johnson et al. (2017). We use an identical method to fit the integrated light ages for the same sample of clusters; the fit results are shown in Figure 12. The best Schechter function has a power-law index of α = −1.91 ± 0.08 and a truncation mass of log (M_c / ⊙ ) = ${4.35}_{-0.12}^{+0.15}$. For the CMD masses, Johnson et al. (2017) found α = −1.99 ± 0.12, with log(M_c /M_⊙) = ${3.93}_{-0.10}^{+0.13}$.

Zoom In Zoom Out Reset image size

Based on the analysis for M31, we find that the inferred value of M_c is 0.4 dex higher when masses are measured using the less-accurate integrated light method rather than the higher-accuracy CMD masses. This result is very striking. A 0.4 dex increase in M_c translates to a factor of ∼2.5 in cluster mass, and biases of this level could significantly impact the trend of M_c with Σ_SFR. It is important to note that the only difference between the two measurements is the method in which the ages and masses were derived.

5.2.2. Examining M_c Biases in Published Mass Function Fits

To test if cluster mass uncertainties are the cause of the observed bias in the previous section, in this section, we run a series of tests fitting mock samples of clusters drawn from a Schechter function, incorporating the effects of mass uncertainty. In our first experiment, we examine the effect of a fixed mass error for each cluster, in the second, we insert single high-mass outliers, while the third experiment incorporates the observed error distributions and Schechter function parameters of galaxies with published CMF fits. Overall, we find that large individual cluster mass uncertainties lead to a significant positive bias for inferred M_c values.

We supplement clusters from M31 with data from M51, M83, and NGC628. The integrated light cluster mass data from these catalogs are from the M83 results in Adamo et al. (2015), and the LEGUS results (Calzetti et al. 2015): Adamo et al. (2017) for NGC628, and Messa et al. (2018) for M51. A detailed description of the cluster identification and mass measurements is given in Adamo et al. (2017). Due to the distance of these galaxies (∼4 Mpc, ∼5× that of M31 and M33), clusters appear only partially resolved, and thus only integrated light measurements are possible. We use the Bayesian mass estimates of the clusters derived using the SLUG code (Krumholz et al. 2015). Their sample includes clusters with masses above 5000 M_⊙.

In all cases, we assume the mass uncertainty estimates to be Gaussian in log(M), where the standard deviation of the Gaussian corresponds to the (84th percentile–16th percentile)/2 of the PDF. We note that while we refer to these errors as mass errors, they are actually errors in log(M/M_⊙). We incorporate these individual cluster mass errors into simulated mass function fits.

For our first experiment, we create a synthetic distribution of 1000 clusters with α, and log(M_c /M_⊙) values of −2 and 5.0, respectively, consistent with the best-fit parameters of typical literature galaxies. We fit the initial distribution as a baseline, then add a uniform Gaussian mass error to each cluster. We then refit a Schechter function to the new distribution of scattered masses. Finally, we compare the fitted M_c parameter to the value of the original distribution. We run 100 trials of this experiment for a fixed mass error of 0.05 dex, and 0.15 dex, which covers the range of error distributions we observe in the literature.

For a mass error of 0.05 dex, we find a median log(M_c /M_⊙) bias of 0.01 dex above the original value, while for the 0.15 dex mass errors, we get a log(M_c /M_⊙) bias of 0.09 dex. When we compare these biases to the median 1σ range for the derived log(M_c /M_⊙) parameters of 0.30 dex, these biases are not outside of normal uncertainty. This base experiment shows that for small errors, the impact on Schechter function parameters is minimal, but that larger errors produce more bias in the Schechter function truncation masses.

During this first experiment, we found a handful of our simulated clusters were scattered to very high masses. To understand the effect that a single high-mass outlier has on the best-fit CMF, we conduct a second experiment. We use the same setup as in our first experiment with no mass errors added. We then add a single, high-mass cluster starting from the maximum mass cluster drawn from the Schechter function (which had log(M/M_⊙) = 5.35), and steadily increasing this mass to log(M/M_⊙) = 6.35. The results are shown in Figure 13. As expected, the single high-mass cluster raises the M_c value inferred, with very significant biases seen for maximum mass outliers more than 1 dex higher than the rest of the distribution. We run this experiment for a range of input M_c values and total number of clusters that reflects values of the galaxies we discuss here (i.e., M31, M33, M51, M83, and NGC628) and find these results to be consistent across input parameters. This shows the significant impact that cluster mass errors can have on the highest-mass clusters. While a single high-mass cluster can impact the inferred M_c value, the majority of the observed bias is not due to single outliers.

Zoom In Zoom Out Reset image size

For each of the three samples presented in Figure 13, we calculate the delta BIC for Schechter parameters relative to the pure power-law model. Not surprisingly, we find that with increased number statistics, the delta BIC favoring a Schechter model is slightly larger. However, regardless of number statistics, the single outlier causes negligible changes to the delta BIC, and in all cases, the delta BIC is larger than 6, providing strong evidence for a Schechter function model over a power-law model according to the Kass & Raftery (1995) guidelines. While the Schechter function is still preferred, for large bias, the inferred M_c values are very poorly constrained with 1σ uncertainties upwards of ∼1 dex. These constraints are worse for the sample with the fewest clusters.

Our third experiment simulates the real distribution of uncertainty in published mass estimates by generating synthetic cluster distributions from Schechter function parameters published in the literature and assigning each synthetic cluster a mass error of an associated real cluster. We then add a Gaussian random mass error to each cluster, and fit a Schechter function to the new distribution of scattered masses running an MCMC fit. Finally, we compare the fitted M_c parameter to the value fitted to the original distribution. We run 100 trials per galaxy. We note that we are able to assign a real cluster mass error to a random synthetic cluster because we do not find a correlation between cluster mass and mass error in any of the cluster samples considered here.

We run this experiment for M83 and NGC628 using the Schechter function fits from Adamo et al. (2015, 2017), respectively, while for M51, we use the results in Messa et al. (2018). We also run the experiment on our set of CMD-derived cluster properties in M31 (Johnson et al. 2017), and M33 (this work), as well as the SED integrated light cluster properties for M31 clusters discussed in Section 5.2.1. For each galaxy, we report the median output log(M_c )–input log(M_c ) for the 100 trials, along with the 1σ confidence interval.

We present the results of this experiment in Figure 14. We find that when we simulate realistic cluster mass uncertainties and fit the CMF, the M_c value recovered is typically higher than the input value. The resulting M_c bias is <1 dex in all cases, with the median bias being comparable to the reported 1σ uncertainties in M_c in all cases. The median bias is highest in NGC 628 and for the integrated light measurements of M31 due to these samples' larger cluster mass errors. NGC 628 also has the smallest number of clusters, while the galaxy with the largest number of clusters (M51) shows a very small bias. We also find that the measured bias in log(M_c ) for M31 clusters discussed in the previous subsection falls within the distribution of simulated differences.

Zoom In Zoom Out Reset image size

We have seen that both large average cluster mass errors and single mass outliers can significantly impact our estimate of M_c . Thus, the distribution of cluster mass errors, not just the average error, is important. To quantify the error distributions of our literature cluster samples, we perform an Anderson–Darling test and find that these distributions can each be described by a log-normal distribution. A log-normal distribution is parameterized by a mean, μ, and standard deviation, σ, as represented by:

$$\begin{eqnarray}&&\displaystyle \frac{1}{x\sigma \sqrt{2\pi }}\exp \left(-\displaystyle \frac{{\left(\mathrm{ln}(x)-\mu \right)}^{2}}{2{\sigma }^{2}}\right).\end{eqnarray} \tag{ 15 }$$

For a log-normal's given μ, σ, we calculate the predicted M_c bias, shown in Figure 15. As in our first experiment, we assume 1000 cluster masses drawn from Schechter function parameters α = −2 and M_c = 5.0 and fit this distribution as a baseline. Here we add cluster mass errors drawn from each log-normal distribution before refitting a Schechter function. We then take the median bias in M_c from 100 trials. These values are the circular points in Figure 15.

Zoom In Zoom Out Reset image size

We find the model to have a clear gradient of increased bias with respect to both log-normal parameters. Thus, either having overall larger errors or having a tail of clusters with large errors can lead to significant biases in the inferred M_c . We also find the model to be consistent with the biases we observe in published fits.

However, we note that our test grid does not exactly match the number statistics of the individual measurements, which could lead to the small differences we see between the inferred values for the individual galaxies and our model grid. We verified that the calculated bias is insensitive to the input M_c .

In practice, the M_c bias could be avoided by incorporating cluster mass uncertainties into the derivation of the CMF with a hierarchical Bayesian model. For large cluster samples with small uncertainties (like for the M33 CMD mass estimates presented here), this more complicated procedure is not necessary. However, it should be considered for smaller samples, or those that rely on integrated light measurements.