THE CHEMICAL SIGNATURES OF THE FIRST STAR CLUSTERS IN THE UNIVERSE

The Galactic stellar halo is a vast ancient repository that takes us back to a time when dark matter collapsed into bound structures and the Galaxy was seeded for the first time. Helmi (2008) and Tolstoy et al. (2009) provide excellent overviews of the many stellar systems and fragments that inhabit the halo. These include field halo stars (Christlieb et al. 2002; Cayrel et al. 2004; Frebel et al. 2005; Cohen et al. 2007), globular clusters (Gratton et al. 2004), dwarf spheroidals (Mateo 1998; Venn et al. 2004), ultra-faint dwarf galaxies (Simon & Geha 2007; Kirby et al. 2008), stellar streams (Ibata et al. 1995; Chou et al. 2010), stellar associations (Walsh et al. 2007), and satellites to dwarf galaxies (Coleman et al. 2004; Belokurov et al. 2009). Already the chemical information arising from the most metal-poor stars is very difficult to unravel (Nomoto et al. 2005; Kirby et al. 2008). We have barely begun to understand what these systems are telling us about the sequence of events that led to the Galaxy (McWilliam et al. 2009; Freeman & Bland-Hawthorn 2002).

The first stars were unique to their time. The first stellar generations changed the universe in many ways; in particular, the chemical properties and the equation of state of the intergalactic medium. But at present there are many unknowns. Did the first stars form in isolation or in groups? Were relatively few massive stars responsible for reionization or was it triggered by the collective effect of massive star clusters? Just what are the processes that govern star formation at extremely low metallicity? Is this exclusively the domain of the most massive stars, or can substantial intermediate and low mass stars form (Tsuribe & Omukai 2006, 2008; Clark et al. 2008)? In other words, did stellar populations observable today exist before reionization (Tumlinson 2010; Okrochov & Tumlinson 2010)? The first star clusters are of great interest in that they shed light on star formation processes in the early universe (Bromm et al. 2002; Abel et al. 2002). There are few if any reliable constraints at the present time.

One of the most interesting developments in recent years is the simultaneous measurement of many elemental abundances for individual metal-poor halo stars or groups of stars (e.g., Beers & Christlieb 2005). Some of these elements, but not all, exhibit intrinsic scatter in the abundance plane [Fe/H] versus [X/Fe] that exceeds the measurement errors. While the scatter is particularly apparent in halo stars (e.g., Roederer et al. 2009), evidence is now emerging that star-to-star abundance variations exist in dwarf galaxies as well (Fulbright et al. 2004; Koch et al. 2008; Feltzing et al. 2009). The observed scatter is likely to increase now that metal poor stars are now detected below [Fe/H] = −3 (Norris et al. 2010; Starkenburg et al. 2010; Simon et al. 2010; Frebel et al. 2010). This has led numerous researchers to argue that the scatter in [X/Fe] is a tracer of an ancient inhomogeneous medium (e.g., Audouze & Silk 1995; Ryan et al. 1996; McWilliam 1997; Ishimaru & Wanajo 1999). If this interpretation is correct, we would expect that a subset of these stars is telling us something fundamental about the first stars and their yields. But not all of the stars are providing us with an unambiguous "reading" of the early enrichment of the primordial interstellar gas. For example, a significant fraction of extremely metal poor stars appear to have undergone mass transfer with a companion (Ryan et al. 2005; Lucatello et al. 2005) which undermines any attempt at inferring the progenitor yields for elements such as CNO and α elements. Binarity may partly explain why the elemental abundances of the most metal poor stars defy a clear explanation at the present time (Joggerst et al. 2010; McWilliam et al. 2009).

We now incorporate a missing ingredient into existing models of stochastic chemical evolution. In the present-day universe, most stars are born in a single burst within compact clusters and stellar fragments, rather than in isolation. This fact is well established in the local universe (Lada & Lada 2003) and is likely to be true at the time of the first stars (Clark et al. 2008). Cluster formation has an important consequence for the distribution of stars in the abundance plane. De Silva et al. (2006, 2007a,b) have shown that both old (∼10 Gyr) and intermediate-age (∼1 Gyr) open clusters are chemically homogeneous to a high degree (Δ[Fe/H] ≲0.03 dex). The open cluster Tombaugh 2 was thought to be a rare example of a chemically inhomogeneous open cluster (Frinchaboy et al. 2008), but this is contradicted by a new study that finds the stellar population to be highly homogeneous (Δ[Fe/H] ≲0.02 dex; Villanova et al. 2010). Apart from a few light elements, the same holds true for globular clusters (Gratton et al. 2004), although some systems (e.g., ωCen) show evidence for more than one burst of star formation. Interestingly, even moving groups can show the same signature of chemical homogeneity (De Silva et al. 2007a; Chou et al. 2010; Bubar & King 2010).

Bland-Hawthorn et al. (2010) provide a condition that ensures chemical homogeneity in a young star cluster based on the surface density of the progenitor gas. They show that chemical homogeneity is expected in open clusters, globular clusters and plausibly clusters more massive than 10⁷ M_☉. This process has not been incorporated into stochastic chemical models to date. Once this is done, we arrive at an important insight that will enhance the interpretation of metal abundance distributions in individual stellar populations. The effects that we highlight can be searched for in existing and in future surveys, as we show.

A key assumption in our present work is that present-day dwarf galaxies are important sites for establishing the yields of the first stars and star clusters. This needs some clarification. While it is likely to be true that the most efficient way to identify metal-poor stars is to target dwarf galaxies, a low value of [Fe/H] is no guarantee that a star is ancient since it may simply reflect environmental conditions (e.g., low star formation efficiency, weak gravity field). Published numerical simulations are unclear on whether the most ancient stars are solely the preserve of the inner bulge (White & Springel 2000; Bland-Hawthorn & Peebles 2006) or spread over the entire Galaxy (Scannapieco et al. 2006; Brook et al. 2007). But it is now well established that the bulge, the halo and all dwarf galaxies comprise stellar populations that are 10 Gyr or older, equivalent to a formation redshift of z ≳ 2 (e.g., Tolstoy et al. 2009). The relative fractions of dwarf populations that formed before, during, or after the reionization epoch is an open question. We are therefore at liberty to explore new observational constraints on the nature of the first stars and star clusters.

