A Novel Method to Automatically Detect and Measure the Ages of Star Clusters in Nearby Galaxies: Application to the Large Magellanic Cloud

The archival data used in this work were acquired from several diverse large surveys, which mapped the Magellanic Clouds at various bands (some examples are shown in Figure 1). Starting from shorter wavelengths, Simons et al. (2014) composed a mosaic using archival data from the Galaxy Evolution Explorer (GALEX; Martin et al. 2005) at the near-ultraviolet (NUV) band (${\lambda }_{\mathrm{eff}}=2275$ Å). The mosaic covers an area of 15 deg² on the LMC. All exposures of the same field were co-added to improve the signal-to-noise ratio. The median exposure time was 733 s, and the 5σ depth varied between 20.8 and 22.7 mags. Unfortunately, this mosaic does not cover the central ∼3 × 1 deg² of the LMC (the bar-region), due to detector limitations. This area was later observed by the Swift Ultraviolet-Optical Telescope (UVOT) Magellanic Clouds Survey (SUMAC; Siegel et al. 2014), which covered an area of ∼4 × 2 deg² around the bar-region, with typical exposures of 3000 s in all three NUV filters of this instrument (UVW1, UVW2, and UVM2).

Zoom In Zoom Out Reset image size

The optical data used here are from the Magellanic Cloud Photometric Survey (MCPS; Zaritsky et al. 2004). These authors observed the central 64 deg² of the LMC with 3.8–5.2 minute exposures at the Johnson U, B, V, and Gunn i filters of the Las Campanas Swope Telescope. Typical seeing was 1.5 arcsecond, with limiting magnitudes that varied, depending on the filter, between 21.5 mag for U and 23.0 mag for i.

