THE s-PROCESS ENRICHMENT OF THE GLOBULAR CLUSTERS M4 AND M22

The assumption that globular clusters (GCs) are simple stellar populations (i.e., populations of stars that formed simultaneously from gas of a uniform chemical composition) has made them ideal laboratories for the study of low-mass stellar evolution (Moehler 2001) and enabled their ages to be accurately determined. This has aided cosmology by setting a lower limit on the age of the universe (Chaboyer et al. 1996; Dotter et al. 2010). However, the simple stellar population model of GCs has been undermined by spectroscopic studies that reveal significant star-to-star abundance variations (≳ 1 dex) in the light elements from C to Al (e.g., Cottrell & Da Costa 1981; Carretta et al. 2009b; Denissenkov & Hartwick 2014). Similar variations found in unevolved stars show that the chemical variations were initially present in the star-forming gas rather than being the result of nucleosynthesis and mixing within the observed stars (Cannon et al. 1998; Gratton et al. 2001). More recently, photometric studies have independently confirmed the existence of multiple populations in the form of split main sequences and subgiant branches in color–magnitude diagrams (e.g., Piotto et al. 2007; Piotto 2009; Milone et al. 2008).

The light element patterns that exist almost exclusively in GCs (i.e., rarely in field stars and open clusters, see Gratton et al. 2000; De Silva et al. 2009) include anticorrelations between the abundances of C and N, Na and O, and sometimes Mg and Al, typically with a C+N+O abundance that is constant within observational errors. The abundance patterns depict a H-burning process at high temperature (>80 MK⁴) combined with dilution by varying amounts of unprocessed material, although the stellar sites where this burning takes place and the mechanism of dilution are presently not well understood (Denisenkov & Denisenkova 1990; Langer et al. 1993; Decressin et al. 2007; Prantzos et al. 2007; D'Orazi & Marino 2010; D'Ercole et al. 2011).

