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

T. Bitsakis, P. Bonfini, R. A. González-Lópezlira, V. H. Ramírez-Siordia, G.

Bruzual, S. Charlot, G. Maravelias, D. Zaritsky

► To cite this version:

T. Bitsakis, P. Bonfini, R. A. González-Lópezlira, V. H. Ramírez-Siordia, G. Bruzual, et al. A Novel Method to Automatically Detect and Measure the Ages of Star Clusters in Nearby Galaxies: Application to the Large Magellanic Cloud. The Astrophysical Journal, 2017, 845, 10.3847/1538-4357/aa8090. insu-03747449

HAL Id: insu-03747449

https://insu.hal.science/insu-03747449v1

Submitted on 3 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

T. Bitsakis¹, P. Bonfini¹, R. A. González-Lópezlira¹,²,³, V. H. Ramírez-Siordia¹, G. Bruzual¹, S. Charlot⁴, G. Maravelias⁵,⁶, and D. Zaritsky⁷

¹Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia, 58089, Mexico; t.bitsakis@ crya.unam.mx

²Helmholtz-Institut für Strahlen-und Kernphysik (HISKP), Universität Bonn, Nussallee 14-16, D-53115 Bonn, Germany

³Argelander Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany

⁴Sorbonne Universités, UPMC-CNRS, UMR7095, Institut d'Astrophysique de Paris, F-75014 Paris, France

⁵Instituto de Física y Astronomía, Universidad de Valparaíso, Av. Gran Bretaña 1111, Casilla 5030, Valparaíso, Chile

⁶Astronomický ústav, Akademie věd České republiky, Fričova 298, 251 65 Ondřejov, Czech Republic
⁷Steward Observatory, University of Arizona, Tucson, AZ 85719, USA

Received 2017 May 2; revised 2017 July 5; accepted 2017 July 5; published 2017 August 11

## Abstract

We present our new, fully automated method to detect and measure the ages of star clusters in nearby galaxies, where individual stars can be resolved. The method relies purely on statistical analysis of observations and Monte-Carlo simulations to define stellar overdensities in the data. It decontaminates the cluster color-magnitude diagrams and, using a revised version of the Bayesian isochrone fitting code of Ramírez-Siordia et al., estimates the ages of the clusters. Comparisons of our estimates with those from other surveys show the superiority of our method to extract and measure the ages of star clusters, even in the most crowded fields. An application of our method is shown for the high-resolution, multiband imaging of the Large Magellanic Cloud. We detect 4850 clusters in the 7 deg² we surveyed, 3451 of which have not been reported before. Our findings suggest multiple epochs of star cluster formation, with the most probable occurring ~310 Myr ago. Several of these events are consistent with the epochs of the interactions among the Large and Small Magellanic Clouds, and the Galaxy, as predicted by N-body numerical simulations. Finally, the spatially resolved star cluster formation history may suggest an inside-out cluster formation scenario throughout the LMC, for the past 1 Gyr.

**Key words:** catalogs – galaxies: star clusters: general – Magellanic Clouds – methods: statistical

**Supporting material:** machine-readable table

## 1. Introduction

Star clusters provide a unique tool to understand the star formation history of Local Group galaxies, given their well constrained distances and our ability to determine their ages with good precision. They can also provide significant information about the initial mass function (IMF), the distribution of star formation within a galaxy, and whether star formation occurs everywhere synchronously (Maragkoudakis et al. 2017).

The Large Magellanic Cloud (LMC) is not only one of the nearest galaxies, it is also an exceptional laboratory of galaxy evolution under the influence of strong gravitational interactions. The LMC and the Small Magellanic Cloud (SMC), constitute an interacting pair, which is also bound to the Galaxy. It is believed that the Magellanic Stream—a 180 kpc long intergalactic filament that extends from the LMC to the SMC and through the Galactic south pole—is a relic of strong tidal encounters between the aforementioned galaxies (see Yoshizawa & Noguchi 2003; Besla et al. 2012).

Various authors have attempted to study the formation of star clusters in the LMC. Bica et al. (2008) compiled a general catalog, using all of the previously known sources, as well as their own findings based on visual inspection on Sky Survey plates. Their catalog contains 9305 extended objects in the Magellanic System (the LMC, the SMC and the stream); among these, about 3700 have been identified as clusters/associations in the LMC. Using the Bica et al. catalog, Glatt et al. (2010) estimated the cluster ages through visual fitting of their color-magnitude diagrams (CMDs) with isochrones generated by the Geneva and Padova codes