Meixner et al. (2006) performed a uniform and unbiased imaging survey of the LMC (called Surveying the Agents of a Galaxy's Evolution, or SAGE), covering the central 7 deg² with both the Infrared Array Camera (IRAC; Fazio et al. 2004) and the Multiband Imaging Photometer (MIPS; Rieke et al. 2004) on-board the Spitzer Space Telescope. The SAGE survey produced mosaics at 3.6, 4.5, 5.8, and 8.0 μm for IRAC, and at 24, 70, and 160 μm for MIPS. The exposure times varied between 43–60 s, depending on the band, and added up to a total of 291 hr for IRAC and 217 hr for MIPS. Our current analysis has been performed on the area covered by SAGE, and we cropped accordingly the mosaics of GALEX, SUMAC, and MCPS surveys.

Using DAOPHOT II (Stetson 1987), Zaritsky et al. (2004) elaborated a photometric catalog, which contains 24.5 million sources in the area covered by MCPS. They also estimated the line-of-sight extinctions of the stars in their catalog to produce an extinction map of the LMC. To that end, they compared the observed stellar colors with those derived by the stellar photospheric models of Lejeune et al. (1997). Thus, they measured the effective temperature (${T}_{\mathrm{eff}}$) and the extinction (A_V) along the line of sight to each star, adopting a standard Galactic extinction curve. They produced two A_V maps, one for hot ($\mathrm{12,000}\,{\rm{K}}\lt {T}_{\mathrm{eff}}\leqslant \mathrm{45,000}$ K) and one for cool ($\mathrm{5,500}\,{\rm{K}}\lt {T}_{\mathrm{eff}}\leqslant \mathrm{6,500}$ K) stars.

Using these extinction maps, we correct the observed colors of the cluster candidate and field comparison stars, after separating them into two categories depending on their color (and thus their ${T}_{\mathrm{eff}}$). Stars having (B − V) ≤ 0.20, which corresponds to A5 or earlier type stars with ${T}_{\mathrm{eff}}\geqslant \mathrm{78,00}$ K, are classified as hot and corrected with the hot star map, whereas for stars with (B − V) > 0.20, we use the cool star map. We also adopt R_V = 3.1 and relative extinctions ${A}_{\lambda }$/A_V from Schlegel et al. (1998).

Detecting spatial clustering of stars in the Galaxy as well as in nearby galaxies, where individual stars can be observed, is a challenging task. Since star clusters are not found in isolation and many times are projected over very crowded fields (e.g., on the central bar of the LMC), it is nearly impossible for a visual search to distinguish cluster members from field stars. An automated method that is based on statistical analysis and Monte-Carlo simulations is arguably more prone to succeed. Schmeja (2011) presented a comparison between four different cluster detection algorithms, which we briefly describe in the following: (1) the star counts, (2) the nearest neighbor, (3) the Voronoi tessellation, and (4) the separation of the minimum spanning tree. The star counts method simply counts stars located in a region of interest and detects overdensities above some local background threshold that is defined by the user. The nearest neighbor method estimates the local source density by measuring the distance of each object to its nth nearest neighbor. A Voronoi tessellation partitions an (x,y) plane with n points into n polygons, such that each polygon contains only one point, and then defines the local source density as the reciprocal of the area of the polygon around this point. Finally, the separation of the minimum spanning tree traces a unique set of lines (called edges), connecting a given set of vertices without closed loops, such that the sum of the edge lengths is minimum; by applying the desired length threshold, star clusters can be defined. As Schmeja (2011) explains, he carried out Monte-Carlo simulations of synthetic star clusters using simple power-law, Gaussian, or King spatial stellar distributions. Then he examined the performance of the aforementioned algorithms and concluded that, while distinct centrally concentrated clusters are detected by all methods, those with low overdensity or highly hierarchical structure are only reliably detected by methods with inherent smoothing, i.e., the star counts and nearest neighbor algorithms. Based on these considerations, and taking into account that the latter requires significantly larger computational time (×200 longer), the star counts method was selected as the optimal approach for our analysis.

Images of reasonable resolution and depth (i.e., allowing to observe and resolve individual stars) are crucial for the identification of the candidate clusters. After detecting the positions of all the stars in an image through a source extraction code (e.g., SExtractor; Bertin & Arnouts 1996), we create a pixel-map, where each star (or extended object) is represented by a single pixel. This is done to remove the effects of seeing and extended sources (e.g., galaxies) from the subsequent analysis. All these pixel stars are assigned the same arbitrary flux value, since their spatial density, and not the photometric properties, is the important parameter for the cluster detection sequence. Once the conversion of the observed image to a pixel-map is done, the star cluster detection code is applied. The main virtue of the code is that it only requires to adjust two parameters: the size of the region of interest and the detection threshold. The region of interest is a box where the density of stars is counted and compared to the local background density, which is estimated over a much larger box, whose size is not critical as long as it samples the local background density well. After accurate testing, we observed that the ideal size of the region-of-interest is approximately that of a cluster at the distance of the LMC (∼60 arcsec in diameter; Nayak et al. 2016), since it minimizes the detection of false associations.

A more critical parameter to set is the detection threshold (${{\rm{\Sigma }}}_{\det }$). As described in Schmeja (2011), large values of ${{\rm{\Sigma }}}_{\det }$ will be able to detect real overdensities only in low background density regions. On the other hand, low ${{\rm{\Sigma }}}_{\det }$ values may be able to extract clusters in high background density regions (e.g., the LMC bar), at the cost of detecting many false associations (results of random projections) in the lower density regions, such as the outskirts of the galaxy. For this reason, our code is using a variable ${{\rm{\Sigma }}}_{\det }$ value that changes as a function of the local background density. To define the relation between ${{\rm{\Sigma }}}_{\det }$ and the background density, we performed Monte-Carlo simulations with artificial clusters having both Gaussian as well as uniform overdensity profiles (accounting for both compact and diffuse clusters), projected over various background values. In Table 1, we list the background densities (both high and low), and the rates of artificial cluster recovery and false association detection, in all the different filters we used during the extraction sequence. Our variable ${{\rm{\Sigma }}}_{\det }$ is able to maintain a constant detection rate.

Using this ${{\rm{\Sigma }}}_{\det }$ calibration, we run the detection code to identify star clusters. At this step, pixel stars, which are considered to be cluster members, are recorded, while pixel stars deemed as background "counts" are excluded from further analysis. Once we are left with the "good" pixel-map (i.e., where each pixel-star is a bona-fide member of a cluster), we need to determine the cluster center and radius. To do so, we use again SExtractor, which is able to detect flux emitted by coherent pixel groups, and provide their centroid and size. We therefore need to "fool" SExtractor into interpreting our good pixel-map as a flux map. We achieve this by smoothing the good pixel-map with a large kernel (similar to that of the point-spread function), to redistribute the arbitrary flux in each pixel-star over the neighboring pixels. This will give a flux-like image over which we can run SExtractor using its intuitive setup. Since the flux-like image is generated from the pixel-map that contains only cluster members, it does not suffer from any ambiguity related to the background. In practice, we recreated the original image, but background subtracted (where "background" in this context means "field stars"). We have to stress here that, due to stellar resolution incompleteness in the most crowded regions of massive clusters (representing <1% of our sample), there might be an overestimation of the final cluster radius of the order of 10%–20%.

To maximize the detection of as many clusters as possible (both young and old), we apply our detection sequence to images at different wavelengths (see Figure 2). For this purpose, we have used imaging from GALEX and Swift in the UV (that probes the hot massive stars), and the Spitzer IRAC, which is dominated by old stellar populations and low-mass stars. Our catalog contains 5459 cluster candidates in an area covering the central 7 deg² of the LMC. For each cluster center and radius (defined by the SExtractor), we find all the stars within by cross-correlating with the extinction corrected MCPS catalog (see Section 2.2). This procedure has produced a photometric catalog of 1.9 million stars located in our candidate clusters.

Zoom In Zoom Out Reset image size

A crucial step in defining a cluster CMD is the proper accounting of the field-star contamination. Here we use a process similar to that described in Mighell et al. (1996). For each cluster candidate CMD, our code also produces a field-star CMD, using all the stars contained in a box of 0.16 deg² around the cluster center but excluding any stars located ≤0.05 deg from the cluster. This is done to sample well the surrounding field-star CMD without including any unaccounted cluster stars located outside the cluster radius. The CMDs are then binned along both axes using the color and magnitude uncertainties—in our case, we used δ(color) = 0.5 mag and δ(magnitude) = 1 mag. The code estimates the membership probability (${p}_{\mathrm{memb}}$) of each star to belong to the given cluster, by considering the number of stars populating the corresponding bins in the cluster and field CMDs. In particular, we adopted the following formula:

$$\begin{eqnarray}&&{p}_{\mathrm{memb}}=1-\displaystyle \frac{{N}_{* ,\mathrm{field}}}{{N}_{* ,\mathrm{cluster}+\mathrm{field}}}\cdot \displaystyle \frac{{A}_{\mathrm{cluster}}}{{A}_{\mathrm{field}}},\end{eqnarray} \tag{ 1 }$$

where ${N}_{* ,\mathrm{field}}$ is the number of field stars (in the 0.16 deg² box), and ${N}_{* ,\mathrm{cluster}+\mathrm{field}}$ is the total number of stars contained within the cluster radius (i.e., a combination of field and cluster stars) in each color–magnitude bin. ${A}_{\mathrm{field}}$ and ${A}_{\mathrm{cluster}}$ are the sizes of the areas from which the field and cluster+field stars are extracted, respectively. Each cluster star candidate is hence assigned a ${p}_{\mathrm{memb}}$, which will later be used by the age determination code. An example of the results of this procedure for the candidate cluster IR1-297 is presented in Figure 3; stars with ${p}_{\mathrm{memb}}\geqslant 0.90$ are represented with blue circles, in green circles are those with 0.60 < p_memb ≤ 0.90, while all other stars (having p_memb ≤ 0.60) are shown as red circles. Most stars in the central denser region of the cluster are assigned a higher membership probability, while green and red circles are mainly encountered in the periphery of the cluster. Therefore, our CMD decontamination is also consistent with the expected radial distribution of the clusters. It is also noteworthy that during this decontamination procedure we excluded from the field region all those stars that might belong to some neighboring clusters, hence keeping clean and unbiased the field-star CMD (e.g., see the empty space in the bottom-right corner of Figure 3). Finally, the code discards as a false detection any cluster candidate with fewer than 20 stars with ${p}_{\mathrm{memb}}\geqslant 0.60$. We imposed this additional step to eliminate any false associations detected during the cluster extraction process; the number of candidates found from star counts thus decreased from 5459 to 4850 (≈11% reduction).

Zoom In Zoom Out Reset image size

To determine the ages of our clusters, we used a Bayesian approach to obtain the most likely theoretical isochrone that reproduces the observed CMD. The method is analytically described in Hernandez & Valls-Gabaud (2008), and Walmswell et al. (2013). In this method, each theoretical isochrone is treated as a probability density function (PDF) on the CMD. Therefore, the probability that a cluster star i located at (${x}_{i}\pm {\sigma }_{x,i}$, ${y}_{i}\pm {\sigma }_{y,i}$) will come from a model isochrone n is

$$\begin{eqnarray}&&{p}_{i}=\int {\rho }_{n}(x,y){U}_{i}(x-{x}_{i},y-{y}_{i}){dxdy},\end{eqnarray} \tag{ 2 }$$

where ${U}_{i}(x-{x}_{i},y-{y}_{i})$ is the error function for the star i and ${\rho }_{n}(x,y)$ is the PDF along the model isochrone n. The error function of the star is taken as a bivariate Gaussian:

$$\begin{eqnarray}&&{U}_{i}(x-{x}_{i},y-{y}_{i})=\displaystyle \frac{1}{2\pi {\sigma }_{x,i}{\sigma }_{y,i}}{e}^{-\left(\tfrac{{(x-{x}_{i})}^{2}}{2{\sigma }_{x,i}^{2}}+\displaystyle \frac{{(y-{y}_{i})}^{2}}{2{\sigma }_{y,i}^{2}}\right)},\end{eqnarray} \tag{ 3 }$$

where the star i is associated with errors ${\sigma }_{x,i}$ and ${\sigma }_{y,i}$. Finally, the total likelihood of the cluster to be drawn from the model isochrone n is given by the product of probabilities over the S stars belonging to the cluster:

$$\begin{eqnarray}&&{L}_{n}=\displaystyle \prod _{i=1}^{S}{p}_{i}\end{eqnarray} \tag{ 4 }$$

The complete set of such likelihoods (L_n) provides the PDF of the age of the given cluster. The assumed cluster age is the maximum of the PDF, with lower and upper uncertainties calculated as the 16th and 84th percentiles, respectively (e.g., see Figure 4). The code we use was originally developed to infer the star formation histories of systems of few observed stars (such as ultra-faint dwarf galaxies and star clusters), and is analytically described in V. H. Ramírez-Siordia et al. (2017, in preparation), where multiple comparisons with real and simulated clusters are performed to test its accuracy, as well as the effects on estimated ages from the uncertainties on metallicity, distance modulus, and extinction. A difference with the original code is that our version uses ${p}_{\mathrm{memb}}$, described in Section 4.1, for the estimation of Equation (2) (${p}_{i}^{{\prime} }={p}_{i}\cdot {p}_{\mathrm{memb}}$). Therefore, the final likelihood of the star to belong to the isochrone will be given by the product of the two probabilities. Moreover, the code is allowed to vary the distance modulus by ±0.25 mag (this corresponds to ≈11 kpc), in order to account for distance variations of the clusters in the LMC (Subramanian & Subramaniam 2009). To ensure a robust age estimation, we fit the (U − V) versus V, (B − V) versus V, and (V − i) versus i CMDs of each cluster, and we combine the results as described in the next section. A comparison between the ages derived from the three different CMDs, yields median ratios of 1.0 to 1.1 with a standard deviation of the ratios ∼5.5.

Zoom In Zoom Out Reset image size

The isochrones we used are a byproduct of an independent project by S. Charlot & G. Bruzual (2017, in preparation).⁸ These authors have assembled complete sets of PARSEC evolutionary tracks computed by Chen et al. (2015) for 16 values of the stellar metallicity, ranging from Z = 0.0001 to Z = 0.06, complemented with the work by Marigo et al. (2013) to follow the evolution of stars through the thermally pulsing asymptotic giant branch (TP-AGB) phase. The isochrone synthesis algorithm described by Charlot & Bruzual (1991) is used to build isochrones from the evolutionary tracks at any age. Whereas the galaxy spectral evolution models by Charlot & Bruzual use a large number of empirical and theoretical stellar libraries to describe the spectrophotometric properties of the stars along these isochrones (e.g., Gutkin et al. 2016; Wofford et al. 2016; Vidal-García et al. 2017), for the purpose of this investigation, we use the BaSeL 3.1 atlas (Westera et al. 2002) to obtain the stellar UBVi magnitudes. The effects of dust shells surrounding TP-AGB stars on their spectral energy distribution is treated as in González-Lópezlira et al. (2010). The BaSeL 3.1 atlas covers uniformly the $({T}_{\mathrm{eff}},\mathrm{log}\ g)$ plane at the LMC metallicity of Z = 0.008 ($[{Fe}/H]=-0.34$), according to studies of Cepheid LMC stars (e.g., Romaniello et al. 2005; Keller & Wood 2006). Clusters with ages ∼2–3 Gyr are expected to have lower metallicities (e.g., Z = 0.006; Piatti et al. 2015), but this difference is insignificant compared to the uncertainties in the estimation of their age. Our final grid of 80 isochrones covers the range of 6.9 ≤ log(age) < 9.7 years.

Our final catalog contains 4850 clusters and is presented in Table 2. Column (1) gives the cluster ID assigned by our detection code. Columns (2) and (3), respectively, report the right ascension (R.A.) and declination (decl.) of the cluster centers, in J2000 decimal equatorial coordinates. The cluster radii are presented in column (4). Finally, columns (5), (6), and (7) contain the best age estimation for each cluster, as well as its lower and upper uncertainty bounds (from the 16th and 84th percentiles of the PDF).

The best age for each cluster is calculated as a combination of the ages resulting from the fitting of the aforementioned three CMDs (i.e., (U − V), (B − V), and (V − i)), after weighting them by both the number of stars that were used in each fit and the total uncertainty of the corresponding age measurement. This is synopsized by the following formula:

$$\begin{eqnarray}&&\mathrm{Best}\ \mathrm{Age}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{3}\tfrac{{N}_{i}\cdot {\mathrm{Age}}_{i}}{(\delta {\mathrm{Age}}_{i})}}{{\sum }_{i=1}^{3}\tfrac{{N}_{i}}{(\delta {\mathrm{Age}}_{i})}},\end{eqnarray} \tag{ 5 }$$