In contrast to the light elements which vary in abundance, GCs are typically homogenous in [Fe/H]⁵ (σ < 0.05 dex; Carretta et al. 2009a) and in the abundances of neutron-capture elements (Z > 30; Gratton et al. 2004; Yong et al. 2006, 2008a; D'Orazi et al. 2010). Exceptions are known, including ω Centauri (Norris & Da Costa 1995; Smith et al. 2000; Johnson & Pilachowski 2010), M22 (Marino et al. 2009), NGC 1851 (Yong & Grundahl 2008; Villanova et al. 2010; Carretta et al. 2011), M2 (Yong et al. 2014), M15 (Sneden et al. 1997; Sobeck et al. 2011), and possibly NGC 2419 (Cohen & Kirby 2012).

Neutron-capture elements refer to elements with atomic number Z >30, because production of these elements is almost entirely by a process of neutron captures and β⁻-decay reactions. Depending on whether the average neutron-capture rates are less than or greater than the average rate of β⁻-decay reactions, the processes are divided into slow (s-process) and rapid (r-process) neutron-capture processes (Burbidge et al. 1957). Although most heavy elements can be synthesized by both processes, elements whose production is dominated by the r- or the s-process in solar system material are commonly referred to as r- and s-process elements.

Due to the large uncertainties involved in numerically modeling nucleosynthesis by the r-process, the r-only component of a heavy element distribution is often inferred from solar system material by subtracting the s-process component, which itself may be determined either theoretically (e.g., Arlandini et al. 1999; Goriely 1999; Sneden et al. 2008) or empirically (e.g., Simmerer et al. 2004).

The s-process takes place at low neutron densities (⩽10¹⁴ cm⁻³; Busso et al. 1999) and operates exclusively on nuclides that are very close to stability, as nuclei that become unstable following neutron capture have time to β-decay back to stability before additional neutron captures occur. In the build-up of progressively heavier elements via the s-process, bottlenecks form around nuclides with "magic" numbers of neutrons (e.g., 50, 82, 126) which form nuclear structures that are more stable against neutron capture than their neighbors (Busso et al. 1999). Three major peaks develop: a light s-peak (Sr, Y, Zr), a heavy s-peak (Ba, La, Ce), and a peak at Pb, with the light peak forming first and the heavier peaks forming later with increasing neutron exposure.

GCs provide laboratories to test and explore our understanding of stellar nucleosynthesis. One cluster that has been studied extensively is M22 (NGC 6656), which exhibits internal variation in [Fe/H] and s-process abundances that are bimodally distributed and neatly separate into two groups (Marino et al. 2009, 2011a; Da Costa & Marino 2011). While there are other well-studied clusters with Fe and s-process variation (e.g., ω Centauri), the simpler chemical evolution history of M22 relative to more complex systems like ω Centauri makes it an attractive system for testing theories about s-process variation in GCs more generally.

Even among GCs that are homogenous in their abundances of Fe and neutron-capture elements there exist puzzles surrounding the chemical evolution of the s-process elements. For example, M4 is a fairly typical mono-metallic metal-poor GC ([Fe/H] =−1.18; Carretta et al. 2009a), except that it has super-solar abundances of s-process peak elements (e.g., Rb, Y, Zr, La, Ba, Pb; Brown & Wallerstein 1992; Ivans et al. 1999). The origin of the s-process elements in M4 and M22 is speculated on in the literature, but often on the basis of individual stellar yields (e.g., Roederer et al. 2011) rather than a full investigation using a chemical evolution model. Very recently, Straniero et al. (2014) presented the first comparison of the s-process distributions of M4 and M22 with the summed contribution from a generation of asymptotic giant branch (AGB) stars at the metallicity of M22 ([Fe/H] =−1.8).

In this paper, we present simplified chemical evolution models of the heavy elements in GCs and predict the abundance variations that arise from s-process production by (1) massive stars with rotation, or (2) a generation of AGB stars that span a range of stellar masses. We then compare our chemical abundance predictions with the observed abundances of stars in M4 and M22. The success or failure of the individual models gives us insight into the stellar sites and timescales of s-process enrichment in GCs, as well as highlighting the shortcomings of current stellar nucleosynthesis models.

We define as massive stars those with sufficient mass to eventually form a collapsing core of Fe and end their lives as core-collapse supernovae. Current estimates for the lower-limit of initial mass required to meet this condition are around 8–12 M_☉, with lower masses required at lower metallicities (Langer 2012; Nomoto et al. 2013; Jones et al. 2013).

In massive stars, neutron-capture nucleosynthesis takes place during presupernova evolution and possibly also during the supernova itself. During convective core He-burning and shell He- and C-burning, neutrons are released via the ²²Ne(α, n)²⁵Mg reaction (Peters 1968; Raiteri et al. 1992; Meyer 1994; The et al. 2007).

The production of ²²Ne occurs via He-burning of ¹⁴N left over from H-burning in the CNO cycle. In models without rotation ²²Ne is secondary since its yield depends on the initial amount present plus any formed from α-capture onto ¹⁴N, which itself is limited by the initial abundance of C+N+O. Hence, there is very little s-process production at low metallicity in nonrotating models. Some production of heavy elements in massive stars does take place (the weak s-process) but this is mainly concentrated around elements of the first s-peak near Y, with virtually no heavy s-elements or Pb being produced (Beer et al. 1992; Pignatari et al. 2010).

In models of massive stars that do include rotation, rotationally induced mixing transports primary ¹²C and ¹⁶O produced in the convective He-core to the H-burning shell, where it is then converted into ¹⁴N via the CN-cycle (Meynet et al. 2006). The primary ¹⁴N is then mixed into and burned in the He core, resulting in an almost-primary production of ²²Ne (Hirschi 2007) that dramatically increases s-process yields at low metallicity.

Pignatari et al. (2008) present the first s-process yields for a rotating massive star with their 25 M_☉ model. They find that rotation increases s-process yields by orders of magnitude and alters the standard weak s-process distribution with a peak of production between Sr and Ba. The high production of heavy s-elements in their model is due to the use of the very low Descouvemont (1993) rate for the ¹⁷O(α, γ)²¹Ne reaction, which is disfavored by recent experiments. Frischknecht et al. (2012) present s-process yields from a set of massive models with updated reaction rates and find that rotation leads to the complete consumption of Fe seeds at metallicities below Z = 10⁻³ and an increase to the production of elements near the Ba peak at the expense of the Sr peak as metallicity decreases.

Figure 1 presents the heavy-element yields of 25 M_☉ massive star models with several initial rotation rates from Frischknecht et al. (2012). This figure demonstrates that under the condition of fast rotation, the s-process production in massive stars at low metallicity begins to include elements that would otherwise be associated uniquely with AGB stars (e.g., Ba, La, and Pb). For this reason, massive rotating stars must be considered as a possible source of the neutron-capture elements in GCs.

Zoom In Zoom Out Reset image size

In this study we use the presupernova yields of neutron-capture elements calculated from a grid of rotating massive stars (including those used to generate Figure 1) with initial masses of 15, 20, 25, and 40 M_☉ at initial metallicities of Z = 10⁻⁵ ([Fe/H] =−3.8) and Z = 10⁻³ ([Fe/H] =−1.8) with α-enhanced initial compositions as described in Frischknecht (2012). The rotation rates of the models are specified by their initial velocity at the equator as a fraction of the break-up velocity (v_crit, the velocity at which centrifugal force balances gravity).

For elements Z ⩽ 26, we use the zero-metallicity explosive yields of Limongi & Chieffi (2012). Although supernova yields presently carry large uncertainties, the effect of varying the Fe yields will be to scale our resulting heavy element distributions up and down while leaving the ratios between elements unchanged. The supernova shockwave will not significantly affect the s-process production and hence the s-process yields are approximated by their presupernova values although the mass yields, and to a lesser extent the s-process distribution, will depend on the mass cut (Tur et al. 2009).

In low- to intermediate-mass (0.8–8 M_☉) stars, the s-process takes place during the thermally pulsing AGB phase of evolution. For further details on AGB stellar evolution and nucleosynthesis, we refer to the reviews by Herwig (2005) and Karakas & Lattanzio (2014).

Figure 2 shows the average composition of the stellar ejecta of AGB models selected from the full grid which includes masses of 2.5, 2.75, 3.00, 3.25, 3.5, 4.0, 4.5, 5.0, 6.0, and 7.0 M_☉ at a metallicity of Z = 0.001 ([Fe/H] =−1.2 scaled solar), and are taken from Fishlock et al. (2014a). This figure displays the transition between the s-process yields from low-mass stars (≲4 M_☉) to intermediate-mass stars as a result of the ²²Ne neutron source becoming active. This transition mass also roughly coincides with our assumed upper limit initial masses for ¹³C pockets in AGB stars of 3 or 3.5 M_☉ at [Fe/H] =−1.2. For the rest of this section, we briefly summarize the operation of the s-process in AGB stars.

Zoom In Zoom Out Reset image size

With increasing initial mass, the maximum temperature in the intershell obtained during a thermal pulse also increases. A consequence is that in stars ≲4 M_☉, fewer thermal pulses are accompanied by a substantial activation of the ²²Ne(α, n)²⁵Mg reaction. Instead, free neutrons for the s-process are mainly released by radiative ¹³C-burning via the ¹³C(α, n)¹⁶O reaction, which is active at temperatures as low as 90 MK (Cameron 1955; Straniero et al. 1995). Producing the required ¹³C has been a challenge for stellar modelers, as the ¹³C abundance left behind by the H-burning shell is too low to allow for sufficient s-processing, and the convective region following a thermal pulse cannot extend into the H-rich region above the ¹²C-rich intershell (Iben 1975b).

Current AGB models achieve s-process nucleosynthesis via a ¹³C pocket in the following manner: protons from the envelope are "partially mixed" beyond the formal convective border into the ¹²C-rich intershell region, thus enabling the CN cycle reactions ¹²C(p, γ)¹³N(β⁺)¹³C (Gallino et al. 1998; Arlandini et al. 1999). The mixing process is required to have only marginal efficiency, otherwise the newly synthesized ¹³C is readily destroyed by further proton captures to make ¹⁴N, which is a neutron poison, i.e., its large neutron-capture cross section makes it an efficient absorber of free neutrons. For models that include a partial-mixing zone (PMZ) by inserting an exponential profile of protons below the inner edge of the envelope convective zone, a ¹³C pocket is formed below a pocket of ¹⁴N (Cristallo et al. 2009; Lugaro et al. 2012). In our post-process AGB nucleosynthesis models (including those used to calculate the yields shown in Figure 2) we include a PMZ of 2 × 10⁻³ M_☉ at the deepest extent of each third dredge-up episode for initial masses ⩽3.0 M_☉ and a 1 × 10⁻³ M_☉ PMZ for selected models at 3.25 and 3.5 M_☉. The technique we use to include a PMZ is identical to Lugaro et al. (2012) and we refer the reader to that paper for more details. We discuss the uncertainties related to ¹³C pockets in Section 4.

Figure 2 illustrates that low-mass stars produce significant quantities of heavy s-peak elements and Pb at low metallicity and confirms previous results in the literature (Travaglio et al. 2001; Van Eck et al. 2001, 2003; Lugaro et al. 2012). This is because the neutron source ¹³C is primary (independent of metallicity) while at low metallicity fewer Fe-seed nuclei (the most abundant heavy element) are available (Clayton 1988). With a large neutron supply per Fe seed, neutrons are preferentially captured by heavier nuclei and the abundance distribution is shifted toward higher atomic numbers.

Figure 2 also shows that the yields of elements heavier than Sr are significantly lower in models with masses ⩾4 M_☉ at a metallicity of [Fe/H] =−1.2. This is because the dominant neutron source in these models is the ²²Ne(α, n)²⁵Mg reaction, which is active at temperatures above about 300 MK (Cameron 1960; Iben 1975a; Goriely & Mowlavi 2000). The high temperatures and ²²Ne nuclei required to activate this source are found near the base of the He intershell in convective zones driven by He-shell flashes. Thus, neutrons are briefly exposed to a relatively large number of Fe seeds at the base of the flash-driven convective zone, and the resulting s-process distribution in intermediate-mass stars is mostly weighted toward the light s-peak near Sr–Y–Zr, with lower yields of heavy s-elements compared to lower mass stars.

In summary, the change from the ²²Ne source operating in convective pulses to radiative ¹³C-burning during the interpulse phase creates a dramatic change in the distribution of heavy elements between models above and below the transition mass of around 3–4 M_☉. The precise mass of this transition is dependent on the choice of the highest mass to include a ¹³C pocket, which is an uncertain parameter that is model and metallicity dependent (Goriely & Siess 2004; Herwig 2004).

The uncertainties that have the greatest effect on the yields of heavy elements are the numerical treatments of convection, mass loss, and reaction rates, as well as the rotation in massive stars and low-temperature opacities in AGB stars (e.g., Marigo 2002; Fishlock et al. 2014b; Constantino et al. 2014).

In stellar models of all masses, convective mixing plays a crucial role in the transport of energy and chemical species. The construction of accurate stellar models requires a method to approximate the effects of convection in 1D stellar evolution codes, as the high computational demands of full 3D hydrodynamical models limit their simulation times to no more than a small fraction of a stellar lifetime (e.g., Stancliffe et al. 2011). The most common numerical treatment of convection is the mixing-length theory (MLT) that depends on the value of an uncertain parameter, α, which is the mixing length in units of the local pressure scale height. The value of α is usually assumed to be constant on the AGB (e.g., the yields shown in Figure 2 use a value of 1.86), even though empirical and theoretical studies both suggest that the value changes with stellar evolution (Lebzelter & Wood 2007; Magic et al. 2014). Larger values of α have been shown to increase the depth of the third dredge-up (Boothroyd & Sackmann 1988), which increases the yields of s-process elements (Cristallo et al. 2009, 2011).

An alternative treatment of convective mixing that has been applied to AGB stars is the full-spectrum of turbulence (FST; Canuto & Mazzitelli 1991; Canuto et al. 1996). FST predicts a higher rate of energy transport than MLT, which leads to increased surface luminosities and higher interior temperatures in stellar models. In the intermediate-mass (4–6 M_☉) models of Ventura & D'Antona (2008) that use FST, temperatures at the base of the convective envelope reach 90–110 MK, which is hot enough for extensive H-burning nucleosynthesis (hot bottom burning). Combined with a luminosity-dependent mass-loss law, the high luminosities of these models drive rapid mass-loss rates that shorten the thermally pulsing AGB phase and reduce the number of dredge-up episodes (Ventura & D'Antona 2005; Ventura et al. 2013). The limited dredge-up in these models leads to a negligible net yield of C+N+O in the stellar wind. Presumably, this would also result in negligible yields of s-process elements, although yields from an FST model with a full s-process network are, to our knowledge, not published at present.