(Lejeune & Schaerer 2001; Girardi et al. 1995, respectively). They concluded that there were various periods of intense star cluster formation over the last 1 Gyr, the two most prominent of which occurred 125 and 800 Myr ago. In addition, Baumgardt et al. (2013) compiled a new catalog with the ages of 307 clusters by selecting the best age estimations from previous publications. Their data suggest a burst of cluster formation about 1 Gyr ago, with little evidence for the existence of similar activity afterwards. More recently, Nayak et al. (2016) presented a semi-automated method to estimate age and reddening of 1072 star clusters in the LMC. They found at least one major cluster formation event 125 Myr ago. There seems, therefore, to be a consensus that the LMC experienced various episodes of star formation, plausibly compatible—given the uncertainties—with the results of N-body simulations on the tidal interactions between the Magellanic Clouds and/or the Galaxy. Recently, such simulations which made use of the latest proper motion observations, suggested that the Magellanic Clouds are probably in their first passage about the Galaxy (e.g., Besla et al. 2007, 2012; Kallivayalil et al. 2013).

Despite all of this progress, the aforementioned studies exploring the age distribution and properties of star clusters in the LMC did not only demand significant time and effort to be conducted, but they also introduced important biases (as described in Sections 3 and 4). Piatti & Bica (2012) have shown that, when not correctly addressing the field-star contamination, one can detect more false associations than real star clusters—owing to stochastic fluctuations. In this work, we introduce the first, fully automated method able to
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of the multiband data we used in our work: (a) a Spitzer/IRAC 3.6 µm mosaic of the central region of the LMC (Meixner et al. 2006), (b) the GALEX/NUV mosaic without the central bar (Simons et al. 2014), and (c) the central region observed by Swift/UVOT (Siegel et al. 2014).</figcaption></figure>

detect and determine the ages of star clusters in nearby galaxies. The detection is performed using a novel spatial clustering procedure that is described analytically in Section 3, while in Section 4 we present the age identification algorithm. Since this method is fully computer-based, it avoids important biases that may be introduced by the limitations of solely visual identification; it is also more suited to the use of clustering thresholds and multi-wavelength data. Using the results of our algorithms, in Section 5, we present the LMC star cluster properties, as well as their age and spatial distributions. Finally, we briefly give our summary and conclusions in Section 6.

Throughout this work, we assume a distance modulus for the LMC of 18.50 mag (Walker 2012).

## 2. The Data

### 2.1. Imaging

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_{eff} = 2275$ Å). The mosaic covers an area of $15 \text{ deg}^2$ 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\sigma$ depth varied between 20.8 and 22.7 mags. Unfortunately, this mosaic does not cover the central $\sim 3 \times 1 \text{ deg}^2$ 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 $\sim 4 \times 2 \text{ deg}^2$ around the bar-region, with typical exposures of 3000 s in all three NUV filters of this instrument (UVW1, UVW2, and UVM2).

The optical data used here are from the Magellanic Cloud Photometric Survey (MCPS; Zaritsky et al. 2004). These authors observed the central $64 \text{ deg}^2$ 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 \text{ deg}^2$

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.

### 2.2. Photometric Catalog and Extinction Corrections

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_{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 ($12,000 K < T_{eff} \le 45,000 K$) and one for cool ($5,500 K < T_{eff} \le 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_{eff}$). Stars having ($B-V) \le 0.20$, which corresponds to A5 or earlier type stars with $T_{eff} \ge 78,000$ 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).

## 3. The Cluster Detection Methodology

### 3.1. *Selecting the Appropriate Extraction Algorithm*

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.

## 3.2. Implementing the Detection Sequence

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., SEtractor; 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 ($\Sigma_{det}$). As described in Schmeja (2011), large values of $\Sigma_{det}$ will be able to detect real overdensities only in low background density regions. On the other hand, low $\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

<table><caption>**Table 1**
Recovery Rates from the Results of Monte-Carlo Simulations</caption><thead><tr><th>Telescope/Filter</th><th>Background Density<br>Pixel Stars deg<sup>-2</sup></th><th>Artificial Cluster Recovery Rate (%)<sup>a</sup></th><th>False Association Detections (%)<sup>b</sup></th></tr></thead><tbody><tr><td>Spitzer IRAC1</td><td>10<sup>5</sup> (high)<sup>c</sup></td><td>97</td><td>92</td></tr><tr><td>Spitzer IRAC1</td><td>10<sup>4</sup> (low)</td><td>97</td><td>10</td></tr><tr><td>GALEX NUV</td><td>2·10<sup>4</sup> (high)</td><td>100</td><td>3</td></tr><tr><td>GALEX NUV</td><td>10<sup>3</sup> (low)</td><td>95</td><td>13</td></tr><tr><td>Swift UVW1</td><td>5·10<sup>4</sup> (high)</td><td>70</td><td>7</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Fraction of artificial clusters recovered by the code.