In Section 2, we present evidence for homogeneous star clusters in the local universe and argue that the same should be true for their high-redshift counterparts. In Section 3, we briefly outline the inhomogeneous chemical evolution models that have been developed to date. In Section 4, we introduce a revised stochastic model of star formation in dwarf galaxies that incorporates the "homogeneity" condition described by Bland-Hawthorn et al. (2010), and show the predictions of the revised model in view of the earlier work. In Section 5 and in the Appendix, we introduce cluster finding algorithms for both large and small data sets, and demonstrate how the effects of clustering may already be visible in existing observations. We also explore the longer term prospect offered by a multi-object echelle spectrograph on an extremely large telescope (ELT). The conclusions are presented in Section 6.

Before revisiting stochastic chemical evolution models, we review the key arguments presented in Bland-Hawthorn et al. (2010) on the conditions under which star clusters are expected to be highly homogeneous in most elements, as observed for local star clusters (Section 1). We then extend these arguments to clusters at low metallicity.

As a star-forming cloud assembles, turbulent diffusion within it will homogenize its chemical composition (Murray & Lin 1990). For clouds whose turbulent motions are primarily on large scales (as is the case for all local molecular clouds—Heyer & Brunt 2004), Bland-Hawthorn et al. (2010) show that the time required for this process to smooth out a composition gradient on the size scale of the cloud is roughly t_cr = L/σ, where t_cr is the cloud crossing time, L is its size, and σ is its velocity dispersion. Smaller-scale gradients are erased more quickly, with the diffusion time varying as the square of the (normalized) characteristic size scale. Since t_cr is comparable to or smaller than the time scale over which star formation takes place, clouds will homogenize as they assemble, and will be pushed away from homogeneity only if supernovae (SNe) occur during the star formation process, before the cloud disperses. The time required for a very massive star to evolve from formation to explosion defines the SN time scale, t_SN ≈ 3 Myr, and we only expect star clusters to be homogeneous if they are assembled on time scales shorter than t_SN.

To determine under what conditions this requirement is satisfied, we must compare t_SN to the cluster formation time scale t_form. There is considerable debate over this time scale, so for our purposes we will adopt the longest, most conservative proposed time scale of 4 t_cr (Tan et al. 2006). Since

$$\begin{equation} t_{\rm cr}= \frac{0.95}{\sqrt{\alpha _{\rm vir} G}}\left(\frac{M}{\Sigma ^3}\right)^{1/4} \end{equation} \tag{ 1 }$$

for a cloud of mass M, column density Σ, and virial ratio α_vir; observed clouds have α_vir ≈ 1.5 (McKee & Ostriker 2007). For convenience, we write the final stellar mass of a cluster as M_* = M, where is the star formation efficiency. Both observational and theoretical arguments suggest ≈ 0.2 independent of M (Lada & Lada 2003; Fall et al. 2010). Thus the condition that t_form ≈ 4t_cr < 2t_SN (where the factor of 2 arises because the typical star forms halfway through the formation process) is satisfied only if

$$\begin{equation} M_{*,5}^{1/4} \Sigma _0^{-3/4} < 2.8, \end{equation} \tag{ 2 }$$

where M_*,5 = M_*/(10⁵ M_☉) and Σ₀ = Σ/(1gcm⁻²), and we have adopted fiducial values of = 0.2 and α_vir = 1.5. Star-forming regions within the galaxy have Σ₀ ≈ 1 g cm⁻² independent of mass (Fall et al. 2010). Globular clusters today have somewhat higher values of Σ, although it is unclear if this reflects the conditions under which they formed, or is the result of dynamical evolution since their formation. Regardless, this analysis suggests that clusters with masses up to ∼10⁵ M_☉ should be chemically uniform and maybe much higher if systems form with the nuclear densities observed in today's globular clusters (Bland-Hawthorn et al. 2010).

There are few observational constraints on the existence of star clusters in extremely metal poor gas. In the nearby universe, clear evidence for star clusters at [Fe/H]≈  −1.7 is observed in the most metal-poor, blue compact dwarf galaxy I Zw 18 (Izotov & Thuan 2004). Globular clusters are known to exist down to [Fe/H]≈  −2.4 (Gratton et al. 2004). We now argue that the above estimate for the uniform mass limit is likely to apply at metallicities as low as [Fe/H] ≈   −5.

Low metallicity can change the star formation process in two ways that are relevant to the chemical signatures of the resulting stars. First, a change in metallicity can affect the way star-forming clouds fragment; if there is no fragmentation down to sub-solar masses below a certain metallicity, then no stars below that metallicity will survive to the present day. Fragmentation of low metallicity gas has received extensive attention in the literature, which we will only summarize here. The main point relevant for our purposes is that, when dust cooling is considered, gas is able to fragment to sub-solar masses even at metallicities at low as ∼10⁻⁶ of the solar abundance (Clark et al. 2008; Schneider & Omukai 2010). Once fragmentation is possible, turbulence naturally generates a mass spectrum of fragments with a slope comparable to the Salpeter slope dn/dm ∝ m^−2.35 (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008), and the properties of the turbulence do not depend on the redshift or the metallicity. Alternately, Clark et al. (2009) have proposed that competitive accretion processes would generate a universal mass spectrum, although doubts have been raised about whether this process in fact operates on both theoretical (Krumholz et al. 2005) and observational (André et al. 2007) grounds. In either case, however, we expect there to be some sub-solar mass stars formed that can survive to the present day.

