Ages and Masses of Star Clusters in M33: a Multi-wavelength Study

We combine archival images for the nearby galaxy M33 (Triangulum Galaxy) from the ultraviolet (UV) to the infrared (IR) to derive ages, masses, and extinctions for the young star cluster population, and compare our physical parameters with published ones. Our goal is to test the robustness of clusters' ages and masses, and possibly improve on existing ones both by expanding the wavelength range of the spectral energy distribution (SED) fits and by using more recent population synthesis models. The rationale for this experiment is to verify the sensitivity of the clusters' physical parameters to observational setups and model choices that span those commonly found in the literature. We derive the physical parameters of 137 clusters, using SEDs measured in eight UV-to-I bands, including H-alpha, from GALEX and ground-based images. We also add the 24 micron image from the Spitzer Space Telescope to help break some age degeneracies. We find that our derived cluster ages show significant differences with earlier determinations, while the masses remain relatively insensitive to the fitting approach adopted. We also highlight an already known difficulty in recovering old, low-extinction clusters, as SED fitting codes tend to prefer younger, higher extinction solutions when the extinction is a free parameter. We publish updated ages, masses, and extinctions, with uncertainties, for all sample star clusters, together with their photometry; given the proximity of M33, this represents an important population to secure for the study of star formation and cluster evolution in spirals.