Another major uncertainty in stellar modeling is the mass-loss rate and its dependence upon stellar parameters. Indeed, massive stars can lose more than half of their mass by the end of core He burning (Chiosi & Maeder 1986). In massive stars with rotation, mass loss transports angular momentum away from the stellar surface (Hirschi 2007, to which we refer for details of the mass-loss prescription used in our massive star models). With the very low mass-loss rates expected at metallicities of Z = 10⁻⁵ ([Fe/H] =−3.8) and below, extremely metal-poor massive stars will evolve differently from observable OB stars (Maeder & Meynet 2000). Adding further complexity, the mass-loss rate would also be increased by the presence of a binary companion. For AGB stars, mass loss is very difficult to determine empirically without an accurate understanding of the dust composition and detailed models of the radiative transfer physics. Because the rate of mass loss controls the amount of time spent on the AGB and the number of thermal pulses, changes to the mass-loss rate have a significant effect on the predictions of stellar yields. In our models, we use Vassiliadis & Wood (1993) mass-loss rates along the AGB, which includes the switch to a superwind phase of extremely rapid mass loss near the tip of the AGB. An alternative is the Bloecker (1995) formula derived from dynamical calculations of the atmospheres of Mira-like stars, which predicts higher mass-loss rates and shorter AGB lifetimes.