Given that small stars form at low metallicity, we can then ask whether our conclusions about cluster chemical homogeneity will continue to apply in this regime. A failure of homogeneity could occur either if clouds did not homogenize during star cluster formation, or if the cloud properties changed such that the cluster formation time became longer than the SN time scale. The former is unlikely because the homogenization time is comparable to the crossing time, and it is implausible that any star formation process could take place in less time. The latter would require that protocluster gas clouds all have surface densities significantly below ∼0.1 g cm⁻² which is highly unlikely to be true. However this would imply that the star-forming clouds had surface densities at or below the mean surface density of observed high redshift galaxies (e.g., Genzel et al. 2006), which is implausible. We therefore conclude that homogeneity should continue to apply as well.⁶

3.1. A Statistical Treatment

Many authors have attempted to explain the declining scatter in elemental abundances with increasing [Fe/H] in terms of star formation within an interstellar medium (ISM) that is enriched by a succession of SN events (Ishimaru & Wanajo 1999; Shigeyama & Tsujimoto 1998; Raiteri et al. 1999; Argast et al. 2000, 2004; Wasserburg & Qian 2000; Karlsson & Gustafsson 2001, 2005; Qian 2000, 2001; Travaglio et al. 2001; Fields et al. 2002). In the [Fe/H] versus [X/Fe] plane,⁷ the scatter in [X/Fe] is due to two or more classes of SN that produce distinct yields of [X/Fe]. In these stochastic models, the scatter converges roughly quadratically to a mean value of [X/Fe] given by the SN yields weighted by the initial mass function (IMF). In its simplest form, this behavior is driven by the number of SN events (n_SN), such that [Fe/H] = [Fe/H]_min + log n_SN, where [Fe/H]_min is the minimum allowed metallicity. In other words, if [Fe/H]_min = −5 consistent with the most metal-poor stars to date (e.g., Frebel et al. 2010), it takes 10⁵ SNe to enrich a gas parcel to solar abundance. An example of this behavior is presented in Bland-Hawthorn & Freeman (2004; Figure 3) where the number of SNe are indicated.

If the interpretation of declining inhomogeneity is broadly correct, we can hope to learn about the yields of the first SNe responsible for the enrichment (Fields et al. 2002). In other words, if the scatter in [X/Fe] is due to a high-yield source (class A) and a low-yield source (class B), an individual star with known [Fe/H] has a fraction of f_A = n_A/n_SN progenitors and f_B = n_B/n_SN progenitors of a total progenitor population of n_SN = n_A + n_B SNe. We can simulate this with a random variable f = f_A drawn from the beta distribution function,

$$\begin{equation} \xi (f;\alpha,\beta) = {{\Gamma (\alpha +\beta)}\over {\Gamma (\alpha)\Gamma (\beta)}}f^{\alpha -1}(1-f)^{\beta -1} \end{equation} \tag{ 3 }$$

defined for f ∈ [0, 1] such that ∫¹₀ξ df = 1 and where Γ is the gamma function. The quantities α and β describe the yields of the two populations, such that $\alpha =\bar{f_A}(n_{\rm SN}$-1) and $\beta = (1-\bar{f_A})(n_{\rm SN}$-1). The mean, variance and higher moments of the normalized distribution ξ(f) depend only on α and β.

Equation (3) has several remarkable properties that are useful for describing the declining influence of chemical inhomogeneities in the early ISM. The mean is given by $\mu = \alpha (\alpha +\beta)^{-1} = \bar{f_A}$ such that $\bar{f_A}$ is the mean of the distribution ξ(f_A). The variance is given by σ² = αβ(α + β)⁻²(α + β + 1)⁻¹ which leads to a scatter that declines as σ ∝ n^−0.5_SN as expected.⁸ Even initially skewed distributions in ξ converge to the normal distribution ($\mu = \bar{f_A}$) in the limit of high n_SN, as expected from the central limit theorem.

Four realizations of the beta distribution (i.e., no clustering) are presented in Figure 1; a total of n_* = 3000 points was used for each simulation distributed in [Fe/H] according to Figure 2 (see Section 5). The log ξ distribution can be renormalized trivially to match a given set of [Fe/H] versus [X/Fe] observations (e.g., Fields et al. 2002). The beta distribution ξ is everywhere continuous in f_A, and therefore f_B, except in the limit of small n_A or n_B. This recognizes the fact that, as stated by Fields et al. (2002), "a given parcel of ISM gas and dust can be enriched to different degrees by the ancestors it had." But in the limit of small n_A or n_B, discreteness effects may exist. Stars with very few prior enrichments, if they can be identified, provide crucial information on the progenitor yields of class A or B sources (e.g., Shigeyama & Tsujimoto 1998; Karlsson & Gustafsson 2001; Ballero et al. 2006), but within the domain of the statistical model, these are rare events and may be difficult to find unambiguously. This is an issue we return to in Section 5. At high [Fe/H], when new sources of metals become dominant (e.g., type Ia SNe, asymptotic giant branch stars [AGBs]), the simple description in Equation (3) breaks down.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To recap, the model described by Equation (3) assumes that a star with a given [Fe/H] has been enriched to different degrees by its type A and type B ancestors which is considered to be a continuous rather than a discrete process. Equation (3) is used to generate a distribution of possible values of f_A at a constant value of [Fe/H] or equivalently n_SN. More complex models show that there is generally no simple relation between f_A and [Fe/H] (e.g., Qian 2001; Bland-Hawthorn & Freeman 2004).

We note already from Figure 1 that spurious groupings (false positives) are inevitable, particularly when the data points are convolved with typical errors of 0.1 dex psf for differential abundance analysis. The importance of the beta distribution in Equation (3) is that it enables us to efficiently generate large numbers of unbiased realizations (≳10³) of the theoretical abundance plane, something that is infeasible with full-blown chemical modeling. These are required to calibrate the statistical properties of our group-finding algorithm in Section 5. This step is necessary if we are to identify significant groups against a rapidly changing background, as observed in Figure 1. Clusters at low [Fe/H] will have a higher level of significance than abundance groupings at higher [Fe/H] (cf. Karlsson et al. 2008) and we need to be able to identify groups with smaller number statistics.