<sup>b</sup> Fraction of false associations with respect to the total number of artificial clusters.

<sup>c</sup> The LMC bar is considered a high background density region; low background density regions are more than 3.5 deg away from the bar.

using a variable $\Sigma_{det}$ value that changes as a function of the local background density. To define the relation between $\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 $\Sigma_{det}$ is able to maintain a constant detection rate.

Using this $\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 SEtractor, which is able to detect flux emitted by coherent pixel groups, and provide their centroid and size. We therefore need to “fool” SEtractor 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 SEtractor 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$^2$ of the LMC. For each cluster center
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Details of the candidate cluster detection process near the 30 Doradus region of the LMC. Green circles correspond to detections made on the *GALEX*/NUV data, while red circles correspond to those made on *Spitzer*/IRAC data. The background image is the *Spitzer*/IRAC 3.6 µm from Meixner et al. (2006).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Example of the decontamination procedure on the candidate cluster IR1-297. In this figure, stars are color-coded based on the membership probability ($p_{\text{memb}}$) associated by our code. We see that stars with $p_{\text{memb}} \ge 0.90$ (blue circles) are concentrated in the center of the cluster, surrounded by stars with $0.60 < p_{\text{memb}} \le 0.90$ (green circles). On the other hand, while several field stars (with $p_{\text{memb}} \le 0.60$; shown in red circles) are projected on the cluster, they are mostly distributed away from its center. Notice the empty region in the bottom-right corner; it corresponds to stars masked-out by the code, since they belong to a neighboring cluster.</figcaption></figure>

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.

## 4. Determination of the Star Clusters' Ages

The determination of the age of a star cluster is a rather challenging process and can be performed using various techniques. One can fit the observed color–magnitude diagram (CMD) of a cluster with sets of theoretical isochrones—this method is arguably the most commonly used—or compare its observed integrated colors (or spectrum) with those derived by theoretical modeling. Although in theory star clusters should present well-defined CMDs, in practice there are a number of factors that introduce significant noise and biases. Field-star

<figure><img src="image_6.png" /><figcaption>**Figure 4.** An example of the *V − i* vs. *i* CMD fitting: cluster candidate IR1-297. Blue: isochrone with the highest likelihood that corresponds to 117 Myr age; green: isochrone in the 16th percentile (71 Myr); magenta: isochrone in the 84th percentile (130 Myr). Only stars with $p_{\text{memb}} > 0.60$ are considered.</figcaption></figure>

contamination, both from the LMC and the Galaxy; photometric limitations and incompleteness; metallicity gradients; multiple main-sequence turn-offs (MSTOs); and blue straggler and binary stars are some major complications.

Previous attempts to determine the ages of the LMC star clusters relied mostly on visual identification techniques. For instance, Pietrzynski & Udalski (2000), Glatt et al. (2010), and Piatti et al. (2015) determined the ages of their clusters through visual identification of the MSTO, after eliminating the field stars by examining the CMD of a region around each cluster. Popescu et al. (2012) carried out a comparison of the integrated broadband photometry of each cluster with models, though without performing a field-star decontamination. More recently, Nayak et al. (2016) used a semi-automated method to estimate the ages. Although the final age determination was decided again by visually fitting the MSTO, a quantitative automated method was used to constrain the range of plausible ages and decontaminate the clusters. Finally, Asa'd et al. (2016) used a code to fit the integrated spectra of star clusters with sets of theoretical models, albeit without removing the contribution of field stars from the spectra. In the next paragraphs, we present a novel age determination method: it is a completely automated, Bayesian CMD fitting algorithm.