Our massive star models include rotationally induced mixing in the form of meridional circulation and shear instabilities which dramatically alter the yields of CNO and s-process elements, depending on the rate of rotation (Frischknecht et al. 2012). The best constraints on the rotation rates of low-metallicity massive stars come from the comparison of chemical signatures in low-metallicity, low-mass stars with the predictions of rotating stellar models. In order to explain the existence of high N/O and C/O ratios at times too early for AGB stars to contribute, Chiappini et al. (2006) infer rotation rates of around 0.5 times the break-up velocity (v_rot/v_crit = 0.5) at [Fe/H] < − 3. Chiappini et al. (2008) claim that rotation is independently supported by the low ¹²C/¹³C ratios of metal-poor stars, which they report are consistent with models having rotational velocities of v_rot/v_crit ≃ 0.5 to 0.6. Fabbian et al. (2009) reach a less-certain conclusion about rotation and interpret high C/O ratios as possible signatures of either Population III stars or rotating Population II stars. The effect of rotation on the s-process yields is illustrated in Figure 1 which shows that rotation is the dominant effect.

For rotating massive star models, another uncertainty with an effect on s-process yield predictions is the competition between the ¹⁷O(α, γ)²¹Ne and ¹⁷O(α, n)²⁰Ne reactions. This is because ¹⁶O is highly effective at capturing free neutrons, which produces ¹⁷O. Neutrons are then either recycled via ¹⁷O(α, n)²⁰Ne or lost via ¹⁷O(α, γ)²¹Ne. The rate of the ¹⁷O(α, γ)²¹Ne reaction is particularly uncertain at the relatively low energies of stellar interiors. The first experimental rates for this reaction were published by Caughlan & Fowler (1988, hereafter CF88) and subsequently disputed by Descouvemont (1993), who predicted on theoretical grounds that the rate should be lowered by roughly a factor of 1000. However, more recent experimental work by Best et al. (2011) supports a rate similar to CF88. Best et al. (2013) report that the ratio between the (α, γ) and (α, n) reactions is best matched by using the CF88 rate divided by ten for ¹⁷O(α, γ)²¹Ne and the Angulo et al. (1999, NACRE) rate for ¹⁷O(α, n)²⁰Ne, the combination of which we will refer to as CF88/10 rates.