But first we must consider more complex treatments of stochastic chemical evolution that will allow us to include the effects of clustering during star formation.

3.2. Inhomogeneous Stochastic Model

4.1. Initial Cluster Mass Function

In order to derive the impact of star clusters on the abundance plane, we must consider the progenitor mass distribution of star clusters. It is now well established that star clusters have a range of masses that extend from a maximum mass (M_max) to a minimum mass (M_min). The form of the birth distribution is known as the initial cluster mass function (ICMF) and is assumed to have the form

$$\begin{equation} {d}N/{d}M_{*} = \chi (M_{*}) = \chi _0 M_{*}^{-\gamma }. \end{equation} \tag{ 4 }$$

The observations may support a universal slope of γ ≈ 2 in all environments, i.e., equal mass per logarithmic bin (Fall et al. 2005, 2009; Lada & Lada 2003; Elmegreen 2010). In this picture, the only parameter that does appear to vary at all between galaxies today is M_max, which can be represented schematically as a Schechter function-like cutoff in the ICMF, although there is no good reason to prefer this functional form over a simple truncation.

A power law is logarithmically divergent at both low and high masses, so it must be truncated somewhere. At low masses, the truncation is due to the discreteness of stellar masses—the smallest cluster is simply one star. At large masses, there must also be a truncation. For this reason, and because the molecular cloud mass function (as opposed to the cluster mass function) is observed to have a non-trivial truncation, there is likely to be a maximum cluster mass, and that it varies depending on galactic environment. In disk galaxies, a possible explanation is that the truncation mass is of order the Toomre mass in the galactic disk (Toomre 1964; Escala & Larson 2008), which is ∼10⁶ M_☉ for the Milky Way, but is significantly larger for present-day starburst/merger galaxies (Larsen 2009) and their high-redshift counterparts (Genzel et al. 2006; Förster-Schreiber et al. 2009). The situation in spheroidal galaxies is much less clear but the common occurrence of globular clusters indicates that massive cluster formation must take place (Elmegreen 2010).

4.2. Stochastic Chemical Evolution

We take as our starting point the stochastic chemical enrichment model presented in Karlsson (2005, 2006) and Karlsson et al. (2008). This model is now updated to include the homogenizing processes that must occur during the formation of star clusters (Section 2). But in order to do this, two additional mixing processes must be accounted for. First, since stars within clusters show no evidence of scatter in chemical abundance ratios, the gas involved in the formation of a cluster must be homogenized prior to, and stay homogenized during, the formation of the individual stars. We include a mixing process to ensure that this is always true during the cloud collapse and cluster formation phase. Second, since massive stars are short lived (τ_⋆ ≲ 20 Myr) they will explode as SNe before the cluster has dispersed. Assuming that the cluster is unbound and that the intrinsic velocity dispersion of the stars in the cluster is ∼1 km s⁻¹, massive stars will, on average, be ∼20 pc apart when they go off as SNe. This is less than the typical size (∼100 pc) of an SNR as it merges with the ambient medium. (e.g., Ryan et al. 1996). Therefore, we assume that the ejecta of all SNe formed within a single cluster enrich the ISM collectively with a mixing mass that scales linearly with the total energy output of the clustered SNe. These two processes, which produce additional averaging of newly synthesized material before a new generation of stars is formed, have not been considered in earlier models.

The mixing volume, V^e_mix, of each chemical enrichment event is given as a power-law expression. The present model is somewhat simplified as we suppress the continuous mixing due to the bulk random motions in the turbulent ISM. Without the time-dependent turbulent diffusion term (cf. Karlsson et al. 2008, Equation (1)), V_mix reduces to

$$\begin{equation} V_{\mathrm{mix}}^{\mathrm{e}}(k) = \frac{4 \pi }{3}\sigma _{\mathrm{E}}(k)^{3/2} = M_{\mathrm{dil}}(k)/\rho, \end{equation} \tag{ 5 }$$

where

$$\begin{equation} \sigma _{\mathrm{E}}(k) = \left(\frac{3M_{\mathrm{dil}}(k)}{4 \pi \rho }\right)^{2/3}, \end{equation} \tag{ 6 }$$

ρ is the density of the ISM and M_dil(k) is the dilution mass of the ejecta of a number k (≡n_SN) of SNe as they merge with the ambient medium (e.g., Cioffi et al. 1988). Here, we assume that the total dilution mass for the ejecta of multiple SNe exploding in a cluster is proportional, on average, to the total energy released by the SNe associated with that cluster, such that

$$\begin{equation} M_{\mathrm{dil}}(k) = \sum _{j=0}^{k}M_{\mathrm{sw}}^j, \end{equation} \tag{ 7 }$$

where M^j_sw is the mass swept up by individual SNe assuming it explodes in isolation. In order to introduce a small amount of randomness, log(M_sw) is drawn from a normal distribution centred on log(M_sw) = 5 with a width of 0.25. Note that this mixing process is denoted a single enrichment event, even though it may involve enrichment by multiple SNe.

In addition to this mixing, there is also a mixing volume associated with the formation of the cluster during the collapse of the molecular cloud, V^f_mix, such that

$$\begin{equation} V_{\mathrm{mix}}^{\mathrm{f}} = M/\rho, \end{equation} \tag{ 8 }$$

where M is the mass of the molecular cloud associated with the star cluster. Within this volume, everything is assumed to be thoroughly mixed before stars are allowed to form.