### 4.1. Decontamination of the CMDs

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 $\delta$(color) = 0.5 mag and $\delta$(magnitude) = 1 mag. The code estimates the membership probability ($p_{\text{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.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Examples of clusters from our catalog presented on the *Spitzer* IRAC 3.6 µm image. The dashed black lines mark the radii, as defined by the star-counts code. (a): cluster NUV-729, age $12.7_{-3.6}^{+9.2}$ Myr; (b) NUV-1711, age $18.2_{-0.4}^{+2.7}$ Myr; (c) IR1-1718, age $53.7_{-6.9}^{+10.8}$ Myr; (d) IR1-1700, age $81.2_{-7.1}^{+9.9}$ Myr; (e) IR1-284, age $549_{-113}^{+226}$ Myr; and (f) IR1-1169, age $4.89_{-2.26}^{+0.72}$ Gyr. Note that clusters (a) and (e) were not included in the catalog compiled by Bica et al. (2008).</figcaption></figure>

<table><caption>**Table 2**

Star Cluster Catalog</caption><thead><tr><th>CID</th><th>R.A.(J2000)<br>(deg)</th><th>Decl.(J2000)<br>(deg)</th><th>Radius<br>(deg)</th><th>log(Age)<br>(years)</th><th>Lower unc.<br>(years)</th><th>Upper unc.<br>(years)</th></tr></thead><tbody><tr><td>IR1-1</td><td>76.6915</td><td>-72.3672</td><td>0.0086</td><td>8.92</td><td>8.83</td><td>8.95</td></tr><tr><td>M2-1065</td><td>78.8299</td><td>-68.7034</td><td>0.0094</td><td>8.02</td><td>7.63</td><td>8.12</td></tr><tr><td>NUV-10</td><td>75.0914</td><td>-71.5908</td><td>0.0097</td><td>8.14</td><td>8.01</td><td>8.22</td></tr></tbody></table>

**Note.** The lower and upper uncertainty bounds are calculated as the 16th and 84th percentiles, respectively. (This table is available in its entirety in machine-readable form.)
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Distribution of star cluster radii in our sample (solid blue line), compared with that from Bica et al. (2008, dashed red line). The dashed-dotted green line is also our distribution of radii, albeit derived in an independent way, using the furthest cluster star with $p_{\text{memb}} > 0.60$.</figcaption></figure>

In particular, we adopted the following formula:

$$p_{\text{memb}} = 1 - \frac{N_{*, \text{field}}}{N_{*, \text{cluster+field}}} \cdot \frac{A_{\text{cluster}}}{A_{\text{field}}} \quad (1)$$

where $N_{*,field}$ is the number of field stars (in the 0.16 deg² box), and $N_{*,cluster+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_{\text{field}}$ and $A_{\text{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_{\text{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_{\text{memb}} \ge 0.90$ are represented with blue circles, in green circles are those with $0.60 < p_{\text{memb}} \le 0.90$, while all other stars (having $p_{\text{memb}} \le 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_{\text{memb}} \ge 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 ($\approx 11\%$ reduction).

## 4.2. The Cluster Age Determination Algorithm

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

$$p_i = \int \rho_n(x, y) U_i(x - x_i, y - y_i) dx dy, \quad (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:

$$U_i(x - x_i, y - y_i) = \frac{1}{2\pi\sigma_{x,i}\sigma_{y,i}} e^{-\left(\frac{(x-x_i)^2}{2\sigma_{x,i}^2} + \frac{(y-y_i)^2}{2\sigma_{y,i}^2}\right)}, \quad (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:

$$L_n = \prod_{i=1}^{S} p_i \quad (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_{\text{memb}}$, described in Section 4.1, for the estimation of Equation (2) ($p'_i = p_i \cdot p_{\text{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 $\pm 0.25$ mag (this corresponds to $\approx 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 $\sim 5.5$.

The isochrones we used are a byproduct of an independent project by S. Charlot & G. Bruzual (2017, in preparation).<fnref n="8" /> 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

<footnote n="8">⁸ The Charlot & Bruzual isochrones are available to the interested user upon request.</footnote>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Comparison of the ages determined from our method (Age\_current) for clusters we have in common with (a) Glatt et al. (2010), (b) Popescu et al. (2012), (c) Piatti et al. (2015), (d) Baumgardt et al. (2013), and (e) Asa'd et al. (2016). The dashed black lines correspond to the one-to-one correlation, while the dotted red ones are the least square fits to the data. The Pearson correlation coefficients (R) are indicated in the upper left corner of each panel.</figcaption></figure>

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<sub>eff</sub>, 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.
<figure><img src="image_14.png" /><figcaption>Figure 8. Age distribution of our LMC star clusters.</figcaption></figure>

### 4.3. Final Star Cluster Catalog

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:

$$ \text{Best Age} = \frac{\sum_{i=1}^{3} \frac{N_i \cdot \text{Age}_i}{(\delta \text{Age}_i)}}{\sum_{i=1}^{3} \frac{N_i}{(\delta \text{Age}_i)}}, \quad (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 $\delta 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.

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.

## 5. Results

### 5.1. Comparisons of Cluster Radii with Other Surveys

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 \le 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_{\text{memb}} > 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.

### 5.2. Comparisons of Cluster Ages with Other Surveys

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$.
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Age distribution of the star clusters found in the LMC bar (left panel) and in the rest of the galaxy (right panel).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Three (top left), four (top right), and five (bottom) component mixture models (dashed green lines), and their individual constituents (solid black lines).</figcaption></figure>

### 5.3. Age Distribution of the LMC Star Clusters

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 ≈310 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
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Spatial distribution of ages for all the star clusters in our sample (in black dots). From top-left to bottom-right we present the positions of star clusters with Age ≤ 20 Myr, 20 < Age ≤ 50 Myr, 50 < Age ≤ 100 Myr, 100 < Age ≤ 250 Myr, 250 < Age ≤ 355 Myr, 355 < Age ≤ 500 Myr, 500 < Age ≤ 750 Myr, and Age > 750 Myr.</figcaption></figure>
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.

Interestingly, the LMC bar (selected loosely by us to lie between R.A. 76h–85h and decl. −70°.5 to −69°.0; 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.

Moreover, only 50% of the clusters we detect are older than ≈300 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<fnref n="9" />). 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.

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).

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.

## 6. Conclusions

We present here a new, fully automated method we have developed to detect and estimate the ages of star clusters in nearby galaxies with resolved stellar populations. The detection is performed using a simple, but robust, algorithm (complemented by Monte-Carlo test simulations) to find overdensities even in the most crowded fields. Our procedure decontaminates the cluster candidate color–magnitude diagrams from field stars and, using a Bayesian isochrone fitting code, estimates their ages. We apply our method on multiband, high-resolution archival images of the Large Magellanic Cloud and we conclude the following.

(a) Out of 4850 clusters in the 7 deg² we surveyed, 3451 have not been reported before. Of those, ~150 (3%) have ages ≤20 Myr. Young clusters can contribute unique information for determining the stellar initial mass function.

(b) The distribution of cluster radii is consistent with the expected sizes of clusters in the LMC, peaking at 5–10 pc.

(c) The results on the age distribution of our star clusters are consistent with various epochs of star cluster formation in the LMC. The most prominent occurred ~310 ago, with secondary ones 10, 160, and 500 Myr ago. All these episodes could be the result of interactions between the LMC, the SMC, and the Galaxy, as suggested by the findings of N-body simulations (e.g., Yoshizawa & Noguchi 2003; Besla et al. 2012). We also show that the age distribution of the LMC bar star clusters is different from that of the rest of the galaxy.

(d) The spatial distribution of the clusters as a function of age suggests that the youngest clusters are located outside of the LMC bar (mostly in the northern and northeastern regions of the galaxy), whereas clusters older than 750 Myr preferentially sit in the bar. The regions with the highest concentration of young clusters (<50 Myr) are those identified by Kim et al. (2003) to host HI supershells. The above results might suggest an inside-out star cluster formation scenario for the LMC, during the past 1 Gyr.

<footnote n="9">⁹ Publicly available at https://people.maths.bris.ac.uk/~mapjg/Nmix.</footnote>
The authors would like to thank Dr. A. Piatti for insightful comments that helped improve this paper, and Dr. H. Oti for helping with the data acquisition in the initial stages of this investigation. The authors wish to thank the anonymous referee for her/his thorough review and valuable comments. T.B. would like to acknowledge support from the CONACYT Research Fellowships program. We gratefully acknowledge support from the program for basic research of CONACYT through grant number 252364. G.M. acknowledges support from CONICYT, Programa de Astronomía/PCI, FONDO ALMA 2014, Proyecto No 31140024, and GA ČR under grant number 14-21373S. R.A.G.L. thanks DGAPA, UNAM, for support through the PASPA program. G.B. acknowledges support for this work from UNAM through grant PAPIIT IG100115. This research made use of TOPCAT, an interactive graphical viewer and editor for tabular data. IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

## ORCID iDs

T. Bitsakis https://orcid.org/0000-0001-5787-8242

S. Charlot https://orcid.org/0000-0003-3458-2275

D. Zaritsky https://orcid.org/0000-0002-5177-727X

## References

Asa'd, R. S., Vazdekis, A., & Zeinelabdin, S. 2016, MNRAS, 457, 2151

Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589

Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227

Baumgardt, H., Parmentier, G., Anders, P., & Grebel, E. K. 2013, MNRAS, 430, 676

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2007, ApJ, 668, 949

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2012, MNRAS, 421, 2109

Bica, E., Bonatto, C., Dutra, C. M., & Santos, J. F. C. 2008, MNRAS, 389, 678

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Charlot, S., & Bruzual, A. G. 1991, ApJ, 367, 126

Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068

Cioni, M.-R. L. 2009, A&A, 506, 1137

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Gieles, M., Larsen, S. S., Bastian, N., & Stein, I. T. 2006, A&A, 450, 129

Girardi, L., Chiosi, C., Bertelli, G., & Bressan, A. 1995, A&A, 298, 87

Glatt, K., Grebel, E. K., & Koch, A. 2010, A&A, 517, A50

Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223

González-Lópezlira, R. A., Bruzual-A., G., Charlot, S., Ballesteros-Paredes, J., & Loinard, L. 2010, MNRAS, 403, 1213

Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757

Harris, J., & Zaritsky, D. 2009, AJ, 138, 1243

Hernandez, X., & Valls-Gabaud, D. 2008, MNRAS, 383, 1603

Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J., & Alcock, C. 2013, ApJ, 764, 161

Keller, S. C., & Wood, P. R. 2006, ApJ, 642, 834

Kim, S., Staveley-Smith, L., Dopita, M. A., et al. 2003, ApJS, 148, 473

Lejeune, T., Cuisinier, F., & Buser, R. 1997, A&AS, 125, 229

Lejeune, T., & Schaerer, D. 2001, A&A, 366, 538

Maragkoudakis, A., Zezas, A., Ashby, M. L. N., & Willner, S. P. 2017, MNRAS, 466, 1192

Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488

Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJL, 619, L1

Meixner, M., Gordon, K. D., Indebetouw, R., et al. 2006, AJ, 132, 2268

Mighell, K. J., Rich, R. M., Shara, M., & Fall, S. M. 1996, AJ, 111, 2314

Nayak, P. K., Subramaniam, A., Choudhury, S., Indu, G., & Sagar, R. 2016, MNRAS, 463, 1446

Pfalzner, S. 2009, A&A, 498, L37

Piatti, A. E., & Bica, E. 2012, MNRAS, 425, 3085

Piatti, A. E., de Grijs, R., Ripepi, V., et al. 2015, MNRAS, 454, 839

Piatti, A. E., Guandalini, R., Ivanov, V. D., et al. 2014, A&A, 570, A74

Pietrzynski, G., & Udalski, A. 2000, AcA, 50, 337

Popescu, B., Hanson, M. M., & Elmegreen, B. G. 2012, ApJ, 751, 122

Richardson, S., & Green, P. J. 1997, J. Royal Stat. Soc., 59, 731

Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25

Romaniello, M., Primas, F., Mottini, M., et al. 2005, A&A, 429, L37

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schmeja, S. 2011, AN, 332, 172

Siegel, M. H., Porterfield, B. L., Linevsky, J. S., et al. 2014, AJ, 148, 131

Simons, R., Thilker, D., Bianchi, L., & Wyder, T. 2014, AdSpR, 53, 939

Stetson, P. B. 1987, PASP, 99, 191

Subramanian, S., & Subramaniam, A. 2009, A&A, 496, 399

Vidal-García, A., Charlot, S., Bruzual, G., & Hubeny, I. 2017, MNRAS, 470, 3532

Vlajić, M., Bland-Hawthorn, J., & Freeman, K. C. 2009, ApJ, 697, 361

Walker, A. R. 2012, Ap&SS, 341, 43

Walmswell, J. J., Eldridge, J. J., Brewer, B. J., & Tout, C. A. 2013, MNRAS, 435, 2171

Werchan, F., & Zaritsky, D. 2011, AJ, 142, 48

Westera, P., Lejeune, T., Buser, R., Cuisinier, F., & Bruzual, G. 2002, A&A, 381, 524

Wofford, A., Charlot, S., Bruzual, G., et al. 2016, MNRAS, 457, 4296

Yoshizawa, A. M., & Noguchi, M. 2003, MNRAS, 339, 1135

Zaritsky, D., Harris, J., Thompson, I. B., & Grebel, E. K. 2004, AJ, 128, 1606