For s-process yields of both intermediate-mass AGB models and massive-star models, the ²²Ne(α, n)²⁵Mg reaction plays a critical role in determining neutron fluxes, and for this reason it has been the subject of a number of studies (Angulo et al. 1999; Jaeger et al. 2001; Koehler 2002; Karakas et al. 2006). Recent rates presented by Longland et al. (2012) have reduced the uncertainties in AGB model abundances caused by uncertainty in these reactions to less than a factor of two.

A major uncertainty for the s-process in low-mass models concerns the formation of a ¹³C pocket. This is because the physical mechanism that leads to ¹³C pockets in stars is yet to be identified. Currently proposed candidates include convective-boundary mixing (Herwig 2000; Cristallo et al. 2004), rotational mixing (Herwig & Langer 2001; Piersanti et al. 2013), or gravity-wave driven mixing (Denissenkov & Tout 2003). Eventually, a deeper understanding of the physics involved might completely eliminate the free parameter that determines the mass of the ¹³C pocket. At present, a variety of constraints have been derived from observations of carbon-enhanced metal-poor stars (Izzard et al. 2009; Bisterzo et al. 2012; Lugaro et al. 2012), planetary nebulae (Shingles & Karakas 2013; Miszalski et al. 2013), and post-AGB stars (Bonačić Marinović et al. 2007; De Smedt et al. 2012).

Aside from the uncertain size of the partial mixing zone and resulting ¹³C pocket, an additional uncertainty relates to the stellar initial masses in which a ¹³C pocket can be formed. With increasing stellar mass, the size of the He-rich intershell region decreases and temperatures at the base of the convective envelope during the third dredge-up increase, inhibiting ¹³C-pocket formation in more massive AGB stars. Goriely & Siess (2004) show that when the third dredge-up takes place with temperatures of around 40–70 MK ¹³C-pocket formation can be suppressed, depending on the details of any diffusive mixing near the convective boundary. At our metallicity of Z = 0.001, the results of Goriely & Siess (2004) suggest that ¹³C-pocket formation could become inhibited above around 3.0–3.5 M_☉ (but see Straniero et al. (2014) for a different view on ¹³C-pocket formation above this mass). To account for this uncertainty on our results, we separately consider two cases in which our nucleosynthesis post process includes a PMZ for all stellar masses up to 3.0 or 3.5 M_☉. This is an approximation in the absence of a physically accurate PMZ included in our stellar model calculations.

In this work we do not consider binary stars, although the presence of a binary companion will also alter the yields with a dependence on the period and mass ratio of the system.

5.1. Differential Abundances and Empirical s-Process Distributions

As an indication of how elemental abundances vary between two stars or stellar populations, it is common to subtract solar bracket [X/Fe] abundances (e.g., Yong et al. 2008a; Roederer et al. 2011). The difference [X/Fe]₂ − [X/Fe]₁ is equal to log₁₀[(X/Fe)₂/(X/Fe)₁], i.e., it measures of the number ratio of X to Fe in system 2 as a factor of the ratio in system 1. In the case that system 1 represents an initial composition that has undergone nucleosynthesis to make the abundances in system 2, a quantity that isolates the net production or destruction of elements is obtained by subtracting the number ratios in linear abundance space, i.e., Δ(X/Fe) = (X/Fe)₂ − (X/Fe)₁, assuming that Fe is either constant or only marginally produced. This quantity is analogous to the net yields of stellar nucleosynthesis models, which are computed by subtracting the abundances in the initial composition from the abundances in the stellar ejecta (e.g., Karakas 2010).

Using a linear abundance subtraction, Roederer et al. (2011, Table 8) calculate an s-process-only residual composition for M22 by subtracting the average X/H number ratios of s-poor from s-rich stars. We use the same technique to derive empirical s-process distributions for M4 and M22, except that we use number ratios relative to Fe. Our own testing confirms that the resulting distributions look very similar regardless of whether abundances relative to Fe or H are used.