INTRODUCTION
Within galaxies, the structures of star formation are a continuous, scale-free hierarchy from parsecs to kiloparsecs (Lada & Lada 2003;Elmegreen 2003;Bressert et al. 2010), which are expected to arise from the self-similar distribution of a turbulence-dominated ISM (Elmegreen & Efremov 1997), mediated by magnetic fields and outflow feedback (e.g., Krumholz et al. 2019b). Star clusters form in the dense regions of the hierarchy (e.g., Elmegreen 2010), and, thus, provide a sensitive probe of the star formation process, while the rest of the hierarchy is quickly dispersed into the stellar field of the galaxy, over timescales of a few tens of Myr (Bastian 2008;Longmore et al. 2014;Grasha et al. 2015Grasha et al. , 2017a. Young star clusters contain the majority of the massive stars, which drive the feedback into the surrounding medium and regulate star formation (Krumholz et al. 2019b;Adamo et al. 2020), but most are not bound because they are not dense enough (e.g. Brown & Gnedin 2021). These clusters are dispersed and randomized over short timescales ( 10 Myr) by a vast array of processes, both internal-such as gas expulsion, stellar evolution, and two-body relaxation,-and external,-such as tidal shear, random motions, interactions with molecular clouds, and secular evolution of the host galaxy (e.g., Gieles & Bastian 2008).
Deriving accurate ages and masses for star clusters is a key step for linking the small scales of the star formation, i.e., the scales of individual stars, with the large scales of whole galaxies. The physical and chemical properties of star clusters constrain scenarios for their formation and evolution, trace the star formation history of galaxies, constrain models of cloud collapse, provide tests for the origin and persistence of spiral arms, and guide cosmological simulations (Shabani et al. 2018;Li et al. 2018;Adamo et al. 2020;Ballesteros-Paredes et al. 2020). The distribution of cluster masses is a sensitive tracer of formation mechanisms (e.g. Whitmore et al. 1999;Larsen 2002;Gieles et al. 2006b,a;Whitmore et al. 2014;Johnson et al. 2017;Adamo et al. 2020), while the age distribution reveal the cluster disruption mechanisms (e.g., Gieles 2009;Bastian et al. 2012;Silva-Villa et al. 2014;Chandar et al. 2016;Adamo et al. 2017;Messa et al. 2018).
M33 (Triangulum Galaxy) is an excellent laboratory for testing techniques for the derivation of physical parameters of star clusters. This nearby spiral galaxy (850 kpc; Ferrarese et al. 2000) is a member of the Local Group and has a large population of clusters. M33 is close enough that its star clusters can be resolved into stars with the Hubble Space Telescope, but distant enough that lower angular resolution facilities can survey it fully. Its modest inclination (∼52 o , from NED 1 ) minimizes the chance of line-of-sight confusion. At the M33 distance, an angular aperture of 1 radius subtends a 4.1×5.2 pc 2 spatial region, which is comparable to the size of star clusters Brown & Gnedin 2021). Sarajedini & Mancone (2007) compiled a list of compact sources found in M33, for a total of 451 objects, 255 of which were confirmed to be star clusters. For these clusters, ages and masses were determined in a series of papers by Ma et al. (Ma et al. 2001, 2002a,b,c, 2004a, using spectral energy distribution (SED) fitting of medium-band optical filter photometry.
Starting from this compilation, we combine broad-band optical imaging with the two UV GALEX images, add the hydrogen recombination line Hα, and include the 24 µm band map from the Spitzer Space Telescope, all from archival holdings. Our goal is to leverage the extended wavelength coverage of our photometric measurements to produce accurate ages and masses for the cluster population in this iconic galaxy. The UV is a sensitive discriminator of young versus old populations and, combined with the optical bands, adds a leverage arm for extinction measurements , while Hα tracks the presence of ionizing massive stars. Finally, the 24 µm band traces the emission from the dust heated by the UV photons of young stars ( 100 Myr, Kennicutt & Evans 2012). Thus the 24 µm emission enables discriminating intermediate age populations (∼10-100 Myr) when multiple solutions are possible but the Hα emission is no longer present and the UV is mostly yielding degenerate results . The metal abundance of M33 is sufficiently high that we expect significant dust and 24 µm emission in correspondence of its young star clusters. U et al. (2009), Bresolin (2011) andToribio San Cipriano et al. (2016) report values between 1/2 solar and ∼20% above solar 2 for both the young stars and the gas within the inner ∼3.5 kpc of this galaxy (0.4 R 25 ), where most of our sources are located ( Figure 1). We compare our results with those of Ma et al. to investigate similarities and differences.
The outline of the paper is as follows. Section 2 describes the cluster compilation of Sarajedini & Mancone (2007). Section 3 discusses the archival imaging data, our processing and photometry. Section 4 describes the models employed for SED fitting and our fitting approach. Section 5 presents the results of the SED fits and our analysis, including comparisons with previous determinations of ages and masses. Section 6 summarizes our findings and conclusions. We present our updated ages and masses and extinctions, as well as photometry, in tabular form for all star clusters in our sample.

CLUSTER SAMPLE
The clusters in our sample come from a compilation of Sarajedini & Mancone (2007) of 451 compact sources in M33, among which the authors identified 255 star clusters using high angular resolution imaging from either the Hubble Space Telescope or ground-based facilities. Of these 255 confirmed star clusters, SM2007 lists both ages and masses for a total of 163, as published in a series of papers by Ma et al. (Ma et al. 2001, 2002a,b,c, 2004a. Ma et al. derived ages and masses by fitting the population synthesis models of Bruzual & Charlot (1993, with the 1996 to the SEDs of the sources. The SEDs consisted of photometry in 8 intermediate-band  Massey et al. 2006) showing 154 of the 163 star clusters from Sarajedini & Mancone (2007) analyzed in this work. The remaining nine clusters are located in the Northern and Southern pointings, and are not shown here. A scale of 6 correspond to ∼1.5 kpc at the distance of this galaxy. Cyan circles mark clusters for which SED fits were performed, while red circles mark clusters excluded from the fits, and for which only photometry is provided. See text for more details.
filters in the wavelength range 3800-10,000Å (Ma et al. 2001). This wavelength range does not enable the authors to constrain uniquely the internal dust reddening of the star clusters, which they either adopt from previous analyses-as in Ma et al. (2001), where the authors use the extinction values of Chandar et al. (1999)-or fix E(B-V) to a constant value (E(B-V)=0.1, Ma et al. 2002a). The values of the ages and masses are given without uncertainties. In their SED fits, Ma et al. allowed the metallicity to remain a free parameter, ranging between 3% solar and 1.4 times solar, with the reference solar metallicity being Z =0.02 (Iglesias et al. 1992).
Although optical photometry for the sources is published in SM2007, we do not use these authors' measurements as we aim to derive uniform-aperture photometry for the star clusters, for the purpose of using internally-consistent measures for our SED fitting experiment.

The Data
Ultraviolet, optical and IR images from a combination of space and ground-based imaging data of the M33 galaxy were retrieved from NED, for a total of 10 bands: Far-UV (FUV), Near-UV (NUV), U, B, V, R, I, Hα (6563Å), 3.6µm and 24µm. The FUV and NUV images trace the youngest stars, with ages 100 Myr, while the ionized gas emission at Hα traces even younger stars, 7-10 Myr, although this age limit increases by about a factor of 3 in the presence of binary stars (Xiao et al. 2018). The IR image at 24 µm traces the emission of the dust heated by the young stars (Draine & Li 2007;Calzetti et al. 2007).
The GALEX satellite (Martin et al. 2005) observed M33 in both the FUV (∼1530Å) and NUV (∼2310Å) channels with a single pointing of its ∼1.6 o ×1.6 o field-of-view (FoV) with angular resolution ∼4 -5 (e.g., Dale et al. 2009). The optical band images were obtained by Massey et al. (2006, for U, B, V, R, and I) and Massey et al. (2007, for Hα) at the Mayall 4-m telescope at KPNO, with the Mosaic Camera, with a 38 ×38 FoV and angular resolution between 0 .8 and 1 . To observe the entirety of M33, images were taken in three different pointings: Northern, Central, and Southern. Mosaics at the wavelengths 3.6 µm and 24 µm were obtained with the Spitzer Space Telescope IRAC and MIPS cameras, respectively, with angular resolution of 1 .7 at 3.6 µm and 6 .5 at 24 µm (Fazio et al. 2004;Rieke et al. 2004).
All images were aligned and resampled to the pixel scale of the optical ones (0 .27 per pixel), to preserve the highest possible angular resolution. This facilitates the identification of the location of the SM2007 sources onto the images, and the assessment of whether neighboring sources are present which may affect the ages and mass determinations. The latter step is made necessary by the significantly lower resolution of the GALEX and Spitzer /MIPS images relative to the ground-based ones.

Photomery
Multi-wavelength photometry was performed on the images using the positions published in the SM2007 catalog. An aperture of 3 in radius was adopted as a compromise between enclosing as much as possible of the UV and IR point spread functions (PSFs) while at the same time minimizing the overlap between photometric apertures of neighboring star clusters. An outer annulus of 1 .1 width around each aperture was used to estimate the background and remove it from the aperture measurement. A consequence of the small aperture size is to exclude the wings of the UV and IR PSFs, which leads to an underestimate of the flux of the sources in these bands. We correct for this effect by applying aperture corrections derived by studying the growth curves of isolated sources in the images, and listed in Table 1. The optical broad-band photometry does not require aperture corrections.
We do not apply aperture corrections to the narrow-band photometry, although it is evident from visual inspection that the ionized gas emission is sometimes more extended than the size of our photometric aperture. This aperture subtends a physical radius of 12 pc at the distance of M33, which corresponds to the Strömgren radius of a 10 4 M , 5 Myr old star cluster, for an assumed electron density n e ∼20 cm −3 . In general, our clusters are either less massive or older, but some overfill the aperture with their ionized gas emission. In order to offset this limitation, after running the SED fits, we visually inspect the continuum-subtracted Hα image in correspondence of each cluster to ensure that the age results are not driven by an underestimated Hα flux. The units of the UV and optical images are count-rates and counts, respectively. We convert those units to physical units of flux density (erg s −1 cm −2Å−1 ) using: the stellar photometry of Massey et al. (2006) with the zeropoints published in Zombeck (1990, p. 100) in the case of the UBVRI images, the stellar photometry of Massey et al. (2007) with the zeropoints published in the KPNO Mosaic instrument manual 3 in the case of the narrowband Hα image, and the conversion factors published in the GALEX instrument manual 4 in the case of the UV images. The list of flux conversion factors is in Table 2. Measured photometry for all confirmed 163 star clusters is listed in Table 6. Fluxes are given prior to correction for foreground Milky Way extinction, which amounts to E(B-V)=0.036 (Schlafly & Finkbeiner 2011).
To directly observe the presence of line emission in our images, we subtract the stellar continuum as traced by the R-band image from the narrow-band filter targeting the Hα line. We rescale the R image prior to subtraction, with the scaling factor derived by computing a ratio of the counts in R to the counts in the narrow-band for emissionline-free stellar sources. The narrow-band filter has a FWHM=80.62Å (Massey et al. 2007) and includes the two [NII](6548,6584Å) metal lines. We use the ratio [NII]/Hα=0.27 published in Kennicutt et al. (2008) to remove the [NII] contamination from the narrow-band filter and obtain an Hα image. The stellar continuum subtracted Hα images are used solely for visual comparison to confirm the ages of the younger star clusters (<10 Myr) through the presence or absence of excess Hα emission. The original, un-subtracted Hα images are used for photometry and SED fitting.
The archival Spitzer images are already calibrated in surface brightness units of MJy/sr. We convert those to fluxes (erg s −1 cm −2 ) within our apertures using standard conversion factors, with the appropriate pixel scale of 0 .27 per pixel. The 3.6 µm photometry is further multiplied by the factor 0.91 to include photometric corrections 5 . The IRAC 3.6 µm image is only used to remove the underlying stellar emission from the 24 µm image, via the formula: F 24,dust [Jy]=F 24 [Jy] -0.035 F 3.6 [Jy] (Helou et al. 2004;Calzetti et al. 2007). The photometry in these two bands is reported in Table 6 in units of flux. We use the 24 µm image mainly to break age degeneracies in the SED fits when the cluster is sufficiently old that the ionized gas emission is no longer present, but still young enough to have UV emission heating the dust.
Visual inspection of all the measured star cluster reveals that 26 out of 163 cannot be retained within our sample for SED fitting for a variety of reasons, as detailed in Table 7. The most common reason is the presence of multiple sources with different colors (likely different ages) within the 3 radius aperture used for the photometry. We thus limit the SED fits to the remaining 137 sources for which we are reasonably confident that one single source is contained within the photometric aperture, or, in case of multiple sources, these show similar colors (likely similar/equal ages).

Measurement Uncertainties
The images at U, B, V, Hα, R, and I are sufficiently deep that flux calibration uncertainties dominate the error budget. These are taken from Massey et al. (2006) for the broad-band filters and from Massey et al. (2007) for the Hα filter. The uncertainties are listed in Table 3 as uncertainties on the logarithm of the flux. For the optical filters, the value quoted is the largest uncertainty of the sample in each filter.
Uncertainties in the FUV, NUV, and in the (stellar continuum-subtracted) 24 µm photometry are dominated by uncertainties in the aperture placement, because of the large PSF size. We calculate these uncertainties by 'wiggling' the position of the aperture by one pixel (0 .27) in each direction and remeasuring the flux within the aperture. The scatter that results is used in the final error budget, listed in Table 3.

Models
The stellar population synthesis models used by Ma et al. (2001) and subsequent papers, i.e., the models of Bruzual & Charlot (1993) Table 3. For the optical filters, we quote the calibration error from Massey et al. (2006) and Massey et al. (2007); the value shown is the largest uncertainty for clusters in the luminosity range of our sample; most star clusters have smaller uncertainty in the optical. UV and 24µm flux errors are calculated through variations of the location of the photometric aperture's center, which dominate the uncertainty budget.
versions exist, such as those of Bruzual & Charlot (2003) and more recent ones, which also include several updates to the stellar populations characteristics (e.g. Gutkin et al. 2016). Wofford et al. (2016) compared the most recent versions of these models with those generated by Starburst99 Vázquez & Leitherer 2005) with the addition of nebular lines from YGGDRASIL (Zackrisson et al. 2011), and the BPASS ones (Stanway et al. 2020), finding that all models yield similar result, with low dispersion in the median values of age, mass, and extinction. In this analysis, we adopt the Starburst99+YGGDRASIL simple stellar population (SSP) models, with the assumption that star clusters can be reasonably approximated by an instantaneous burst (single) population (Wofford et al. 2016). The models are deterministic, meaning that the stellar Initial Mass Function (IMF) is assumed to be fully sampled from the lowest to the highest stellar masses. This is not the case for star clusters with masses<3,000-5,000 M (e.g., Cerviño et al. 2002), implying that the age and mass determinations for clusters below this mass limit will carry an uncertainty larger than the formal one we quote. While the correct approach for such low mass clusters is to use a Bayesian fitting method that returns full posterior probability distributions for a stochastically sampled IMF (e.g., SLUG Krumholz et al. 2019a), this approach is beyond the scope of our current work.
Starburst99 ) spectral synthesis models were generated with a Kroupa (2001) IMF in the range 0.1-120 M and metallicity Z = 0.02. This metallicity value is higher than the mean value of the M33 region where most of our clusters are located, ∼0.5-1 solar, but we choose it for the following three reasons: (1) the model 'ingredients' are inconsistently sampled for metallicity below solar (Vázquez & Leitherer 2005); (2) the next model metallicity available is about 1/2 solar, at the lower end of the range for our clusters; and (3) ages and masses derived via multi-band SED fits are relatively insensitive to variations of a ∼2.5 factor in metallicity of the models, as shown in Figures 13 and 14 (see, also, de Grijs et al. 2005, and SM2007). We do expect star clusters older than a few 100 Myr to have lower metallicity than younger ones, as also found by Ma et al. (2001). However, our study is mostly interested in star clusters younger than ≈1-3×10 8 yrs, which is the age range where we expect the addition of UV, Hα, and dust emission (24 µm) tracers to have the most impact on the clusters' physical parameters determination.
The models utilize the Padova tracks with AGB treatment, to better represent intermediate and old stellar populations (Girardi et al. 2000;Vázquez & Leitherer 2005) and are generated for the age range 1 Myr-13 Gyr; the age steps are 1 Myr in the 1-15 Myr range, 10 Myr in the 20-100 Myr range, 100 Myr in the 200-1000 Myr range, and 1 Gyr in the 2-13 Gyr range. The addition of nebular emission lines via YGGDRASIL (Zackrisson et al. 2011) includes a covering factor of 50% for the gas emission, to account for leakage of ionizing photons out of HII regions (e.g., Calzetti et al. 2021). The models assume non-rotating, single stars. There is, however, increasing evidence that stellar population models implementing rotating, binary stars are better fits for data of star clusters, especially at young ages and low metallicities (Stanway et al. 2020). Variations at the young ages can be as large as a factor 2-3 (Wofford et al. 2016). However, given that M33 star clusters span a range of ages much larger than this scatter and that the galaxy is fairly metal rich, we will not include effects of rotating and binary stars, although this should be noted as a point for future investigations. We note that the Starburst99+YGGDRASIL models do not include Very Massive Stars (>150 M ) either, which are known to be present in young star clusters (e.g. Crowther et al. 2010;Smith et al. 2016) and are included in more recent populations synthesis models (e.g. BPASS, Eldridge et al. 2017;. Our goal in choosing the Starburst99+YGGDRASIL models is to remain close to the astrophysical assumptions used in the stellar population synthesis models used by Ma et al. (2001) and subsequent papers.
The effects of dust on the stellar population SEDs are applied via a starburst attenuation curve (Calzetti et al. 2000), with the equation: where k(λ) is the starburst curve, and E(B-V) is the color excess applied to the SEDs. We generate models with the color excess range E(B-V) = 0-1, with step 0.02. As we use the directly measured cluster photometry of Table 6, the small foreground Milky Way extinction with E(B-V) = 0.036 is corrected with the Milky Way curve of Fitzpatrick (1999). The dust-attenuated SEDs are then convolved with the transmission curve of each of the UV and optical filters to produce synthetic luminosities in 8 bands, that are normalized to the default mass of Starburst99, 10 6 M . The synthetic photometry so obtained is then used to fit the observed SEDs of the star clusters. Ma et al. (2001) and subsequent works by these authors estimated the ages and masses of the star clusters in M33 by using a least square fit between observed and synthetic photometry (Kong et al. 2000). We adopt the same approach, described by:

Fitting Approach
is the synthetic flux density of the models described in the previous section in the i-th band, and σ(λ i , n) is the uncertainty in the observed flux density; the model populations are a function of the age t, the color excess E(B-V), and the cluster mass M. The index i runs from 1 to 8, which is the number of bandpasses in our SEDs. The goal is to minimize the reduced χ 2 red =χ 2 /(N f reed − 1), where N f reed =5 is the number of degrees of freedom for 8 independent datapoints and 3 parameters. Uncertainties in the fits are obtained by investigating the distribution of χ 2 red values around the minimum, as described in Calzetti et al. (2021); we also study the χ 2 red distributions to ensure that our best fits are not the result of local minima. Ma et al. (2001) explored solutions for three different metallicities -Z = 0.0004, 0.004, and 0.02, finding that older star clusters are better fit by models with lower metallicity. However, as already discussed above, these authors could not constrain the amount of intrinsic extinction in the clusters, due to the limited wavelength range of their SEDs, and either adopted a constant value of E(B-V)=0.1 or values derived by previous authors. Conversely, we employ a single metallicity value, but extend the wavelength range of the clusters' SEDs to better separate age from extinction. Our goal is to understand how changing the wavelength range of the SEDs and the fitting approach affects the values of ages <1-3×10 8 yr, where the age-extinction degeneracy is generally pronounced. The inclusion of UV colors provides, theoretically, a better handle of the effects of extinction, and, when accompanied by the U−B color, provides a sensitive discriminant for separating age from extinction. With the addition of the Hα emission we can better identify clusters younger than ∼7-8 Myr, a regime where the above colors are often of inconclusive diagnostics (see discussions in Boquien et al. 2009;Calzetti 2013;Calzetti et al. 2015). Furthermore, We use the 24µm emission as a discriminant for multiple-age solutions, i.e., for the presence of multiple minima in the distribution of the χ 2 red values. For clusters with two viable ages with significant difference (e.g. 8 Myr and 80 Myr), a visual check of the 24µm emission at the location of the cluster is used as the deciding factor. The 24 µm dust emission has been shown to follow closely the 8 µm emission in HII regions (Relaño & Kennicutt 2009), and the latter has been shown to decrease dramatically with increasing age of the star cluster out to at least 300-400 Myr (Lin et al. 2020). Sharma et al. (2011) showed that most of the 24 µm bright sources in M33 are generally younger than 10 Myr.

Ages, Masses, and Extinctions
The ages, masses, and extinctions of the 137 star clusters in M33 obtained from the photometry in Table 6 via SED fitting are listed in Table 5, together with the ages and masses compiled by SM2007. Examples of SEDs and their best fits are shown in Figures 2, 3, and 4, for five cases each with ages ≤10 Myr, 10 Myr-100 Myr, and >100 Myr, respectively. As expected, the youngest star clusters usually display Hα in emission and strong UV fluxes relative to the optical ones, the oldest (>100 Myr) clusters have red SEDs with no Hα emission, while the clusters with ages intermediate between these two often display blue SEDs, with significant UV emission, but no or little measurable Hα emission. A few individual cases are shown in Figures 5, 6, and 7 for the young, intermediate, and old stellar population cases; here we display, in addition to the observed SED and its best fit, image cutouts centered on the cluster, in the light of Hα (both continuum-subtracted and non) and 24 µm, to show the cluster's appearance at these wavelengths. Young clusters show clear Hα and 24 µm emission, which decreases in intensity with increasing age. In Figure 8, the case of a star cluster with extended Hα is shown; in this case, the ionized gas emission has developed a 'ruptured bubble' morphology and our photometric aperture only captures a fraction of the Hα emission from this cluster, which explains the faintness of the Hα line in the SED plot. The distributions of ages, masses, and extinctions are shown in graphical form in the two panels of Figure 9, where the masses and E(B-V) are plotted as a function of the ages in the left and right panels, respectively. Statistics for the star clusters are also given in Table 4.  Most (90%) of the clusters younger than 10 Myr have masses below 3,000 M and only one is more massive than 5,000 M (Table 4), indicating that cluster masses in this age range are uncertain due to stochastic sampling of the stellar IMF (Adamo et al. 2017). The general impact of stochastic sampling on our age and mass determinations for these young clusters is that the clusters could be moderately more massive and a factor of a few older than what we find; the discrepancies between deterministic and stochastic models decrease for older ages and higher mass clusters (Krumholz et al. 2015, , their Figure 14). Older age clusters are usually more massive, as seen in the left hand-side panel of Figure 9, where the model line marks the standard lower luminosity limit for the detections. As expected, the masses we determine at a given age are above the model line (SM2007; see, also, Adamo et al. 2017). The trend for older clusters to be more massive than younger ones at the low mass end is simply due to fading with age, while the observation that they are more massive also at the high end of the mass distribution is the well known size-of-sample effect (Hunter et al. 2003). Figure 10 compares our derived masses with those reported by SM2007. The masses scatter around the 1-to-1 line (in log scale) up to ∼10 4 M , and our masses become systematically smaller than those reported in SM2007 at higher values. We attribute this discrepancy to the systematically younger ages we derive from our SED fits, which yields lower masses at a given luminosity. Overall, however, the masses derived in this work and those reported by SM2007 track each other reasonably well. The most massive clusters in this sample are about 2×10 5 M , which is typical of spiral galaxies (Larsen 2009;Adamo & Bastian 2018).
A comparison of the ages derived here with those reported in SM2007 reveals, conversely, major discrepancies, as shown in Figure 11. There are several instances in which clusters that were classified as 'old' (older than ∼100 Myr) in SM2007 are found to be around or younger than 10 Myr with our SED fits and visual checks of the Hα images. At the opposite end, many clusters that are given ages 10 Myr in SM2007 are found to have much older ages by our SED fits, >30-40 Myr. One important characteristic of our results is that we find less than a handful of clusters with ages >10 9 yr, while SM2007 list many. Whitmore et al. (2020) discusses in detail that there appear to be a limitation in the ability of SED fitting routines to identify old clusters correctly: in the presence of a 'red' SED, fitting routines often prefer a highly-extincted, younger age solution to an extinction-free, older age solution. A work-around to this issue is to impose a limit to the highest value of E(B-V) that an older star cluster can have, in order to force an older age solution (as discussed in Whitmore et al. 2020). This approach was adopted by Ma et al. (2004b) for a set of the The four small panels to the right show for the cluster, beginning at the top left and moving clockwise: the Hα emission, the continuum-subtracted Hα emission, an RGB color composite (red = R band, green = Hα band, and blue = U band), and the stellar-subtracted 24µm emission. The magenta circle shows the size of the aperture used for the photometric measurements. The scale bar below the aperture represents 3". The flux scale of the images is linear, however the upper and lower limits differ between bands. M33 clusters, where they impose a constant value E(B-V)=0.1 to derive their ages. We note this limitation, which implies that several of our clusters in the 300 Myr-1 Gyr age range could be much older than what we find.  Artificially 'too-young' ages produce color excess values that are unusually high, or distributions that are uncharacteristically skewed towards high values. The right panel of Figure 9 shows the distribution of the best-fit E(B-V) as a function of cluster age. The maximum E(B-V) value we recover is 0.6 mag, and >2/3 of the clusters have E(B-V)<0.2 mag, as expected for modestly reddened environments. As a reminder, the quoted values of E(B-V) include the contribution of the foreground Milky Way dust, E(B-V)=0.036 mag. The distribution is notably 'flat',  with several clusters yielding ages ≥300 Myr and E(B-V)≥0.4 mag. One such cluster is shown at the top of Figure 4. In this case, the UV-optical SED are not able to provide constraints on the age-dust degeneracy, and the UV emission from the stellar population is weak, implying a weak 24 µm dust emission. Figure 12 shows the distribution of the 24 µm/FUV flux ratio as a function of the FUV-NUV color. The 24 µm/FUV flux ratio measures the fraction of UV light absorbed by dust and re-emitted in the infrared; the FUV-NUV color, here given as the ratio of the flux density at the two wavelengths, measures the amount of dust reddening suffered by the UV SED (Meurer et al. 1999;Calzetti et al. 2000). The star clusters in M33 are consistent with modest-to-negligible reddening and processing of UV light Figure 11. Comparison of our SED-fit ages with those of the SM2007 catalog.
into the infrared, in general agreement with our findings for the E(B-V) values. The red (more negative) FUV-NUV colors at very low values of the 24 µm/FUV ratio are easily explained with effects of age in the clusters SEDs (Calzetti et al. 2005;Cortese et al. 2006). From this analysis, we conclude that, while many of the clusters with ages ≥300 Myr in our sample could be older if they are less dusty than what we find from our SED fits (see discussion in Whitmore et al. 2020), the data at hand do not allow us to unequivocally conclude this. However, our age determinations for the clusters younger than ≈300 Myr are reasonably robust.

Colors versus Age and Extinction
The B-V color becomes redder for increasing age, as discussed by many authors, including SM2007. Our measurements follow this general trend, although we note a significant scatter towards red colors at fixed age ( Figure 13). A large scatter in the direction of redder colors at fixed age is also observed for the NUV-I color, although there is still a general trend of redder colors for increasing age (Figure 14). This color, and similar colors (e.g., NUV-K) that use a large wavelength baseline to magnify effects, was introduced as an age indicator for galaxies once the effects of internal extinction are removed (e.g., Muñoz-Mateos et al. 2007;Cortese et al. 2008). We verify here that the same colors can be used as effective age indicators also in the case of star clusters. As already discussed by other authors, one downside of the long wavelength baseline is the sensitivity of the NUV-I color to dust extinction, as shown both by the scatter noted above ( Figure 14) and by the trend of redder colors for increasing E(B-V) values we observe for our clusters (Figure 15). In this figure we show the locus of a 200 Myr star cluster for increasing values of E(B-V) as a red line. Most clusters in our sample scatter around this line, with <10 clusters deviating from this trend by showing much redder NUV-I colors.
In Figures 13 and 14 we show model expectations for both metallicities Z = 0.02 (red lines) and Z = 0.008 (blue lines). Notably, the differences between the two model lines are small, and well within our measurement and parameter fit uncertainties; this further justifies the use of Z = 0.02 models for our SED fits. As already mentioned above, while the colors follow the general trend of the age lines, the data scatter mostly above these lines; this is particularly evident for the NUV-I color (Figure 14), and is the expected effect of dust extinction.

SUMMARY AND CONCLUSIONS
Leveraging the extended wavelength baseline provided by the GALEX UV together with optical and Hα bands, we have used SED-fitting to derive improved ages, masses, and extinction values for 137 confirmed star clusters in M33.  (Meurer et al. 1999). This correlation has been previously observed for galaxies and HII regions, and is extended here to clusters. The model line, in red, is from Calzetti et al. (2005). Most of the clusters in M33 are consistent with low-to-negligible dust attenuation. The large scatter in the FUV-NUV colors at constant 24µm/FUV value are due to the age differences in the stellar populations of the clusters, with redder colors for older populations. Figure 13. The B-V flux ratio v. age using our measurements, without any extinction correction. Correcting the flux ratio for the effect of the Milky Way foreground extinction would make it bluer (move downward) by 0.036 dex, which is significantly smaller than the scatter in the data.
Identified from the compilation of Sarajedini & Mancone (2007), all clusters also have masses and ages derived in a series of papers by Ma et al. (2001Ma et al. ( , 2002aMa et al. ( ,b,c, 2004a. The 24 µm emission detected by the Spitzer Space Telescope provides an additional discriminant for breaking the age-extinction degeneracy below 100 Myr. We exclude 26 clusters from the original list of 163 for a number of reasons that affect their SEDs, such as presence of multiple sources with different colors within the measurement aperture (Table 7). We compare our results with the earlier ones.
In terms of cluster ages, we find overall agreement with the earlier results, but also noticeable differences. About half dozen clusters that were given ages older than 100 Myr are younger than ∼10 Myr based on our method. We also find ages 100 Myr for about a dozen clusters that were previously given ages 10 Myr. We note a dearth of clusters with ages >1 Gyr in our sample, possibly due to a well-known problem with SED fits which tend to prefer younger ages at higher extinctions over older, less extincted solutions (Whitmore et al. 2020). As a result, we observe a clumping of clusters around ∼200-300 Myr of age.  Unlike ages, our SED-derived masses track the earlier determinations reasonably well. The main discrepancy is that we do not find cluster masses 2×10 5 M , while SM2007 reports a little over a dozen clusters above this mass. We attribute much of the discrepancy to the dearth of clusters older than ∼1 Gyr in our fits. We should remark that the mass range we find for the star clusters in M33 is consistent with the mass range found for star clusters in other spiral galaxies at comparable age intervals (Adamo & Bastian 2018).
All our star clusters have moderate-to-small internal extinction, all with E(B-V)<0.6 mag, and about 2/3 with E(B-V)<0.2 mag. We do not observe a trend between age and extinction, although the fast timescale for the diffusion of star clusters in galaxies (Grasha et al. 2017a,b) would imply that they reach regions of low dust content within several tens of Myr. Aforementioned age-extinction degeneracies for intermediate age clusters (≥200-300 Myr) could explain why we observe an extinction value as high as ∼0.6 mag in a ∼300 Myr old cluster; an older age with a lower extinction may be an alternate, viable explanation.
In addition to providing the list of physical parameters for the 137 clusters whose SEDs we fit, we also publish our 10-band photometry for all 163 clusters, which includes: GALEX FUV and NUV, optical U, B, V, R, Hα, and I, and infrared centered at 3.6 µm and 24 µm. Our photometry table can provide a resource to test alternative methods for deriving clusters' physical parameters that may overcome some of the limitations discussed in this work.  Table 5 continued on next page  Table 5 continued on next page  Table 5 continued on next page  Table 5 continued on next page

±0.03
Older best fit age retained.
Note-SED-fitting derived ages, masses, and colors excesses, E(B-V), together with 1 σ uncertainties for 137 confirmed star clusters in M33, listed together with the ages and masses compiled in Sarajedini & Mancone (2007. The E(B-V) values in this Table are the sum of the internal dust reddening in the star cluster plus the foreground reddening of the Milky Way, E(B-V)=0.036 (Schlafly & Finkbeiner 2011).
In some instances of two or more sources present in the aperture, the best fit results are retained if the sources appear of comparable age (similar behavior in colors), as indicated in the Comments field; in this case, the quoted mass is the sum of all masses for the sources in the aperture. For sources yielding with two best-fits to the age/extinction, the reported value reflects an assessment based on the 24 µm emission. The retained best fit is given in the Comments field.  Table 6 continued on next page  Table 6 continued on next page  Table 6 continued on next page  Table 6 continued on next page