In our models, we explore the range 1 ≲ γ ≲ 2.5 (e.g., Kroupa & Boily 2002), between clusters of mass (M_min, M_max) = (5 M_☉, 5 × 10⁴ M_☉), which is our fiducial mass range (Larsen 2009; Portegies Zwart et al. 2010). We adopt a simple scaling relation between the cluster mass and the mass of the parent molecular cloud, such that M = M_*/, where we adopt a star formation efficiency of = 0.2 (see Section 2). Similarly, the stellar initial mass function (IMF) is governed by

$$\begin{equation} {d}n/{d}m = \phi (m) = \phi _0m^{-\alpha }. \end{equation} \tag{ 9 }$$

For simplicity, α = 2.35, which recovers the Salpeter IMF. The range of stellar masses is set to 0.1 ⩽ m/M_☉ ⩽ 100 and ϕ₀ = 6.03 × 10⁻² M_☉^1.35.

Formally, the average number of clusters contributing to the chemical enrichment in a random point in space can be expressed by the parameter μ_e, here given by

$$\begin{equation} \mu _{\mathrm{e}}(t) = \int _0^t \sum _{k=0}^{k_{\mathrm{max}}} a_k\times V_{\mathrm{mix}}^{\mathrm{e}}(k) u_{\mathrm{cl}}(t^{\prime }){d}t^{\prime } \end{equation} \tag{ 10 }$$