Figure 3 shows our calculated s-process-only residuals of M22 〈s-rich〉 − 〈s-poor〉 and 〈M4〉 − 〈M5〉 relative to the solar abundances (Asplund et al. 2009). To visually emphasize small differences (that are significant within the errors), we plot on a linear scale. The distance from the zero point is related to the amount of dilution with s-poor material, while the shape of the distribution is relatively independent of this uncertain parameter and primarily depends on the relative abundances in the stellar ejecta. In agreement with Roederer et al. (2011), we interpret the empirical s-process distributions of M4 and M22 as representing enrichment by material of a similar but not identical composition. We suggest that the s-process distributions of these two clusters are distinguishable as representing the results of different nucleosynthetic sites or stellar mass ranges. We now discuss the observations, starting with M4.

Zoom In Zoom Out Reset image size

5.2. The s-Rich Globular Cluster M4

5.3. The Two Populations in M22

We present abundance evolution results calculated using a new code, Evel ChemEvol to solve the equations of chemical evolution for a single-zone (for an review, we refer to Pagel 2009). Our testing with the AGB yields and self-pollution scenario described by Fenner et al. (2004) confirms that the code correctly reproduces the abundance results of an existing chemical evolution code. For more details of the output validation tests, see the Appendix.

Our simplified chemical evolution model includes a single short burst of star formation as a first-order attempt at understanding the enrichment of GCs. The final abundance outputs of the chemical evolution model represent the IMF-weighted (Kroupa et al. 1993) sum of ejecta from a range of stellar masses with yields that are interpolated from a grid of stellar models.

Our derivation of an s-process-only component from the observational abundances and the similar subtraction of the initial abundances from the final abundances of the models (or the subtraction of the final abundances of two different models) enables us to compare our chemical evolution predictions with both cluster systems simultaneously, although the initial composition will affect the ratios of elements in the stellar yields (e.g., [ls/hs]).

For each stellar mass in the range from 15 to 40 M_☉, massive star yields are interpolated from our grid of stellar models with initial rotation rates of 0.0 and 0.4 as a fraction of v_crit (Frischknecht et al. 2012). For the particular initial mass of 25 M_☉, we have also have yields from stellar models with rotation rates of 0.4 and 0.5 v_crit with and without alternative reaction rates (CF88/10). From the yields of the 25 M_☉ stellar models, we calculate a set of factors (one per chemical species) that approximate the effect of these alternative parameters on the yields of the other models in the grid with different initial masses.

Table 1 and Figures 4–6 show the quantitative results of our chemical evolution models with rotating massive stars and AGB stars. The first two rows show the observation results in terms of [ls/hs]⁶ and [Pb/hs] ratios of the s-process residuals. A timescale is given for the models with AGB stars, which is the stellar lifetime of the lowest included mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Rotating Massive Stars

6.2. AGB Stars

We test chemical evolution models that predict the output of a single generation of low-metallicity AGB stars, with the results provided in Figure 6 and Table 1. We vary the lower limit of the stellar mass range as a free parameter because this corresponds to the uncertain age difference between the s-process polluters and the s-process-rich stars (minus the gas cooling time). Because of the uncertainty over the upper mass limit for AGB stars to have a ¹³C pocket, we separately test chemical evolution models in which the 3.25 and 3.5 M_☉ yields are calculated from models with and without a PMZ of 1 × 10⁻³ M_☉.

For M4, the [ls/hs] and [Pb/hs] ratios are bracketed from above and below by models with AGB yields that have lower limit masses of 3.00 and 3.25 M_☉, respectively. From the stellar lifetimes, this corresponds to a minimum enrichment timescale 239–290 Myr. As the 3.00 M_☉ stellar model includes a PMZ and the 3.25 M_☉ model does not, this indicates a small contribution from stars with a ¹³C pocket. If the models up to 3.5 M_☉ include a PMZ, the [ls/hs] and [Pb/hs] ratios of M4 are bracketed by 3.5 and 4.0 M_☉ lower-limit models, corresponding to a 140–200 Myr minimum enrichment timescale. With the uncertain upper mass limit for the ¹³C pocket formation, the minimum enrichment timescale for M4 is likely around 140–290 Myr.