where i = 1–3 correspond to the results of the (U − V), (B − V), and (V − i) CMD fits, respectively. N_i is the number of stars fitted in the i-th CMD, Age_i is the best age estimation in that CMD, and δAge_i is the difference between the upper and lower uncertainties of Age_i. In Figure 5, we present some characteristic examples of clusters from our catalog, ordered by increasing age. Although most of them were detected in the same band in which they are displayed (i.e., the Spitzer IRAC 3.6 μm), clusters (a) and (b) were only found in the UV bands, thus stressing the importance of using multiple bands in the identification procedure.

Zoom In Zoom Out Reset image size

The cross-correlation of our new cluster catalog with those of Bica et al. (2008) and Werchan & Zaritsky (2011) shows that we identify 3451 new clusters in the 7 deg² we surveyed.

In Figure 6, we compare the distribution of the radii of our clusters (estimated the Gaussian full width at half maximum of each cluster's counts), shown with a solid blue line, with that of the the semimajor axes of the clusters from the sample of Bica et al. (2008), presented with a dashed red line. Although they do not represent the exact same parameter, both distributions peaking around 5–10 pc, with our maximum being more pronounced. Our method fails to detect clusters with radii smaller than 4 pc, a range in which Bica et al. (2008) found 3% of their clusters. At the distance of the LMC, 4 pc correspond to only 8 arcseconds. Since the spatial resolution of the images we used for the cluster identification varies from 2 to 5 arcseconds (0.5–1.25 pc), clusters with $r\leqslant 4$ pc will be hard to identify unambiguously. To ensure that our code is not subject to some selection bias, we estimate the radii of our clusters in an independent way. In particular, we use the results of the membership probability (presented in Section 4.1) to find, for each cluster, the distance of the furthest cluster star with ${p}_{\mathrm{memb}}\gt 0.60$, and then we adopt this as an alternative cluster radius. We over-plot these values in Figure 6, with a dotted-dashed green line and find that it is almost identical to our initial radii estimation. Another important difference with Bica et al. (2008) is the absence of the minor peak at ∼21–22 pc. Such associations could originate after the expulsion of gas during cluster formation, which usually results in rapid mass loss and eventually their dissolution (e.g., Pfalzner 2009). Hence, the detection of these objects is very improbable due to their short lifetimes and low luminosities. Alternatively, this might also be a selection bias introduced by having set an upper limit (i.e., the region-of-interest) on the cluster radii during the detection process.

Zoom In Zoom Out Reset image size

In Figure 7, we compare our age estimations for clusters in common with other surveys. Glatt et al. (2010) estimated ages by visually identifying the clusters' MSTOs. They also attempted to decontaminate their CMDs, again by visual means. In fact, the visual identification may underestimate the ages of some clusters, since the existence of blue stragglers or field stars projected on the extrapolation of the main sequence to brighter magnitudes might confuse the observer. Only a statistical method that fits the integrated CMD can overcome such limitations. There is nonetheless, as one can see in panel (a) of Figure 7, a good correlation between the two works (with Pearson R = 0.75). In panel (b) of the same figure, we compare our age estimations with those from Popescu et al. (2012), who derived the ages for their clusters by comparing their observed colors with those from Monte-Carlo simulations. Although their method is more sophisticated than visual estimations, unfortunately, it fails to remove the contamination of the field stars and therefore does not appear to correlate with our results better than the previously considered work by Glatt et al.; it yields R = 0.71. Piatti et al. (2015) used high-resolution VISTA Magellanic Cloud (VMC) survey near-infrared CMDs to estimate the ages of ∼300 clusters in the LMC bar and the 30 Dor region. They visually fitted theoretical isochrones on the (Y − K) versus K CMDs, using a sophisticated field-star decontamination technique (Piatti et al. 2014) that takes into account density variations in the field. Comparing their age estimations with ours yields R = 0.76. Baumgardt et al. (2013) selected the best age estimates from previous publications to compile their sample. Their results correlate slightly better with our measurements, having R = 0.79. Finally, Asa'd et al. (2016) estimated the ages of 27 massive LMC clusters by fitting their integrated cluster spectra with theoretical models (Bruzual & Charlot 2003). Although their sample is very small, their results are in excellent agreement with ours, with R = 0.99.

Zoom In Zoom Out Reset image size

With cluster ages and positions available, one can identify prominent peaks in their formation history. This can provide important information about the LMC–SMC–Galaxy interactions in space and time. In Figure 8, we present the age distribution of our clusters. Since histogram peaks can result as artifacts of the binning scheme, we adopted our bin size using the unbiased Freedman–Diaconis rule; its value is 0.075 dex (well above the resolution of our isochrones, which is 0.03 dex). The histogram shows the existence of at least one significant peak $\approx 310\,\mathrm{Myr}$ ago, with secondary ones 160 and 500 Myr ago. Moreover, a smaller increase appears at ≈10–20 Myr. Such enhancements in the cluster formation can be associated with an LMC–SMC direct collision (about 100–300 Myr ago) predicted by the models of Besla et al. (2012), or tidal interactions between the Clouds and the Galaxy. According to Besla et al. (2007), the LMC is still in its first passage about the Galaxy, and the related tides might have triggered star formation during the past gigayear. Similarly, Harris & Zaritsky (2009) studied the star formation history of the LMC, by comparing the MCPS CMDs with theoretical models; they suggested various peaks of star formation (at 12, 100, 500 Myr, and 2 Gyr), many of which are also consistent with our analysis.

Zoom In Zoom Out Reset image size

Interestingly, the LMC bar (selected loosely by us to lie between R.A. 76^h–85^h and decl. −705 to −690; see left panel of Figure 9) harbors fewer peaks than anywhere else in the galaxy (right panel of the same figure): there, cluster formation appears only at ∼10 and 310 Myr ago.

Zoom In Zoom Out Reset image size

Moreover, only 50% of the clusters we detect are older than $\approx 300\,\mathrm{Myr}$. Although it is true that the MSTOs of such clusters are well below our detection limit, their scarcity is probably real, owing to cluster dissolution. Star clusters are destroyed through various mechanisms, such as residual gas expulsion (e.g., Baumgardt & Kroupa 2007), two-body relaxation, and external tidal fields (e.g., Baumgardt & Makino 2003), as well as via tidal heating from shocks, and the harassment from the passage of giant molecular clouds (e.g., Gnedin & Ostriker 1997; Gieles et al. 2006). According to Baumgardt et al. (2013), as many as 90% of the clusters with ages greater than 200 Myr can be dissolved per dex of lifetime: this could explain the decline of the cluster population older than ∼500 Myr in our sample.

Albeit cluster dissolution might explain why the age distribution resembles a Gaussian, one can see that additional populations might be drawn by fitting its exact shape using multiple components. For that purpose, we used a code for Bayesian analysis of univariate Gaussian mixtures (NMIX⁹ ). This code implements the approach presented in Richardson & Green (1997). The results suggest that the above distribution can be fitted either by three, four, or five components (see Figure 10); these solutions have Bayes K-factors between them of ∼1, while those including more or fewer components fail to provide good fits (when comparing the successful models to each one of the rejected univariate distributions, we get K-factors > 5). In Figure 10, we show the three successful mixture models (in dashed green lines,) as well as their individual components (solid black lines). It is worth noticing here that the use of variable width in normal (Gaussian) components is a reasonable approximation, since it is expected that the age uncertainties will dominate over the length of the individual star cluster formation events. According to the results of NMIX, the LMC has experienced various epochs of cluster formation over the past gigayear, with the most probable 10, 310, and 400 Myr ago; these events seem to agree with the visually identified peaks presented above. Additionally, possible secondary peaks exist 70 Myr, 250 Myr, and 1.1 Gyr ago.

Zoom In Zoom Out Reset image size

Important information can also be drawn from the spatial distribution of the cluster ages. In Figure 11, we present that distribution (in decimal equatorial coordinates), with ages separated into eight different bins of varying sizes. As shown, younger clusters (≤50 Myr) are mostly distributed in the northern and northeastern regions of the LMC (where the arm is located), as well as around its bar; they tend to gather in the regions indicated by Kim et al. (2003) to host HI supershells (see Figure 6 of that article), and in the star-forming region 30 Dor. On the other hand, older clusters (250–500 Myr) appear very concentrated in the LMC bar, and relatively uniformly distributed outside of that. These results are consistent with an inside-out formation scenario for the LMC—at least for the past 1 Gyr. Cioni (2009) concluded likewise, based on the observed flattening of the metallicity gradient in the LMC. She argued that, as the galaxy builds-up inside-out, star formation moves toward the outer parts of the galaxy with the result that those regions become enriched by metals and the metallicity gradient flattens out (see also Vlajić et al. 2009, and references therein).

Zoom In Zoom Out Reset image size

The information presented above can be used not only to study the integrated star formation history of the LMC, but also to constrain dynamical simulations of the interactions among the LMC, the SMC, and the Galaxy. It is also important to infer the properties of the clusters, their relation with the environment, and how they might compare to clusters in other galaxies. In a forthcoming paper, we will apply our novel technique to study the cluster populations of the SMC and other nearby galaxies.