Simultaneous Constraints on the Star Formation History and Nucleosynthesis of Sculptor dSph

In this work, we primarily use literature abundances derived from medium-resolution spectroscopy with the DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003) on the Keck II telescope. A number of previous works have obtained spectra of red giant branch stars in several globular clusters and classical dSphs. We compile several abundance ratios ⁶ from these catalogs: [Fe/H] and [αあるふぁ/Fe] ([Mg/Fe], [Si/Fe], and [Ca/Fe]) abundances from Kirby et al. (2010), [C/Fe] from Kirby et al. (2015), [Ni/Fe] from Kirby et al. (2018), [Mn/Fe] from de los Reyes et al. (2020), and [Ba/Fe] from Duggan et al. (2018). We also use supplemental data from the DART survey (Tolstoy et al. 2006), which used ESO VLT/FLAMES to obtain high-resolution (R ≳ 20,000) spectra of RGB stars in dSphs. The DART abundance ratios for Sculptor dSph are presented in North et al. (2012), who measured [Mn/Fe], and Hill et al. (2019), who measured all the other abundances used in this work.

In particular, we use the DART data to modify the [Ba/Fe] abundances for our analysis. Although the majority of barium is produced in the s-process in AGB stars (see, e.g., Table 10 of Simmerer et al. 2004), r-process nucleosynthesis also contributes to the production of barium—particularly at low metallicities—and the sites and yields of r-process nucleosynthesis are poorly constrained (Cowan et al. 2021). For this reason, we opt not to include the r-process in our GCE model. To accommodate this choice, we remove the r-process contributions to the barium yields by using measurements of europium (Eu), which is almost entirely produced by the r-process (Simmerer et al. 2004). We use the available Eu measurements from the DART survey (Hill et al. 2019) to compute [Ba/Eu], which is an indicator of the ratio of s-process to r-process contributions. We follow the procedure outlined in Duggan et al. (2018) to convert [Ba/Eu] to the fraction of barium produced from the r-process:

$$\begin{eqnarray}&&{f}_{r}=\displaystyle \frac{\tfrac{{N}_{s}(\mathrm{Eu})}{{N}_{s}(\mathrm{Ba})}-{10}^{{[\mathrm{Ba}/\mathrm{Eu}]}_{\odot }-[\mathrm{Ba}/\mathrm{Eu}]}}{\tfrac{{N}_{s}(\mathrm{Eu})}{{N}_{s}(\mathrm{Ba})}-\tfrac{{N}_{r}(\mathrm{Eu})}{{N}_{r}(\mathrm{Ba})}}\end{eqnarray} \tag{ 1 }$$

where N_s (X) and N_r (X) are the solar s-process and r-process number abundances of element X, obtained from Table 10 of Simmerer et al. (2004).

Since the number of available Eu measurements is relatively small, we compute a simple statistical correction by fitting a line to [Ba/Eu] as a function of metallicity. This is shown in Figure 3, where we derive the best-fit line ⁷ [Ba/Eu] = 0.52[Fe/H] + 0.42, which can be used to determine the [Ba/Eu] ratio for all stars with metallicity measurements. The s-process-only barium yield, [Ba/Fe]_s , can then be computed using the fraction f_r (Equation (1)):

$$\begin{eqnarray}&&{[\mathrm{Ba}/\mathrm{Fe}]}_{s}=[\mathrm{Ba}/\mathrm{Fe}]+\mathrm{log}(1-{f}_{r}).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Table 1 lists the full catalog of all abundances, including the original [Ba/Fe] abundances as well as s-process-only [Ba/Fe]_s. The estimated uncertainties include both statistical and systematic uncertainties, as reported in the original catalog papers. We note that we use [C/Fe] abundances that Kirby et al. (2015) have corrected for astration (the depletion of carbon by stars moving up the red giant branch; Suntzeff 1981; Carbon et al. 1982; Smith & Briley 2006) using the corrections proposed by Placco et al. (2014).

We now consider the simple GCE model used to fit the data described above. Conceptually this model is similar to that used by Kirby et al. (2011), and we refer the reader to that work for more details about the individual equations.

The model treats each dwarf galaxy as a chemically homogeneous, open-box system. In this system, the gas-phase abundance of each element is tracked over a discrete grid of time steps (Δでるたt = 1 Myr). Gas inflows and ejecta from CCSNe, Type Ia SNe, and AGB stars can contribute to the gas-phase abundance of each element, while gas outflows and star formation remove gas-phase elements. The full model is therefore described by

$$\begin{eqnarray}\begin{array}{rcl}{\xi }_{j}(t) & = & {\int }_{0}^{t}\left(-{\dot{\xi }}_{j,\mathrm{SF}}+{\dot{\xi }}_{j,\mathrm{II}}+{\dot{\xi }}_{j,\mathrm{Ia}}+{\dot{\xi }}_{j,\mathrm{AGB}}\right.\\ & & +\left.\ {\dot{\xi }}_{j,\mathrm{in}}-{\dot{\xi }}_{j,\mathrm{out}}\right)\end{array}\end{eqnarray} \tag{ 3 }$$

where ξくしー_j (t) is the gas-phase abundance of element j.

The other terms in Equation (3) describe the processes that contribute or remove gas-phase elements. Star formation is described by a Schmidt-like power law (Schmidt 1959; Kennicutt 1998):

$$\begin{eqnarray}&&{\dot{\xi }}_{j,{\rm{SF}}}(t)={A}_{\star }{\left(\displaystyle \frac{{M}_{j,{\rm{gas}}}(t)}{{10}^{6}\,{M}_{\odot }}\right)}^{\alpha }.\end{eqnarray} \tag{ 4 }$$

Following Kirby et al. (2011), the rate of pristine gas inflows is parameterized with a fast increase and slow decline (e.g., Lynden-Bell 1975):