where $a_k = N_k/\sum _{k=0}^{k_{\mathrm{max}}}N_k$ is the fraction of clusters in which k SNe explode, while u_cl(t') is the formation rate of star clusters, closely related to the star formation rate and V^e_mix(k = 0) = 0. The value of k_max is set by the IMF and the upper mass limit of star clusters. On average, 270 massive stars explode as SNe in a cluster of mass 5 × 10⁴ M_☉. Such a cluster will not form more than 320 SNe at the 3σ level. We set k_max = 400 to be on the safe side.

Now, if we make the simplifying assumption that star clusters are randomly distributed in space (i.e., not themselves clustered), the probability of finding a region in space enriched by κ events, i.e., κ clusters producing one or more SNe, at time t is given by the Poisson distribution

$$\begin{equation} P(\kappa,\mu _{\mathrm{e}}(t))=e^{-\mu _{\mathrm{e}}(t)}\mu _{\mathrm{e}}(t)^{\kappa }/\kappa ! \end{equation} \tag{ 11 }$$

To follow the chemical enrichment in a typical dwarf galaxy, we assume a simulation box of 1.6 kpc on a side, with an initial, constant particle density of n₀ = 1 cm⁻³, corresponding to a gas density of ρ₀ = 2.06 × 10⁻²⁴ g cm⁻³. The initial mass of baryons in the box is thus 3.1 × 10⁷ M_☉ which is sufficient to make the largest star clusters considered in this work. The initial metallicity is set to Z = 0. In the box, star clusters are allowed to form from collapsing molecular clouds. The molecular clouds are distributed randomly within the box and the masses of the corresponding clusters are distributed according to Equation (4). During the collapse of a molecular cloud, the volume, V^f_mix, of gas corresponding to the mass of the cloud is made chemically homogeneous. The stars of the cluster will all have chemical abundances equal to those of the parent cloud. The number of massive stars, k, in the cluster that will explode as SNe is, again, determined by the Poisson statistics P(k, μ_SN), where the mean number of SNe in a cluster of a given mass is given by

$$\begin{equation} \mu _{\mathrm{SN}} = \epsilon M f_{\mathrm{SN}}/\overline{m}, \end{equation} \tag{ 12 }$$

where is the star formation efficiency, M is the mass of the molecular cloud, f_SN = 1.9 × 10⁻³ is the fraction of massive stars exploding as SNe, and $\overline{m}=0.35 \,{M_{\odot }}$ is the mean stellar mass in a single stellar population. The last two parameters are determined by the IMF. For large values of μ_SN, the Poisson distribution approaches a normal distribution. We make use of this fact to simplify the calculations for μ_SN>32. In clusters where one or more SNe are formed, the ejecta of all SNe, with masses distributed according to Equation (9), are homogeneously mixed with the metals already present within the volume V^e_mix. The Fe-core collapse SN yields, in particular of Ca and Fe, are taken from Nomoto et al. (2006) while the yields of Eu, representing r-process elements and assumed to be formed in O–Ne–Mg core collapse SNe in the mass range 8 ⩽ m/M_☉ ⩽ 10, are taken from Argast et al. (2004). In order to match the simulations with observations of metal-poor stars in the Galactic halo, the yields of Eu are doubled.

In order to account for the amount of mass locked up in low-mass stars and stellar remnants and the mass lost due to star formation driven galactic outflows, a fraction, f_out = 0.46, of the mass of each star-forming molecular cloud is subtracted from the total mass of the system. The gas density is decreased accordingly. Infall and mass lost via interaction with the surroundings, such as tidal and ram pressure stripping are not considered here. When the total mass and gas density of the system has been re-calculated, the next cycle begins with the formation of a new cluster. Since our focus is in the metal-poor regime, the yields of type Ia SNe and AGB stars are not considered here.

4.3. Results

The results of the clustered simulations for [Ca/Fe] and [r/Fe] with respect to [Fe/H] are presented in Figures 3 and 4, respectively. For all models, the number of simulated stars is n_* ∼ 10⁷, equivalent to a stellar mass of 3 × 10⁶ M_☉ and a luminosity of 10⁶ L_☉ assuming a Salpeter IMF; the luminosity is two times higher for a Kroupa IMF. To begin with, we do not consider the effect of distance on the observed number counts; this is treated in Section 5.2. The modeled number of stars is much larger than we can expect to obtain in nearby dwarfs, but the results provide sufficient resolution to understand the impact of the ICMF parameters. The overall distribution of the stars with metallicity is shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is a clear progression moving from high γ to low γ in the occurrence of clustered abundance signatures. The high γ limit is indistinguishable from the beta (non-clustered) models. This is not unexpected since, in the limit of γ = 2.5, almost all stars are formed in small star clusters such that unique groupings in abundance space are poorly represented.

As γ decreases, the clustering increases markedly but at the expense of the dispersion in [X/Fe]. In the limit of γ = 1, the vertical spread has vanished over all [Fe/H] with little or no abundance spread in [X/Fe]. While γ = 1 is smaller than what is observed, it is in fact a useful surrogate for demonstrating the impact of a higher maximum cloud mass M_max for larger γ values. For example, in Figure 5, we show the result of running γ = 1.5 and γ = 2.0 models with M_max ≈ 10⁶ M_☉. We see how the high mass cutoff at high γ mimics the behavior of a lower value of γ: the vertical scatter is greatly reduced, and groupings are seen to extend to lower [Fe/H].

Zoom In Zoom Out Reset image size

In Figure 6, we show the cumulative fraction of stars and the integrated light as a function of stellar absolute magnitude for the model dwarf galaxy. In Table 1, we give a rough breakdown of how many stars are expected for the dwarf galaxy as a function of stellar apparent magnitude, distance and metallicity. We consider V = 16 and V = 18 to be the bright and faint limit of an 8 m class experiment; these magnitude brackets become V = 20 and V = 22 on an ELT (Bland-Hawthorn et al. 2010).

Zoom In Zoom Out Reset image size

Table 1. Log of the Expected Numbers of Stars as a Function of V mag for an Old, Metal Poor Dwarf Spheroidal (see Figure 2) with 10⁷ Stars (Stellar Mass ≈ 3 × 10⁶ M_☉) at Distances of 10, 30 and 100 kpc (Column 1)

[Fe/H]	Salpeter				Kroupa
	−5:−4	−4:−3	−3:−2	−2:−1	−5:−4	−4:−3	−3:−2	−2:−1
V = 16
10	0.7	1.9	3.0	3.5	1.0	2.2	3.3	3.8
30	⋅⋅⋅	0.8	1.9	2.4	⋅⋅⋅	1.1	2.2	2.7
100	⋅⋅⋅	-	0.4	1.0	⋅⋅⋅	-	0.7	1.3
V = 18
10	1.1	2.3	3.4	3.9	1.4	2.6	3.7	4.2
30	0.7	1.9	3.0	3.5	0.9	2.2	3.2	3.8
100	⋅⋅⋅	0.6	1.7	2.2	⋅⋅⋅	0.9	2.0	2.5
V = 20
10	2.2	3.4	4.5	5.0	2.5	3.7	4.8	5.3
30	1.0	2.2	3.3	3.8	1.3	2.5	3.6	4.1
100	0.1	1.3	2.4	2.9	0.4	1.6	2.7	3.2
V = 22
10	2.8	4.0	5.1	5.6	3.0	4.3	5.3	5.8
30	2.0	3.3	4.3	4.8	2.3	3.5	4.6	5.1
100	0.9	2.1	3.2	3.7	1.2	2.4	3.4	4.0

Notes. Columns 2–5 are the star counts for a Salpeter IMF in four metallicity bins; columns 6–9 are the counts for a Kroupa IMF in the same bins determined from the MDF in Figure 2. The four metallicity bins are [Fe/H] = (−5: − 4), (−4: − 3), (−3: − 2), (−2: − 1); dashes indicate that no stars are expected.

Download table as:  ASCIITypeset image

In the next section, we analyze these simulations using a group finder in order to provide an objective assessment of the amount of measurable clustering in the limit of high and low n_* that we can expect. We apply this analysis for different [Fe/H] cutoffs since, as we have seen, the abundance plots become crowded as [Fe/H] increases.

5.1. Group Finding for Large n_*

5.2. Group Finding for Small n_*

Clustering should be present even in the limit of only a few data points. To emphasize this fact, we have simulated the abundance measurements for a dwarf galaxy at a distance of 30 kpc (see Figures 7 and 8) as observed on 8 m class and 30 m class telescopes respectively. We adopt a stellar mass of 3 × 10⁵ M_☉ typical of a faint dwarf galaxy. This object has about 10⁶ stars, a luminosity of 10⁵L_☉ and an absolute V mag of M_V = −7.6 assuming a Salpeter IMF. The luminosity is a factor of two higher for a Kroupa IMF. The star counts are consistent with the model values in Table 1 when scaled to the adopted lower mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the 8 m experiment (Figure 7), we assume measurement errors of 0.1 dex in both [r/Fe] and [Fe/H]. For γ = 1.5 and γ = 2.0, the effects of clustering are evident and this holds true if we double the measurement error. The impact of this measurement uncertainty over the full simulation is shown in Figure 9. We conclude that, once more extensive studies are made of the nearest dwarf galaxies, the effects of clustering may become evident even before the advent of ELTs. We envisage projects which focus on relatively few stars that appear grouped in abundance space in low resolution spectroscopic data. Long integrations and differential analysis will allow the null hypothesis to be tested that the stars have identical abundances in all elements. For the 30 m experiment (Figure 8), we have assumed a general improvement in the atmospheric models and the experimental errors, and therefore adopt errors of 0.05 dex. The effects of clustering, which are easily seen, remain clearly visible even after a twofold increase in the measurement errors in both axes. This simulation is a powerful statement of the importance of multi-object echelles on ELTs.

Zoom In Zoom Out Reset image size

For a fixed q_den and data dimensionality, the number of spurious groups due to Poisson noise increases linearly with the total number of data points n_*. Therefore, EnLink can be used in the limit of small n_* with q_den set to 3, i.e., a minimum of one more than the number of data dimensions. But with so few data points, the full power of EnLink is not being exploited such that more rudimentary statistical techniques may be better. However, EnLink is particularly efficient in treating data sets with more than two dimensions, even in the limit of small n_*. The significance of groups increases dramatically when we apply EnLink on an abundance space with more dimensions. The normalizing distribution in Equation (3) is easily extended to higher dimensions. In Figure 10, we apply EnLink to the 3D ${\cal C}$-space ([Fe/H], [Ca/Fe], [r/Fe]) for three different measurement errors (0.05, 0.1, 0.2 dex) typical of contemporary observations on 8 m telescopes at V = 18. Even in the presence of large errors, clustering signals are seen in all cases.

Zoom In Zoom Out Reset image size

This paper has explored the prospect of probing the mass scales of the first star clusters. We stress that we have used a lower cluster-mass limit of M_min = 5 M_☉ (typically 10 stars) which is very conservative. A more reasonable value may be an order of magnitude higher (cf. Bland-Hawthorn et al. 2010). For γ = 2 and a fixed number of simulation particles, the lower threshold has the effect of reducing cluster membership by a factor of 10, and increasing the "background" by the same factor. This lowers the overall clustering signal by an order of magnitude. But we adopt the conservative lower mass limit in order to account for a possible dwarf stellar population that are not born in clusters. We were unable to find any observational constraints on the diffuse versus clustered population in dwarf galaxies (cf. Lada & Lada 2003).

In our conservative analysis, we find that there is an intimate connection between properties of the ICMF and the amount of potentially detectable clustering in the abundance plane. A flat ICMF and/or a ICMF with a high mass cutoff produces strong clustering in the [Fe/H] versus [X/Fe] abundance plane (${\cal C}$-space). While our models are inevitably oversimplified, the phenomenon should be detectable on 8–10 m telescopes (e.g., Figures 7 and 8). This gains support from existing observations of open clusters, globular clusters and moving groups (Castro et al. 1999; Shen et al. 2005; Randich et al. 2006; Sestito et al. 2007; De Silva et al. 2009; Chou et al. 2010; Bubar & King 2010). A clean "clustering" signature in ${\cal C}$-space, particularly at low metallicity, is important for a number of reasons. First, it indicates the presence of massive star clusters in the early universe and conceivably provides a constraint on the mass of the first systems. Second, it provides a clean signal of the progenitor abundances in the cloud prior to cluster formation. This abundance measurement is averaged over a substantial amount of gas and is therefore not subject to mixing anomalies (Karlsson & Gustafsson 2001) or mass transfer in binaries (e.g., Suda et al. 2004; Ryan et al. 2005; Lucatello et al. 2005).

Strong clustering would indicate a highly flattened ICMF, or a high mass cutoff. This could herald the onset of the formation of massive star clusters in dwarf galaxies (e.g., Bromm & Clarke 2002). If the star formation efficiencies were low at that time, this may require supermassive gas clouds (≳10⁷ M_☉) to have formed even at the earliest times (Abel et al. 2000), possibly consistent with the regular occurrence of massive star-forming clumps at high redshift (Genzel et al. 2006; Förster-Schreiber et al. 2009; Elmegreen & Elmegreen 2006). Conversely, if such clustering was not observed, then we would infer that the slope of the early ICMF is steep, or the maximum cluster size is relatively small compared to the present day. But the observed scatter would need to be consistent with the non-detection of clustering.

In future investigations, from the perspective of chemical signatures, we will look at the degeneracy between steep ICMFs with a high mass cutoff and flatter ICMFs with a lower mass cutoff (see Section 4.2). We will also look at a wider class of chemical elements that are considered in numerical simulations of the first stars. We will look at the improvement in cluster identification with more chemical elements, particularly in the limit of small n_*. The clustered abundance signatures will provide unique insight into the most ancient star clusters. These signatures are signposts of the chemistry immediately before the onset of star formation, untainted by mass transfer in close binary pairs or incomplete mixing anomalies.

Finally, it is an extraordinary fact that we can probe back to the first billion years from observations of the local universe. We can say with absolute certainty that stars existed at this time. These were responsible for the first chemical elements (Ryan-Weber et al. 2006) and for reionizing the fog of hydrogen that permeated the early Universe (Fan et al. 2002). Precisely when the first star clusters formed is unknown. It seems likely, however, that gas was able to fragment at very high density even at primordial abundance levels (Clark et al. 2008). It may be possible to directly probe these environments in an era of the Atacama Large Millimetre Array and the James Webb Space Telescope. But we believe that some of the most important insights, particularly with regard to progenitor yields, will undoubtedly come from near-field cosmology. To this end, it will be necessary to equip the next generation of ELTs with wide-field multi-object spectrographs that operate at high spectroscopic resolution (R ≳ 20,000).

J.B.H. is supported by an Federation Fellowship from the Australian Research Council (ARC), which also funds T.K.'s research position. S.S. is funded under the ARC DP grant 0988751 which partially supports the HERMES project. M.R.K. is supported by the National Science Foundation through grant AST-0807739; by NASA through the Spitzer Space Telescope Theoretical Research Program, provided by a contract issued by the Jet Propulsion Laboratory; and by the Alfred P. Sloan Foundation through a Sloan Research Fellowship. J.B.H. is indebted to Merton College, Oxford for a Visiting Research Fellowship, and to the Leverhulme Trust for a Visiting Professorship to Oxford. J.B.H. and T.K. are grateful to the BIPAC Institute, Oxford, for their hospitality.

In the EnLink clustering scheme, each group is characterized by a maximum density ρ_max and a minimum density ρ_min, and these are used to define its significance S as

$$\begin{equation} S=\frac{\ln (\rho _{\rm max})-\ln (\rho _{\rm min})}{\sigma _{\ln (\rho)}}, \end{equation} \tag{ A1 }$$

where σ_ln(ρ) is the standard deviation associated with the density estimator and is a constant for a given dimensionality and q_den.

For Poisson-sampled data the distribution of density as estimated by the code using the kernel scheme is log-normal and the variance satisfies the relation

$$\begin{equation} \sigma _{\ln (\rho)} = \sqrt{{\cal V}_d ||W||_2^2/q_{\rm den}}, \end{equation} \tag{ A2 }$$

where q_den is the number of neighbours employed for density estimation, ${\cal V}_d$ the volume of a d-dimensional unit hypersphere and ||W||²₂ the L₂ norm of the kernel function (Sharma & Johnston 2009). For q_den = 6 and d = 2, σ_ln(ρ) evaluates to 0.4714.

Thus EnLink has two free parameters: (1) the significance threshold of the group S_T, and (2) the number of nearest neighbours q_den used for density estimation. The variable q_den should be set to less than the minimum desired size of the groups because the density is smoothed over a scale of q_den nearest neighbours such that the significance and hence the probability of detecting groups of size less than q_den is drastically reduced. Additionally, q_link = min(10, q_den − 1) neighbours are used for linking the groups, which means groups whose density peaks lie within q_link nearest neighbours of each other cannot be separated. Since we are interested in identifying groups with less than 10 data points, we set q_den = 6; the case for n_* ≲ 6 is discussed in the next section.

Our initial analysis is for data points with [Fe/H] < −2.5. This includes a crowded, high-density region of the abundance plane which skews the clustering analysis considerably. However, this may be an important regime if the first chemical elements arise from pair instability SNe (Karlsson et al. 2008). The significance of spurious groups is distributed almost as a Gaussian distribution and hence the significance threshold S_T is used to suppress these groups. An improved version of the empirical formula (Sharma & Johnston 2009) for the expected number of spurious groups, valid for S_T>1, is given by

$$\begin{equation} G(S>S_T)=\left(1-{\rm erf}\left(\frac{S_T f_{dq}}{\sqrt{2}}\right)\right)\frac{0.4n_*}{q_{\rm den}} \end{equation} \tag{ A3 }$$

where $f_{dq}=0.5 \sqrt{d(1-2.3/q_{\rm den})}$ is a small correction term and n_* is the number of data points.⁹ In Figures 11 and 12, we plot the significance distribution S of identified groups for data sets with different values of n_* and γ as labelled on the plots.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We set S_T such that the expected number of spurious groups (false positives) with S>S_T is about 2 in all data sets which we calibrate from the beta models. S_T is set to 1.0, 2.35, 3.5 and 4.35 for data with n_* as 100, 300, 1000 and 3000 respectively. This number is increasing because, for a fixed q_den and data dimensionality, the number of spurious groups due to Poisson noise increases linearly with the total number of data points n_*.

Each panel shows the outcome of five EnLink analyses: the cluster-less beta model (solid black line); the simulated models A1 and B1 having unrealistic abundance errors 0.01 dex shown as solid red and blue curves respectively; the simulated models A5 and B5 having abundance errors 0.05 dex shown as dotted red and blue lines respectively. In each panel, the mean number of identified groups G along with its dispersion σ_G is also given. The mean and dispersion were calculated using 100 random realizations of each data type. The plotted distributions are also averaged over these 100 random realizations. We also provide a measure of the statistical significance of detecting clusters in a data set as follows. If G_B is the number of groups as predicted by the smooth beta model, the statistical significance of clustering in a model is given by

$$\begin{equation} C_S=\frac{\langle G\rangle - \langle G_B\rangle }{\sqrt{\sigma _{G_B}^2+\sigma _{G}^2}} \end{equation} \tag{ A4 }$$

The third column shows the value of C_S for each data set.

It can be seen in Figure 11 and 12 that the expected number of spurious groups for the beta models are nearly independent of the value of γ, which is a consequence of the adaptive metric scheme. Specifically, EnLink uses the concept of a locally adaptive metric, which is used for calculating densities and nearest neighbours of data points (a refinement over earlier developments by Ascasibar & Binney 2005) to increase the efficiency of detecting clusters in multi-dimensional spaces. If instead clustering is performed using a Euclidean metric on raw data, the significance distribution of spurious groups due to Poisson noise is found to vary with the distribution of points in the abundance space, e.g., the beta models with different values of γ.

In the limit of few data points, the locally adaptive metric scheme is equivalent to using a metric which is given by the inverse of the dispersion along each dimension. As a check we also performed the analysis using this simpler scheme and found equivalent results, thereby demonstrating the robustness of our group-finding analysis. For the sake of accuracy, we use the beta models to evaluate the significance of clustering, but strictly speaking this is not required and any Poisson sampled cluster-less data would also suffice. In fact, the empirically derived formula given by Equation (A3) can also be used directly to predict the number of spurious groups.

Next we compare the parameter C_S for the different cases. It can be seen that as n_* and γ increase, C_S increases also. First we look at cases with unrealistic measurement errors of 0.01 dex. It is clear that the γ = 2.5 models are nearly undetectable for all values of n_* and closely resemble the beta models. For γ = 2.0, a minimum value of n_* = 300 is needed for C_S ≳ 1. For lower values of γ, signatures of clustering are visible with as few as 100 points. Increasing the abundance errors to 0.05 dex significantly affects our ability to detect clusters. When sufficient number of data points are present, the A5 and B5 models perform better for γ = 1.5 and 2.0 as compared to γ = 1. This is because the clusters in these cases are more in number, are spread over a larger area, are less crowded, and hence easier to detect.

Since clustering is most prominent for lower values of [Fe/H], we also investigated data sets with a lower metallicity cutoff, [Fe/H] < −3.0. These results are shown in Figures 13 and 14. The parameter C_S is found to increase in general for all cases, unlike what is seen in the [Fe/H] < −2.5 analysis (Figures 11 and 12). It is striking how much better EnLink performs on these simulations which is to be expected given the broader intrinsic dispersions of [r/Fe] in the models. In light of these results, another way to look at Figure 4 is the significance S of clustering over the abundance plane. In Figure 15, we show the S contours for the four different values of γ with lower significance regions to high [Fe/H]. Clustered abundance data points near the mean value of [r/Fe] are more likely at low [Fe/H] for low values of γ. Clustered abundance data points away from the mean are favoured by higher values of γ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size