Although our AGB yields are not an exact match to the metallicity of M22 ([Fe/H] =−1.2 versus −1.8 in M22's s-poor group), we explore the similarities between our chemical evolution results and the observational data. The chemical evolution of heavy elements in the s-rich group of M22 requires even lower mass stars than M4, however a simultaneous match of [ls/hs] and [Pb/hs] ratios is not found in our results. The [ls/hs] ratio in M22 is bracketed by 2.75 and 3.0 M_☉ lower-limit models, while the [Pb/hs] ratio is bracketed by models with lower mass limits of 3.00 and 3.25 M_☉. If the stellar masses up to 3.5 M_☉ include a PMZ, we find that M22's [ls/hs] ratio is between those of the 3.00 and 3.25 M_☉ lower-limit models, while the predicted [Pb/hs] of these models is still too high to match the data. In both of our test cases for the upper limit mass of ¹³C pocket formation, the dual contribution from stars with a ¹³C pocket as well as stars with a ²²Ne neutron source are required. With our assumption that ¹³C pockets transition from fully developed to negligible between initial masses of 3.0 and 3.5 M_☉, we predict a lower limit on the polluter masses of 2.75–3.25 M_☉, which corresponds to a minimum enrichment timescale of 240–360 Myr.

A common method for comparing measured abundances with the predictions of stellar models is to use the yield results of a single stellar model rather than a grid covering a range of stellar masses that has been weighted by an initial mass function. In Table 1 and Figure 7, we present single-mass yield results for comparison with our chemical evolution results. The slope of the IMF means that the lowest contributing mass will have the largest contribution to the final abundances, however the single models with a ¹³C pocket importantly lack the significant production of light s-elements that is due to intermediate-mass AGB stars. This difference is apparent in the high Rb/Sr ratio and overall higher abundances of Rb, Sr, Y, Zr, and Mo by the 3.0–7.0 M_☉ model shown in Figure 6, as compared with same ratio from the single 3.0 M_☉ model shown in Figure 7.

Zoom In Zoom Out Reset image size

We have used the shape of the s-process distributions in M4 and M22 to identify the s-process polluter mass range and the corresponding maximum stellar lifetime, which places a lower limit on the timescale of s-process enrichment. As well as s-process enrichment, M4 and the two groups in M22 also feature anti-correlated variations in O and Na, which we do not attempt to explain. Although we match the s-process distribution of M22's s-rich group with the ejecta of AGB stars, the coexistence of an Fe variation in M22 likely required some fraction of the ejecta from massive stars to be kept within the cluster to form new stars with a higher Fe abundance.

The lower-limit mass range of 2.75–3.25 M_☉ in our best-fitting models for M22 corresponds to a stellar lifetime in the range 300 ± 60 Myr. Assuming that the time for the ejecta to cool and form new stars is relatively small, this value is consistent with the 300 Myr upper limit derived from isochrone fitting of the subgiant branch region by Marino et al. (2012) and Joo & Lee (2013). A match between the intergroup age difference and the lifetime of the minimum polluter mass supports a scenario in which the metal-rich group in M22 has been self-enriched with material ejected from stars coeval with the present-day s-poor group. Alternative scenarios in which the two metallicity groups in M22 (which have independent light element anti-correlations) are the result of a merger of two separate GC systems or the second generation is formed from s-process rich material accreted from outside the cluster are also plausible. However, under both of these alternative scenarios, the close match between timescales of pollution and the age difference between the stellar groups would be a coincidence.

Further evidence for a lower mass limit of ≈3 M_☉ and an enrichment timescale of ≈300 Myr for M22 is the measured 0.6 dex spread in F abundances (Alves-Brito et al. 2012; D'Orazi et al. 2013b). D'Orazi et al. (2013b) report F abundances that correlate with O, are anticorrelated with Na, and increase between the two groups. The authors suggest that the s-rich group has been enriched by the ejecta of stars with masses between 4 and 5 M_☉, as these stars would destroy (rather than produce) F while O is destroyed in the early stages of GC formation. However, these measurements could be heavily affected by systematic errors as is claimed by de Laverny & Recio-Blanco (2013), who argue that a reliable detection of the HF line in M22 stars is unlikely due to errors in radial velocity correction, continuum subtraction, and the removal of telluric absorption lines.

The matching of M4 − M5 to a model of AGB ejecta opens up the question of how the formation of M4 differed to that of M5. A scenario similar to M22 in which an s-poor generation of stars pollutes the interstellar medium from which a second generation forms is ruled out by observations of constant s-process abundances in M4, which do not feature the same bimodality found in M22. A more likely scenario is that M4 and M5 formed out of material in an inhomogeneous early Galactic halo. James et al. (2004) show that Ba and Eu abundances plotted as a function of [Fe/H] for mono-metallic GCs (including M4 and M5) fall within the spread of halo field star values, suggesting that they share a common origin or a similar enrichment process.

Our inferred enrichment timescales for M4 and M22 are roughly a factor of two larger than the 150 ± 50 Myr for both clusters inferred by Straniero et al. (2014). They require a larger minimum contributing mass (4.0 ± 0.5 M_☉) due to their inclusion of a prescription for core–envelope convective boundary mixing (Cristallo et al. 2009) that predicts small ¹³C pockets in AGB models with masses as high as 4.5 M_☉. The predictions of Straniero et al. (2014) and those in this paper are both consistent within the uncertainty of the age spread in M22 derived from isochrone fitting. Our results support their conclusion that neutron captures from both ¹³C pockets and the ²²Ne source operating in convective pulses are required to explain the s-process enrichment of M4 and M22.

While our massive star models could not reproduce the s-process enhancements seen in M4 and in M22, there are still large uncertainties on the yields of s-process elements from rotating massive star models. For example, the yields of Pignatari et al. (2008) show a ratio between Y and Ba of approximately unity (see their Figure 2). This suggests that while AGB stars produce the best fit with our adopted stellar yields, other sets of yields may change our conclusions as to the nature of the polluters of heavy elements in GCs.

We consider the effect of a possible r-process difference between the s-rich and s-poor samples of up to [r/Fe] =0.4. Using the solar system r-process fractions of Simmerer et al. (2004), the effect would be to increase [Y/Fe] by 0.15 dex, [Zr/Fe] by 0.07 dex, [Ba/Fe] by 0.09 dex, [La/Fe] by 0.14 dex, and [Ce/Fe] by 0.11 dex. The net result for the [ls/hs] ratio would be a change of less than 0.01 dex. The ratio [Pb/Fe] would increase by 0.12 dex, resulting in a [Pb/hs] change of less than 0.01 dex. We conclude that our results hold independently of a possible r-process difference between M4 and M5 or the two populations in M22. A dilution by pristine material would shift the [X/Fe] ratios in the s-process residual to closer to zero, but would not affect the relative abundances between elements.

Our models predict [Pb/hs] ratios that are too high to match the observations of M22. A similar phenomenon is reported by De Smedt et al. (2014) for metal-poor ([Fe/H] <−1) post-AGB stars in the Magellanic clouds, which they refer to as the "lead discrepancy." If the Pb measurements are correct, then a solution to the lead discrepancy will likely require a better understanding of the mixing that leads to the formation of a ¹³C pocket, possibly by modeling it as an advective process, rather than more typical diffusive treatment. One form of extra mixing that is not included in our AGB models (or in most AGB stellar models) is the mixing due to rotation. The study of AGB models with rotation by Piersanti et al. (2013) hints at a possible solution to the lead discrepancy, as they find that rotation reduces the [Pb/hs] ratio in the stellar yields.

Rb is overproduced in our best-fitting model in comparison with the observational data for M4, while separate Rb abundances for the two groups in M22 are not available in the literature. An overproduction of Rb in AGB stellar models is also noted by D'Orazi et al. (2013a) under the assumption that AGB stars are responsible for the light element variations in M4. Their 6 M_☉ model with the mass-loss rate from Vassiliadis & Wood (1993) and a mixing-length parameter of α = 1.75 produces too much variation in neutron-capture elements, for too little variation in Na. Their solution is to use the higher mass-loss rates of Bloecker (1995) and a boosted mixing-length parameter (α = 2.2), which improves the fit to the abundances in M4 by increasing the temperature at the base of the convective envelope and reducing the cumulative dredge-up of s-process elements into the envelope.

The opposite case of Rb underproduction is found when stellar models (5–9 M_☉) are compared with AGB stars in the Galaxy and Magellanic clouds (van Raai et al. 2012; Karakas et al. 2012), although recent work by Zamora et al. (2014) suggests that the inferred Rb abundances may be systematically overestimated due to the presence of circumstellar envelopes. The implementation of a delayed superwind to increase Rb yields explored by Karakas et al. (2012) would likely worsen the discrepancy in our results, unless there was a simultaneous reduction in Rb production by the less massive (<5 M_☉) AGB models.

Future stellar models at the correct metallicity for M22, and more generally improvements to the numerical treatment of mixing and mass loss, might help to reduce some of the discrepancies with observations of the s-process abundances in globular clusters. Even with the current models, the application of similar techniques to other clusters with s-process variation such as M2 (Yong et al. 2014) and ω Centauri would enable us to characterize the full range of enrichment timescales and polluter masses among the anomalous GCs.

This research has made use of NASA's Astrophysics Data System. L.J.S. and A.I.K. thank Chris Sneden for helpful discussions about spectroscopic uncertainties. A.I.K. was supported through an Australian Research Council Future Fellowship (FT110100475). R.H. acknowledges the support from Eurocore project Eurogenesis and ERC Starting grant No. 306901. R.H. and A.I.K. acknowledge support from the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan.

To validate the output of the new chemical evolution code Evel ChemEvol used in this study, we use the stellar yields and GC self-pollution scenario described by Fenner et al. (2004). A metal-free initial composition is first polluted with the ejecta of massive stars up to a metallicity of [Fe/H] =−1.4. Subsequently, star formation takes place on a timescale of 10⁷ yr. The ejecta from stars <6 M_☉ is kept within the system, while ejecta from more massive stars is lost.

Our chemical evolution results in Figures 8–10 for Evel ChemEvol correspond almost exactly to Figures 1, 3, and 4 of Fenner et al. (2004). Small differences in the output can be explained by differences in the stellar lifetime function and the treatment of the massive star pollution phase, which were not specified in detail in the original paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results of this comparison give us confidence that the Evel ChemEvol is producing the correct abundance outputs and can be used to explore new chemical evolution scenarios.