$$\begin{eqnarray}&&{\dot{\xi }}_{j,\mathrm{in}}(t)={A}_{\mathrm{in}}\displaystyle \frac{{M}_{j,\mathrm{gas}}(0)}{{M}_{\mathrm{gas}}(0)}\left(\displaystyle \frac{t}{\mathrm{Gyr}}\right)\exp \left(\displaystyle \frac{-t}{{\tau }_{\mathrm{in}}}\right)\end{eqnarray} \tag{ 5 }$$

where ${M}_{j,{\rm{gas}}}(0)/{M}_{{\rm{gas}}}(0)$ is the initial mass fraction of element j, indicating that the inflows are primordial. ⁸ We assume that gas outflows are predominantly caused by supernovae, so the outflow is assumed to be linearly proportional to the supernova rate:

$$\begin{eqnarray}&&{\dot{\xi }}_{j,\mathrm{out}}(t)={A}_{\mathrm{out}}\displaystyle \frac{{M}_{j,\mathrm{gas}}(t)}{{M}_{\mathrm{gas}}(t)}\left({\dot{N}}_{\mathrm{II}}+{\dot{N}}_{\mathrm{Ia}}\right)\end{eqnarray} \tag{ 6 }$$

From Equations (4), (5), and (6), we define the variables {A_⋆, αあるふぁ, A_in, τたう_in, A_out} as free parameters in the model.

The ejecta from SNe and AGB stars at a given time step (${\dot{\xi }}_{j,\mathrm{II}}$, ${\dot{\xi }}_{j,\mathrm{Ia}}$, ${\dot{\xi }}_{j,\mathrm{AGB}}$) depend on the numbers of SNe/AGB stars occurring at that time, which are determined by convolving the past SFH with a DTD. For a given type of astrophysical event, the DTD describes the expected event rate as a function of τたう, where τたう is the delay time after a δでるた-function burst of star formation. The DTD for Type Ia SNe is observed to be a power law with index ∼–1 (e.g., Maoz et al. 2010):

$$\begin{eqnarray}\begin{array}{c}{{\rm{\Psi }}}_{{\rm{Ia}}}=\left\{\begin{array}{ll}0 & {t}_{{\rm{delay}}}\lt 0.1\,{\rm{Gyr}}\\ ({10}^{-3}\,{{\rm{Gyr}}}^{-1}\,{M}_{\odot }^{-1}){\left(\displaystyle \frac{{t}_{{\rm{delay}}}}{{\rm{Gyr}}}\right)}^{-1.1} & {t}_{{\rm{delay}}}\geqslant 0.1\,{\rm{Gyr}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 7 }$$

The exact parameterization of the Type Ia DTD is still an open question, and we discuss the effects of modifying the values in Equation (7) on our results later in Section 4.2.2.

For CCSNe and AGB stars, the DTD is primarily set by the stellar IMF, since stellar lifetimes depend strongly on stellar mass. For ease of comparison with Kirby et al. (2011), we use a Kroupa et al. (1993) IMF; we consider the effect of changing the IMF in Section 4.2. We further assume that all stars with birth masses between 10 and 100 M_⊙ explode as CCSNe at the end of their lifetimes, and that all stars with birth masses ≤10 M_⊙ eject mass through AGB winds within the final 1 Myr time step of their lifetimes (Marigo & Girardi 2007). Stellar lifetimes are then parameterized as a function of mass using the following equations (Padovani & Matteucci 1993):

$$\begin{eqnarray}&&{\tau }_{\star }(M)=(1.2{\left(M/{M}_{\odot }\right)}^{-1.85}+0.003)\,\mathrm{Gyr}\end{eqnarray} \tag{ 8 }$$

for stars with M > 6.6 M_⊙, and

$$\begin{eqnarray}&&{\tau }_{\star }(M)={10}^{\tfrac{0.334-\sqrt{1.790-0.2232(7.764-{\mathrm{log}}_{10}(M/{M}_{\odot }))}}{0.1116}}\,\mathrm{Gyr}\end{eqnarray} \tag{ 9 }$$

for stars with M ≤ 6.6 M_⊙. Figure 4 compares these equations (black dashed line) with numerical results from the Binary Population and Spectral Synthesis (BPASS) code (Eldridge et al. 2017; Stanway & Eldridge 2018), showing that Equations (8) and (9) are consistent with stellar evolution models.

Zoom In Zoom Out Reset image size

To determine the number of SNe/AGB stars at each time step, Kirby et al. (2011) computed the full convolution of the past SFH with the DTD (see their Equations (7), (10), and (13)). In this work, we instead track the elemental abundances in a forward-looking array; at each time step, we compute the number of stars that will produce SNe or AGB winds in the future, and we add the nucleosynthetic yields from these SNe/AGB stars (see next section) to the appropriate future times in the array. This approach, similar to that of the One-zone Model for the Evolution of GAlaxies code (OMEGA; Côté et al. 2017), eliminates most numerical integration from the model. Consequently, the computation of the GCE code is approximately an order of magnitude faster, making it possible to run it many times for a Markov Chain Monte Carlo sampler (see Section 3.1).

A key component of our GCE model is the set of nucleosynthetic yields from supernovae and AGB stars. A number of models have predicted yield sets, and we summarize a subset of these models in Table 2. The physical assumptions and computational limitations inherent in these models can produce significant uncertainties in their predicted yields (see Figures 12–14 in the Appendix).

Rather than selecting uncertain yield sets, we instead choose to parameterize the yields with representative analytic expressions, which we will allow to vary in our final fit. First, we fit analytic functions of mass and metallicity to the yields plotted in Figures 12, 13, and 14. For each nucleosynthetic source (CCSNe, AGB stars, and Type Ia SNe), we then vary the input yield sets in the GCE model to determine which elemental abundances are the most sensitive to variations in the yields. For example, for CCSNe, we run the GCE model twice using the same set of fiducial parameters {A_in, τたう_in, A_out, A_⋆, αあるふぁ} (for simplicity, we use the parameters measured by Kirby et al. 2011), varying only the input CCSN yields—either those from Nomoto et al. (2013) or those from Limongi & Chieffi (2018). We find that the abundances of C, Mg, and Ca predicted by the GCE model are sensitive to the input CCSN yields—that is, the predicted [X/Fe] abundances change by >0.2 dex at the peak of the metallicity distribution function (MDF) (−1.5 < [Fe/H] < −1.0). We therefore define parameters in our analytic functions for C, Mg, and Ca that can be varied to match all of the input yield sets.

The Appendix describes the final analytic functions and parameters in more detail. We obtain seven parameters that represent variations in the nucleosynthetic yields: the Fe yield from Type Ia SNe (Fe_Ia), the exponent of the C yields from CCSNe (expC_II), the normalization of Mg and Ca yields from CCSNe (normMg_II, normCa_II), the normalization of the C yield from AGB stars (normC_AGB), and the normalization and peak of the Ba yield from AGB stars (normBa_AGB, meanBa_AGB). These parameters are then allowed to vary in the GCE model.

We use the chemical evolution model described in the previous section to simultaneously match the MDF as well as the abundance trends of [Mg/Fe], [Si/Fe], [Ca/Fe], [C/Fe], and s-process-only [Ba/Fe]_s (see Section 2.1) as a function of [Fe/H]. We choose not to include [Mn/Fe] and [Ni/Fe] in the model fitting. Manganese and nickel are iron-peak elements that are produced in the same nucleosynthetic events as iron. However, unlike iron, these elements—particularly manganese—are likely more sensitive to the physics of Type Ia supernova than to the SFH (e.g., Seitenzahl et al. 2013; Seitenzahl & Townsley 2017). Rather than use [Mn/Fe] and [Ni/Fe] to fit the model, we instead use these abundances to validate our model and probe additional physics in Section 5.2.

Following Kirby et al. (2011), we treat the chemical evolution model as tracing a path _j (t) in the six-dimensional {[Fe/H], [Mg/Fe], [Si/Fe], [Ca/Fe], [C/Fe], [Ba/Fe]_s} space. The probability of a star forming at any time t is given by dP/dt = ${\dot{M}}_{\star }$(t)/M_⋆, where M_⋆ is the final stellar mass. The likelihood of a star i forming along the path defined by the chemical evolution model is therefore given by the line integral of dP/dt along the path _j for the total duration of the model t_final:

$$\begin{eqnarray}&&{L}_{i}={\int }_{0}^{{t}_{\mathrm{final}}}\left(\prod _{j}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{i,j}}\exp \displaystyle \frac{-{\left({\epsilon }_{i,j}-{\epsilon }_{j}(t)\right)}^{2}}{2{\left({\sigma }_{i,j}\right)}^{2}}\right)\displaystyle \frac{\dot{{M}_{\star }}(t)}{{M}_{\star }}{dt}\end{eqnarray} \tag{ 10 }$$

where _i,j is the jth observed elemental abundance ratio for star i, and σしぐま_i,j is the corresponding uncertainty. The final time step, t_final is not a free parameter. Rather, it is the last time step before the galaxy runs out of gas.

The total likelihood for a model is therefore proportional to the product of L_i for all N stars:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}L & = & \displaystyle \prod _{i}^{N}{L}_{i}\times \left(\displaystyle \frac{1}{\sqrt{2\pi }\delta {M}_{\star ,{\rm{obs}}}}\exp \displaystyle \frac{-{({M}_{\star ,{\rm{obs}}}-{M}_{\star ,{\rm{model}}})}^{2}}{2{(\delta {M}_{\star ,{\rm{obs}}})}^{2}}\right.\\ & & {\left.\times \displaystyle \frac{1}{\sqrt{2\pi }\delta {M}_{{\rm{gas}},{\rm{model}}}}\exp \displaystyle \frac{-({M}_{{\rm{gas}},{\rm{model}}})}{2(\delta {M}_{{\rm{gas}},{\rm{obs}}})}\right)}^{0.1N}.\end{array}\end{array}\end{eqnarray} \tag{ 11 }$$

Here, the additional terms require that the final stellar and gas masses of the model (M_⋆,model and M_gas,model) match the observed masses (M_⋆,obs and M_gas,obs) within the observational uncertainties. The exponent 0.1N is chosen to weight these terms relative to the abundance distributions, to prevent them from dominating the likelihood while ensuring that the models end up with approximately the correct stellar and gas masses. We use observed stellar masses and uncertainties from Woo et al. (2008). We assume the observed gas mass M_gas,obs is zero for all dSphs, and we choose an arbitrary uncertainty of δでるたM_gas,obs = 10³ M_⊙ to ensure the model converges. This is similar to observed upper limits on gas measurements (e.g., Putman et al. 2021, who find an upper limit of (3.2 ± 0.4) × 10³ M_⊙ on Sculptor's H i mass).

To estimate the values of the 12 free parameters {A_in, τたう_in, A_out, A_⋆, αあるふぁ, Fe_Ia, expC_II, normMg_II, normCa_II, normC_AGB, normBa_AGB, meanBa_AGB} that minimize the negative log-likelihood (−ln L), we used the emcee Python module (Foreman-Mackey et al. 2013) to implement a Markov Chain Monte Carlo (MCMC) ensemble sampler. Table 3 describes the inputs and outputs of this MCMC sampling: the priors, initial values from linear optimization, and the best-fit values for each parameter.

Table 3. MCMC Free Parameters

Parameter	Description	Prior	Initial Value	Best-fit Result
Ain	Normalization of gas infall rate (109 M⊙ Gyr−1)	${ \mathcal U }(0,5)$	1.07	${0.53}_{-0.08}^{+0.09}$
τたうin	Gas infall time constant (Gyr)	${ \mathcal U }(0,1)$	0.16	${0.27}_{-0.02}^{+0.03}$
Aout	Gas lost per supernova (${M}_{\odot }\,{\mathrm{SN}}^{-1}$)	${ \mathcal U }(0,20)$	4.01	${4.79}_{-0.19}^{+0.19}$
A⋆	Normalization of star formation law (106 M⊙ Gyr−1)	${ \mathcal U }(0,10)$	0.89	${0.79}_{-0.13}^{+0.14}$
αあるふぁ	Power-law index of star formation law	${ \mathcal U }(0,2)$	0.82	${0.72}_{-0.07}^{+0.09}$
FeIa	Fe yield from Type Ia SNe	${ \mathcal U }(0,0.9)$	0.80	${0.58}_{-0.03}^{+0.03}$
expCII	Exponent of C yield from CCSNe	${ \mathcal U }(0,2)$	1.0	${1.32}_{-0.03}^{+0.03}$
normMgII	Normalization of Mg yield from CCSNe	${ \mathcal U }(0,2)$	1.0	${1.41}_{-0.09}^{+0.10}$
normCaII	Normalization of Ca yield from CCSNe	${ \mathcal U }(0,0.5)$	0.01	${0.24}_{-0.05}^{+0.05}$
normCAGB	Normalization of C yield from AGB stars	${ \mathcal U }(0.4,5)$	0.60	${1.98}_{-0.36}^{+0.44}$
normBaAGB	Normalization of Ba yield from AGB stars	${ \mathcal U }(0,1)$	0.33	${1.08}_{-0.21}^{+0.31}$
meanBaAGB	Mass of peak Ba yield from AGB stars (M⊙)	${ \mathcal N }(2,0.25)$	2.0	${2.80}_{-0.16}^{+0.16}$

Note. For all parameters, the best-fit values are reported as the median (50th percentile) values, with uncertainties based on the 16th and 84th percentiles.

Download table as:  ASCIITypeset image

For all free parameters except meanBa_AGB, we assumed uniform priors with lower limits at 0 to avoid unphysical negative values and upper limits chosen based on the range of parameters determined by Kirby et al. (2011). For the parameter meanBa_AGB, which dictates the mass of the AGB stars that produce the most barium, we assume a normal prior with a mean of 2 M_⊙ and standard deviation 0.5 M_⊙. This is because meanBa_AGB describes when the s-process begins to contribute meaningfully to the abundance of barium: low meanBa_AGB means that lower-mass AGB stars produce most of the s-process barium, so [Ba/Fe]_s begins to increase at later times and higher [Fe/H]. As a result, any measurements of [Ba/Fe]_s at low [Fe/H] will have outsized leverage on the value of meanBa_AGB. We therefore enforce a Gaussian prior to keep this parameter within physically reasonable limits. Initial values of the parameters were chosen by performing a simple linear optimization of $-\mathrm{ln}L$. We sampled 10⁶ steps using 32 ensemble members or "walkers" initialized about these values, discarded the first 10⁴ "burn-in" steps, and used the remaining steps to sample the posterior distribution of the parameters.

We use the GCE model described in the previous section to fit the stellar abundances of Sculptor dSph. For Mg, Si, Ca, and C, we simultaneously fit our model to both the medium-resolution abundances from DEIMOS (Table 2) and the high-resolution DART abundances (Hill et al. 2019; North et al. 2012) (filled blue and empty orange points, respectively, in Figure 5). When fitting [Fe/H], we use only the DEIMOS yields because the DART sample is significantly smaller than the DEIMOS sample (N_DART = 89 compared to N_DEIMOS = 376). When fitting s-process Ba abundances, we use only the DART yields because the statistical correction used to remove the r-process contribution (Section 5.3) was based on the DART yields of Eu and Ba.

Zoom In Zoom Out Reset image size

Figure 5 shows the best-fit abundance trends from the GCE model, illustrating that the model fits the stellar abundances reasonably well. We also plot the outputs from the best-fit GCE model, showing how the components of the galaxy change over time, in Figure 6. We note that while the rates of Type Ia SNe (dotted line in bottom panel) appear relatively low—particularly when compared to AGB stars at late times—this is because the average yield per Type Ia SN is much larger (typically by at least two orders of magnitude; see Figures 12 and 14) than the average yield per AGB wind. Despite their small numbers, Type Ia SNe dominate Sculptor's chemical evolution at late times.

Zoom In Zoom Out Reset image size

The corresponding cumulative SFH from this best-fit model is shown in Figure 7 as a thick black line; to give a rough sense of the uncertainty in this model, thin black lines represent random realizations of the posterior distribution. We find that Sculptor has an ancient stellar population: the best-fit SFH is a single burst of star formation with a relatively short duration of ∼0.92 Gyr. We note that chemical evolution models measure relative rather than absolute ages, so the exact location of our measured SFH on the time axis is uncertain. Bromm & Yoshida (2011) point out that the halos thought to be the hosts of the first galaxies are predicted to form ∼500 Myr after the Big Bang; we therefore assume our SFH begins when the universe is 500 Myr old.

Zoom In Zoom Out Reset image size

As shown in Figure 7, we also compare our results to previous measurements of dSph SFHs. Nearly all measurements of Sculptor's SFH have been based on CMD fitting. Da Costa (1984) and Monkiewicz et al. (1999) fit a few model isochrones of varying ages and metallicities to CMDs and found that most of the stars in Sculptor must be relatively ancient (∼13 ± 2 Gyr old). Dolphin (2002) reanalyzed the data from Monkiewicz et al. (1999) by interpolating over a grid of synthetic CMDs to get a "true model CMD" rather than fitting discrete isochrones, coming to the same conclusion that Sculptor must be "entirely ancient." On the other hand, de Boer et al. (2012) argued that Sculptor has a much more extended SFH than these previous measurements would indicate. Using a new technique to simultaneously fit the CMD and the MDF, de Boer et al. (2012) found that Sculptor had a continuous period of star formation with a duration of ∼6–7 Gyr (dashed pink line in Figure 7). Savino et al. (2018) extended this technique to include horizontal branch stars in fitting the CMD and MDF, finding an extended SFH with a prominent tail of star formation at younger ages.

However, these results were based on relatively shallow photometry that did not reach significantly below the main-sequence turnoff (MSTO). Other measurements from deeper photometry appear to confirm the original picture of Sculptor as an ancient galaxy that formed all of its stars in a short burst. Weisz et al. (2014) measured CMDs using the Hubble Space Telescope, obtaining photometry with ∼30% completeness at ∼2 mag below the MSTO. From these, they determined that Sculptor formed ∼90% of its stars >10 Gyr ago (dashed green line and shaded region in Figure 7). Similarly, Bettinelli et al. (2019) used DECam to obtain a CMD down to ∼2 mag below the MSTO and found that Sculptor had a single burst of star formation with a full width at half maximum of ∼2.2 Gyr (dashed purple line in Figure 7).

As Figure 7 shows, SFHs measured from GCE models are also ancient, in qualitative agreement with the results from deep photometry. Yet GCE models tend to produce shorter absolute star formation durations (dotted lines) than the photometrically derived SFHs (dashed lines). This discrepancy likely arises from the limitations of photometric methods; although photometry is excellent at determining absolute ages, as discussed in Section 1 the age resolution of CMD fitting degrades for old and metal-poor populations. Abundance-derived estimates of the star formation duration may therefore better resolve the relative spread in ages within the ancient population of Sculptor.

Our results are largely consistent with previous SFH measurements using GCE models. The SFH we derive in this work (solid black line in Figure 7) has a duration of 0.92 Gyr. This is slightly shorter than the SFH found by Vincenzo et al. (2016), who used a one-zone model and found that 99% of the stars in Sculptor dSph formed within the first 2.16 Gyr of its evolution (dotted orange line in Figure 7). However, their model only aimed to fit the stellar metallicity distribution function and did not include information about individual abundances. Kirby et al. (2011), on the other hand, used similar methods to our work and traced both the stellar MDF and several abundance ratios. They found that Sculptor finished forming stars within 1.1 Gyr; as Figure 7 shows, this SFH (dotted blue line) is entirely consistent with the random realizations of the posterior distributions from our model (black lines).

This work expands on that of Kirby et al. (2011) by using more data—they used only DEIMOS abundances and did not use any [C/Fe] or [Ba/Fe] abundance information. Our model incorporates s-process abundances, as well as additional abundances measured from the VLT DART survey. It also has additional free parameters beyond those used in Kirby et al. (2011), which we use to fit the analytic functions describing nucleosynthetic yields. The consistency between our results and those of Kirby et al. (2011) suggests that including these additional parameters in our model does not significantly impact our main results; instead, as we will discuss later in Section 5.1, these parameters provide additional useful information about predicted nucleosynthetic yield sets.

In this section, we discuss the simplifying assumptions on which our chemical evolution model depends, and their potential impact on our results. These assumptions can broadly be classified into three categories: assumptions inherent to the construction of the model, assumptions in the model inputs, and other potential sources of systematic errors. In Figures 8 and 9, we illustrate some of the effects of changing these assumptions.

Zoom In Zoom Out Reset image size

4.2.1. Model Construction

In constructing our simple one-zone model, we have made a number of inherent assumptions. For example, the primary assumption in our model is the "one-zone" assumption of instantaneous mixing. This approximation is reasonably well-founded for dSphs; Escala et al. (2018) found that in simulated dwarf galaxies, a well-mixed ISM due to turbulent metal diffusion successfully reproduces observed abundance distributions of dSphs.

Other model assumptions may have greater impacts on the measured SFH. For example, the parameterization of gas inflow (Equation (5)) strongly influences the shape of the SFH because star formation depends on the gas mass, which is in turn predominantly set by the gas inflow. Fortunately, the inflow parameterization is constrained by the metallicity distribution function. This is shown in Figure 8, which illustrates the output MDFs of models using other common inflow parameterizations. To compare the goodness-of-fit of these models, we compute the Akaike information criterion (AIC); a lower AIC implies less information loss, so a given model is a "better" fit than our fiducial model if Δでるた(AIC) = AIC_model – AIC_fiducial is positive. A constant gas inflow, given by the equation

$$\begin{eqnarray}&&{\dot{\xi }}_{j,\mathrm{in}}(t)={A}_{\mathrm{in}}\displaystyle \frac{{M}_{j,\mathrm{gas}}(0)}{{M}_{\mathrm{gas}}(0)}\end{eqnarray} \tag{ 12 }$$

where A_in is a free parameter, is entirely unable to match the MDF of Sculptor dSph (Δでるた(AIC) ≫ 0).

An exponentially declining inflow,

$$\begin{eqnarray}&&{\dot{\xi }}_{j,\mathrm{in}}(t)={A}_{\mathrm{in}}\displaystyle \frac{{M}_{j,\mathrm{gas}}(0)}{{M}_{\mathrm{gas}}(0)}\exp \left(\displaystyle \frac{-t}{{\tau }_{\mathrm{in}}}\right)\end{eqnarray} \tag{ 13 }$$

with free parameters A_in and τたう_in, performs better than the constant inflow model, but not as well as our fiducial model (Δでるた(AIC) = 3). As shown in the top panel of Figure 9, using the exponentially declining inflow produces a single burst of star formation that ends after 0.69 Gyr (green dashed line). This is shorter than our fiducial model SFH (total duration of 0.92 Gyr) by a factor of ∼25%, likely because higher gas inflow at early times accelerates star formation and more quickly depletes gas. While this demonstrates that the parameterization of the gas inflow may have a significant effect on the predicted shape and duration of the SFH, more complex parameterizations of the gas inflow are somewhat disfavored for Sculptor dSph. As noted by de Boer et al. (2012) and Escala et al. (2018), among others, a bursty SFH might produce a wider spread in [αあるふぁ/Fe] at fixed [Fe/H] than observed in Sculptor (although see, e.g., Marcolini et al. 2008).

Zoom In Zoom Out Reset image size

Our model also ignores environmental effects like tidal or ram pressure stripping, which may cut off gas inflows and/or contribute to the removal of gas. Simulations (e.g., Kazantzidis et al. 2017) have shown that these environmental effects can, depending on galaxy orbital parameters, contribute significantly to gas loss in dSphs. To test this, we apply a simple model of ram pressure stripping by adding a constant to the gas outflow (Equation (6)). This parameterization is similar to that used in the analytic model of Kirby et al. (2013), who found that ram pressure stripping successfully reproduced the metallicity distribution function of Sculptor dSph. Following this parameterization, we apply an additional constant outflow starting at [Fe/H] = −1.5 (based on the best-fit parameterization found by Kirby et al. 2013, who found that ram pressure stripping in Sculptor began at [Fe/H] < −1.41):

$$\begin{eqnarray}&&{\dot{\xi }}_{j,\mathrm{out}}(t)=\displaystyle \frac{{M}_{j,\mathrm{gas}}(t)}{{M}_{\mathrm{gas}}(t)}\left[{A}_{\mathrm{out}}\left({\dot{N}}_{\mathrm{II}}+{\dot{N}}_{\mathrm{Ia}}\right)+{C}_{\mathrm{ram}}\right].\end{eqnarray} \tag{ 14 }$$

Here, we allow C_ram to be a free parameter setting the removal of gas due to ram pressure stripping. We fit it along with the other parameters in our GCE model, ⁹ finding a best-fit value of ${C}_{\mathrm{ram}}={2.76}_{-1.77}^{+1.53}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The resulting SFH is shown as a orange dotted line in the top panel of Figure 9, which is almost exactly the same as the fiducial model; the overall star-forming duration of 0.90 Gyr is ∼2% different from that predicted by the fiducial model and has almost the same goodness-of-fit (as measured by the AIC).

There are a number of other model assumptions that may affect our measurement of the SFH. For example, we do not account for reionization, which may heat infalling gas and delay star formation. Furthermore, our parameterization of the outflow rate as linearly proportional to the supernova rate (Equation (6)) assumes that supernovae are the only factors in determining outflow rates and that all supernovae (both Type Ia and core-collapse) contribute equally to outflows. However, recent hydrodynamic simulations indicate that other factors, such as a galaxy's gas fraction (e.g., M. Orr et al. 2021, in preparation) and the clustering of supernovae (e.g., Fielding et al. 2018), can strongly affect whether supernovae are able to produce galactic outflows. Equation (6) also assumes that the galactic gravitational potential remains constant, but dark matter accretion, environmental effects, and stellar feedback might all affect the underlying gravitational potential of a galaxy. Fully addressing these assumptions would require a more sophisticated model, which we defer to future work.

4.2.2. Model Inputs

4.2.3. Additional Systematics

As described in Section 2.3, the nucleosynthetic yields from CCSN and AGB models can vary widely. We aimed to work around these uncertainties by parameterizing the yields with analytic functions (Appendix). For a number of these yields we defined free parameters, allowing the model to vary the yields in order to best match the observed abundance trends. By examining the behavior of these free parameters, we can determine whether certain yield sets are preferred over others.

Among the CCSN yields, we varied the yields for C, Mg, and Ca. The best-fit yields, shown as blue lines in Figure 13, are a better match to the yields compiled by Nomoto et al. (2013) rather than those from the models of Limongi & Chieffi (2018). There are a number of differences between these two yield sets that could contribute to this discrepancy. At least one major difference is in the assumed explosion "landscape" (i.e., what masses of progenitor stars explode): Limongi & Chieffi (2018) assumed that all stars with initial masses above 25 M_⊙ implode, so that no yields are produced from explosive nucleosynthesis. As a result, many of their predicted CCSN yields (the dotted lines in Figure 13) are zero for high-mass progenitors. Our model, which finds nonzero yields at masses >25 M_⊙ (i.e., extremely prompt nucleosynthesis) for the parameterized elements C, Mg, and Ca, is therefore more consistent with the Nomoto et al. (2013) yields, which assume that all massive progenitor stars explode. This may suggest that an explosion landscape different from the simple one assumed by Limongi & Chieffi (2018) is needed to match observations (see, e.g., Griffith et al. 2021).

Our GCE model also exhibits other features—high Mg yields and relatively low Ca yields at low progenitor masses—in the best-fit yields that appear to be consistent with the Nomoto et al. (2013) yields. These features could result from a number of model assumptions. For example, Limongi & Chieffi (2018) evolve massive stars from pre-main sequence to pre-supernova, fix the explosion energy such that all exploding models produce exactly 0.1 M_⊙ of Fe, and include the effects of rotation. Nomoto et al. (2013), on the other hand, evolve massive stars from pre-supernova to explosion, assume fixed supernova energy E_SN = 10⁵¹ erg, and do not include rotation. Determining the exact reasons why our model appears to agree with the C, Mg, and Ca yields of Nomoto et al. (2013) is beyond the scope of this paper. We also note that our result differs slightly from the result of Nuñez et al. (2021, in preparation), who independently infer CCSN nucleosynthetic yields for a number of elements using damped Lyαあるふぁ systems. They find that while the C yields from Nomoto et al. (2013) ¹⁰ are consistent with the observed [C/Fe] ratio, the Limongi & Chieffi (2018) yields are typically more consistent with other observed abundance ratios, including [C/O].

For the AGB yields, we varied the yields of C and Ba, shown in Figure 14. The best-fit yields (blue lines) are more consistent with the Stromlo AGB yields from Karakas & Lugaro (2016) and Karakas et al. (2018), which predict enhanced C and Ba values at low progenitor masses. The AGB yields predicted by the FRUITY models (Cristallo et al. 2015) do not predict such a large enhancement. The large C and Ba yields may result from the Stromlo models including a deeper third dredge-up mixing, which brings more He-shell material to the stellar surface. This mixing predominantly affects elements made by neutron capture such as barium, along with the products of partial He-shell burning such as carbon.

Manganese and nickel are both iron-peak elements predominantly produced in Type Ia SNe. As discussed in Section 3.1, Mn is particularly sensitive to the physics of Type Ia SNe. The production of the only stable isotope of manganese, ⁵⁵Mn, depends strongly on the density—and therefore the mass—of the progenitor white dwarf (Seitenzahl et al. 2013; Seitenzahl & Townsley 2017). White dwarfs near the Chandrasekhar mass (M_Ch ≈ 1.4 M_⊙) are expected to produce solar or supersolar [Mn/Fe] and [Ni/Fe], while sub-M_Ch models tend to produce subsolar [Mn/Fe] and [Ni/Fe] (see, e.g., Figure 12).

Kirby et al. (2019) and de los Reyes et al. (2020) found that in Sculptor dSph, subsolar [Ni/Fe] and [Mn/Fe] abundances indicate that sub-M_Ch Type Ia SNe likely dominate. Both found that the Ni and Mn yields in Sculptor were most consistent with ∼1 M_⊙ sub-M_Ch models from Leung & Nomoto (2020). However, these comparisons were based on an analytic model that made a number of simplifying assumptions (e.g., that CCSNe are the only nucleosynthetic sources at early times, that CCSN yields are metallicity-independent, that the only contributions to stellar abundances are from supernovae). Our GCE model can be used to make more sophisticated comparisons.

In Figure 10, we use the best-fit parameters from Table 3 and vary the Type Ia yields of Mn and Ni. We find that the observed [Mn/Fe] and [Ni/Fe] trends are consistent with sub-M_Ch models, supporting the hypothesis that sub-M_Ch Type Ia SNe likely dominate in Sculptor dSph. In particular, the Leung & Nomoto (2020) model of a ∼1 M_⊙ CO white dwarf with a helium shell appears to best fit our data. This is in broad agreement with the findings of Kobayashi et al. (2020), who used a one-zone chemical evolution model to study dSphs and found that a significant fraction of sub-M_Ch Type Ia SNe are needed to reproduce the observed iron-peak abundances. Similar results have also been obtained in the Milky Way (see, e.g., Palla 2021, and references therein).

Zoom In Zoom Out Reset image size

Although we have attempted to pick a representative subset of Type Ia models, there are many other yield sets that we do not consider here. For the models that we have used in our comparison, we have included the effect of metallicity-dependent yields; however, our GCE model does not make any assumptions about the masses of Type Ia SNe that explode, so we are unable to include mass-dependent yields. For simplicity, we also do not consider the effects of multiple simultaneous channels of Type Ia SNe. In contrast, see, e.g., Kobayashi et al. (2020), who found that Type Iax SNe, potentially the pure deflagrations of hybrid C+O+Ne white dwarfs, are required to match the observed [Mn/Fe] abundances. Furthermore, we do not attempt to correct the observed Mn or Ni abundances for the effects of stellar atmospheres not being in local thermodynamic equilibrium (LTE). These non-LTE effects may be especially significant for manganese (e.g., Bergemann et al. 2019).

As shown in Figure 5, our simple GCE model is largely able to capture the behavior of the barium abundances produced by the slow neutron-capture process [Ba/Fe]_s . However, as discussed in Section 2.1, the rapid neutron-capture process (r-process) produces a non-negligible amount of barium. It also produces the majority of europium. Because our initial model does not include any r-process contribution, it significantly underpredicts the total Ba and Eu abundances (solid cyan lines in Figure 11).

Zoom In Zoom Out Reset image size

Despite recent work identifying neutron star mergers (NSMs) as a key site of the r-process (see, e.g., Rosswog et al. 2018; Shibata & Hotokezaka 2019, and references therein), the details of r-process nucleosynthesis are still uncertain. A full analysis of r-process nucleosynthesis is beyond the scope of this work, but our GCE model can place some constraints on r-process timescales and yields. In particular, we can constrain the rough timescale of r-process nucleosynthesis: does the r-process primarily occur in prompt events (such as high-mass or rapidly rotating CCSNe) or relatively more delayed events (such as NSMs)?

To test this, we use our best-fit GCE model (parameters from Table 3) and modify the input yields to simulate prompt or delayed r-process events. To model prompt r-process events, we add a contribution to the CCSN yields for Ba and Eu. We use Ba yields predicted by Li et al. (2014) (see the Appendix for details) and assume a universal r-process ratio of [Ba/Eu] ∼ −0.7 (Sneden et al. 2008) to compute Eu yields. To model delayed r-process events, we assume a relatively steep NSM DTD in the form of a t^−1.5 power law (e.g., Côté et al. 2019; Simonetti et al. 2019). Specifically, we initially adopt the following NSM DTD from Simonetti et al. (2019):

$$\begin{eqnarray}\begin{array}{c}{{\rm{\Psi }}}_{{\rm{Ia}}}=\left\{\begin{array}{ll}0 & {t}_{{\rm{delay}}}\lt 0.01\,{\rm{Gyr}}\\ ({10}^{-4}\,{{\rm{Gyr}}}^{-1}\,{M}_{\odot }^{-1}){\left(\displaystyle \frac{{t}_{{\rm{delay}}}}{{\rm{Gyr}}}\right)}^{-1.5} & {t}_{{\rm{delay}}}\geqslant 0.01\,{\rm{Gyr}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 15 }$$

We assume that each NSM produces the same amount of r-process elements: M_Ba = 2.3 × 10⁻⁶ M_⊙ and M_Eu = 2.3 × 10⁻⁷ M_⊙, from Li et al. (2014).

In the left panels of Figure 11, we plot the results of this test, showing the model predictions for [Ba/Mg] and [Eu/Mg] as a function of [Fe/H] assuming different r-process parameterizations: no r-process contributions (solid cyan line), contributions from either prompt or delayed channels (solid blue and purple lines, respectively), and contributions from both panels (solid black line). We plot [X/Mg] rather than [X/Fe] because Type Ia SNe contribute significantly to Fe at late times and could complicate our interpretation. Observations from both DEIMOS and DART are plotted for comparison. A successful model should be consistent with the observed trends: subsolar [Ba/Mg] below [Fe/H] ≲ −1.5 that increases to near-solar at higher metallicities, and near-solar [Eu/Mg] that is roughly constant as a function of metallicity.

We first consider whether prompt r-process events alone can reproduce these observed trends. Our initial test (left panels of Figure 11) finds that the prompt r-process channel (solid blue line) appears to significantly underpredict both [Ba/Mg] and [Eu/Mg]. This may suggest that our assumed CCSN yields of Ba and Eu (respectively from Li et al. 2014 and Cescutti et al. 2006) are too low. We therefore attempt to produce better consistency with observations by increasing the CCSN yields of Ba and Eu, as shown by the dotted lines in the right panels of Figure 11. In order to better match the Ba abundances, the r-process yields from typical CCSN-like events must be drastically increased by a factor of 20. This may correspond physically to some combination of increasing the rate of CCSN-like r-process events and increasing the Ba and Eu yields expected from these events. This enhancement of prompt r-process contributions is able to increase the predicted [Ba/Mg] trend at low metallicity, making it consistent with observations at −2.5 ≲ [Fe/H] ≲ −2 (blue dotted line in the upper right panel of Figure 11), but it underpredicts the observed [Ba/Mg] at −2 ≲ [Fe/H] ≲ −1.5.

Additionally, Frebel et al. (2010) measured abundances in the most metal-poor star in Sculptor and found abundance ratios of [Fe/H] ∼ −3.8 and [Ba/Mg] ≲ −1.58. This is well below the [Ba/Mg] ∼ −1 at low metallicities predicted by the prompt r-process channel with enhanced yields. Since r-process events are likely rare (e.g., Ji et al. 2016), it is possible that this star simply formed before any prompt CCSN-like events could produce barium. However, the prompt channel alone also underestimates the observed [Eu/Mg] trend (blue dotted line in the lower right panel of Figure 11), despite the drastic enhancement in r-process yields. We therefore tentatively conclude that prompt CCSN-like r-process events alone may not be able to explain the observed trends, and at least some contribution from a delayed NSM-like process is likely needed.

We now consider the delayed r-process channel. The left panels of Figure 11 show that the delayed r-process channel (solid purple lines) is more consistent with observations than the prompt r-process channel (solid blue lines). The delayed channel also appears to dominate r-process nucleosynthesis; prompt r-process events do not contribute significantly to the combined prompt+delayed model (solid black lines). Yet this delayed r-process model, based on the Simonetti et al. (2019) NSM DTD (Equation (15)), underpredicts the observed [Eu/Mg] trend. It also slightly overpredicts [Ba/Mg] at low [Fe/H] and underpredicts [Ba/Mg] at −2 ≲ [Fe/H] ≲ −1.

In order to better match observations, we increase the normalization of the delayed r-process contributions by a factor of five. Physically, this can be done either by increasing the normalization of the Simonetti et al. (2019) NSM DTD or by increasing the Ba and Eu yields expected from individual NSMs. Both of these modifications are plausible because both the NSM DTD and the nucleosynthetic yields from NSMs are relatively uncertain. We also increase the NSM t_min from 10 to 50 Myr. As shown by the dashed lines in the right panels of Figure 11, this allows the delayed r-process channel to begin contributing to r-process yields at later times than the CCSNe, so that [Ba/Mg] begins increasing at a higher [Fe/H]. This more-delayed NSM DTD produces a steeper [Ba/Mg] trend as a function of metallicity, which is more consistent with the observed trend. However, it also produces a steeper [Eu/Mg] trend as a function of metallicity, which is less consistent with the observed flat [Eu/Mg] trend (Skúladóttir et al. 2019). We conclude that relatively delayed (with minimum delay times ≳50 Myr) r-process events alone may not be able to reproduce observed abundance trends in Sculptor, and that a combination of prompt and delayed events is needed.

Our result is consistent with previous results. Duggan et al. (2018) found that in multiple Local Group dSphs, the positive trend of [Ba/Fe] as a function of [Fe/H] appears to be steeper than the positive trend of [Mg/Fe] as a function of [Fe/H]. They argued that the primary r-process source of Ba must therefore be delayed relative to CCSNe (the primary source of Mg) in order to produce the [Ba/Fe] trend in dwarf galaxies. In a separate analysis of Sculptor dSph, Skúladóttir et al. (2019) use the flat [Eu/Mg] trend to argue the opposite: that the primary source of Eu must not be significantly delayed relative to the primary source of Mg. As our GCE model shows, both prompt and delayed r-process channels may be needed to reproduce the observed Ba and Eu trends in Sculptor. This agrees with investigations of the r-process in the Milky Way and in the universe (e.g., Matteucci et al. 2014; Cescutti et al. 2015; Wehmeyer et al. 2015; Côté et al. 2019; Siegel et al. 2019), which find that a combination of delayed NSMs and prompt CCSN-like events can successfully reproduce the observed trend and scatter of [Eu/Fe] in the Milky Way.