© ESO, 2016

1. Introduction

Over the last decade new observing facilities have allowed accurate measurements of elemental abundances in a large number of individual stars not only in the Milky Way but also within its satellites, specifically the Local Group dwarf spheroidal galaxies (dSphs) and ultra-faint dwarfs (UFDs).

These measurements cover galaxies of very different mass and star formation histories providing crucial constraints on models of galaxy evolution (see Tolstoy et al. 2009, for a review). Not only should these models reproduce the observed dynamical properties of galaxies, but they should also be able to account for their chemical properties. Unfortunately, these chemical constraints are often neglected despite the fact that they can lead to erroneous conclusions.

While the mean metallicities at a given mass or luminosity is normally close to observations, the abundance ratios and their corresponding scatter is often poorly reproduced. This discrepancy suggests difficulties in accurately reproducing the number of supernovae or their associated feedback which are the dominant source of observable metals.

The yield of α-elements in Type II supernovae (SNeII) strongly depends on the progenitors mass. Tsujimoto et al. (1995) and Woosley & Weaver (1995) found that the ratio [ Mg / Fe ] in SNII ejecta decreases by two orders of magnitude between stars with masses above 50 M_⊙ and stars lighter than 10 M_⊙. As the lifetime of massive SNII ranges from ~3 to 30 Myr, one naturally expects that the local interstellar medium (ISM) polluted by massive stars is different from that polluted later by lighter stars even if the stars originated in the same star-forming region. It was pointed out long ago by Audouze & Silk (1995) that a small number of exploding supernovae with varied masses sampling an initial mass function (IMF) introduce some scatter in abundance ratios. Low scatter in the observed abundances must then be a sign of continual pollution by successive supernovae and/or the sign of efficient mixing, which must be reproduced in simulations.

In the Milky Way, above [ Fe / H ] > −3.5, using a semi-analytical model, Karlsson (2005) found that the probability of finding a star enriched by less than ten supernovae is very low, explaining a possible mixing of the different yields of stars with varied abundances. However, despite this mixing, stars with [ Mg / Fe ] < 0 are still predicted to exist (Karlsson & Gustafsson 2005). The problem naturally worsens for dwarf galaxies, where the rate of exploding supernovae is much lower. To efficiently mix the ISM, not only dynamical effects like ISM turbulence but also the motion of stars before they explode need to be considered. Both are present with a certain degree of precision inside N-body simulations. Recently, focusing on a small volume (~32 pc³) Feng & Krumholz (2014) showed that the star-to-star variation in abundances within an open cluster may be considerably reduced owing to the mixing of inhomogeneous gas during the process of star formation. Similar conclusions have been recently obtained by Petit et al. (2015).

Different chemo-dynamical simulations have been performed since the 1990s to study the metal enrichment in cosmological contexts (Steinmetz & Mueller 1994; Mosconi et al. 2001; Lia et al. 2002; Scannapieco et al. 2005; Oppenheimer & Davé 2008; Wiersma et al. 2009; Few et al. 2012), in spiral galaxies (Raiteri et al. 1996; Berczik 1999; Friedli et al. 1994; Friedli & Benz 1995; Lia et al. 2002; Samland & Gerhard 2003; Stinson et al. 2006; Kobayashi & Nakasato 2011; Tissera et al. 2012; Few et al. 2012; van de Voort et al. 2015), in ellipticals (Kawata & Gibson 2003; Kobayashi 2004, 2005; Martínez-Serrano et al. 2008), and in dwarf irregular, spheroidal or elliptical galaxies (Carraro et al. 2001; Ricotti & Gnedin 2005; Marcolini et al. 2006, 2008; Kawata et al. 2006; Valcke et al. 2008; Revaz et al. 2009; Okamoto et al. 2010; Sawala et al. 2010; Schroyen et al. 2011; Revaz & Jablonka 2012; Schroyen et al. 2013). However, while a few authors present [ α/ Fe ] vs. [ Fe / H ] diagrams, very little discussion exists on the scatter of α-elements.

Interestingly, it appears that in early works the scatter in abundance ratios was not an issue, because of the artificial mixing induced by the poor resolution used at that time (Raiteri et al. 1996; Berczik 1999). As pointed out by Mosconi et al. (2001), increasing the resolution of smoothed-particle hydrodynamics (SPH) models increases the scatter in the chemical properties. This is indeed witnessed in recent high-resolution models of dSphs (Sawala et al. 2010; Revaz & Jablonka 2012). In their zoom-in simulation of a Milky Way-like galaxy, van de Voort et al. (2015) predicted a subsolar [ Mg / Fe ] population that is not observed. Similar cases also exist in elliptical models (Kawata & Gibson 2003; Kawata et al. 2006; Marinacci et al. 2014).

This scatter in abundance ratios can be reduced by fixing the size of the region in which metals are expelled (Kobayashi & Nakasato 2011), because on average, particles receive more ejecta with varied yields. A similar mixing is achieved by increasing the number of SPH neighbours as the volume that receives ejecta is increased, according to Revaz & Jablonka (2012). Wiersma et al. (2009) introduced a natural SPH smoothing over the volume into which metals are ejected.

An alternative solution consists in inserting a diffusion mechanism for the heavy elements through the ISM. Following a suggestion from Groom (1997), Carraro et al. (1998) added a diffusive equation in their SPH implementation. However, the impact of this term on the distribution of abundances was not studied. Martínez-Serrano et al. (2008) proposed an implementation of the classical diffusion equation in their SPH code. A similar approach has been followed by Greif et al. (2009) but only applied to simulation of supernovae remnants. This technique was improved by Shen et al. (2010) in the context of the intergalactic medium, where the diffusion coefficient results from a turbulent mixing model. Similar diffusivity has been used in the spiral galaxy simulations of Brook et al. (2012). and has resulted in a considerable reduction of the oxygen abundance dispersion, thereby overcorrecting the problem.

It is clear that the impact of the different possible schemes used to distribute stellar ejecta into the ISM needs to be thoroughly examined. This is particularly important for the SPH technique, where metals issued from a supernovae are usually distributed among neighbouring gas particles, using a weight proportional to the SPH kernel. Consequently particles at different radii receive differing amount of metals which may generate an artificial scatter. Other schemes could help to reduce this scatter, however, the final impact on stellar abundances is not straightforward and detailed examinations are needed.

A supplementary source of scatter comes from the modelling of the initial mass function (IMF) as resolutions increase. Numerical techniques usually assume that a stellar particle represents a single stellar population (SSP). However, with the mass in the highest resolution simulations approaching that of small clusters, serious concerns arise. These questions include: how to deal with stellar particles of masses smaller than what is required to ensure a full sampling of the IMF? and What happens if the stellar particle mass goes below 1000 M_⊙, approaching the mass of the most massive stars? It is crucial to answer these questions and estimate the extent to which a poor IMF sampling may bias the final stellar abundances and the star formation history.

The aim of this paper is to examine how chemo-dynamical N-body simulations can reproduce the observed scatter of [α/Fe] at galactic scales in dSph and Milky Way-like galaxies. We start by presenting facts concerning the scatter observed in the stellar abundance ratios. We then describe the schemes used in our code GEAR, including recent improvements concerning the IMF sampling, spreading of elements, and mixing schemes. We apply our method to very simple cases, where only one or two supernovae explode and pollute an initially homogeneous box. In addition to the verification that our code correctly reproduces the Sedov-Talor solution, these tests allow us to understand the bias that metal spreading techniques induce and to test the efficiency of mixing schemes. In a following step, we simulate Fornax and Sextans like dSphs. Based on the prediction of star formation histories and abundance scatters, the limitations of different IMF sampling schemes and the necessity of mixing schemes are discussed. Finally, we show that our choice of parameters for dSphs also reproduces the star formation rate and α-element abundances in a Milky Way-like galaxy.

2. Observational facts

As our aim is to properly reproduce the abundance ratio scatter seen in observations, we start by presenting observational facts. Figure 1 shows the compilation of 25 different samples in the Milky Way and its satellites. The relation between [ Mg / Fe ] and [ Fe / H ] arises from more than 1700 individual stars with abundances derived from high-resolution spectroscopy (R ≥ 20 000). All values have been scaled to the solar abundances of Asplund et al. (2009). Three major conclusions can be drawn: (i) Stars with [ Fe / H ] < −2.5 form a plateau around [ Mg / Fe ] ~ 0.4 with an error weighted standard deviation of about 0.2−0.3 dex. This is illustrated by the yellow shaded region of Fig. 1 corresponding to the 1σ deviation around the mean [ Mg / Fe ] value, computed for all stars with a metallicty lower than −2.5. (ii) Very few metal-poor stars ([ Fe / H ] < −2.5) are found at subsolar [Mg/Fe] values. (iii) With the exception of the peculiar cases of ET0381 (Jablonka et al. 2015) and SDSSJ0018-0939 (Aoki et al. 2014), none of these subsolar cases are found at [Mg/Fe] below −0.5. This three major features guide our discussion and choices in the following.

3. The code GEAR and its improvements

Our tests and analyses are conducted with the code GEAR, which was developed by Revaz & Jablonka (2012). It is a fully parallel chemo-dynamical Tree/SPH code based on Gadget-2 (Springel 2005). The code GEAR has been successfully used in the context of dwarf spheroidal galaxies, ensuring good performance, numerical convergence, conservation of the total energy budget and reproduction of the main observable properties of dSphs. In addition to the hydrodynamics of the gas, GEAR includes the complex treatment of baryonic physics, namely gas cooling, star formation, chemical evolution, and Type Ia (SNeIa) and Type II supernova feedback. Numerical codes are constantly improved and we list and comment hereafter the recent changes compared to the initial version, excluding the initial mass function, and mixing schemes which we discuss in more depth in Sects. 4 and 5.

3.1. Pressure-entropy formulation of SPH

We have upgraded the SPH formulation of GEAR using the pressure-entropy formulation described in Hopkins (2013). This new formulation alleviates difficulties traditional methods have had in the treatment of fluid mixing instabilities. This is particularly useful in the context of galaxy formation and evolution to correctly treat the dynamics of the ISM when, for example, the galaxy is subject to ram pressure stripping (Nichols et al. 2015).

3.2. Individual and adaptive timesteps

The individual and adaptive timesteps scheme now precisely follows the algorithm proposed by Durier & Dalla Vecchia (2012), which extends the timestep limiter of Saitoh & Makino (2009). The algorithm can be summarised as:

• Timestep limiter: each particle ensures that its timestep is shorter than or equal to a multiple (here taken as four) of any neighbouring particle.

• Timestep update: an inactive particle (one whose timestep does not coincide with the current one) becomes instantaneously active if it is touched by the feedback energy of any other particle.

• Timestep criterion: the default acceleration timestep criterion of Gadget-2 involving the ratio between the gravitational softening and the acceleration is supplemented by a restriction including the SPH smoothing length (see Eq. (B3) in Durier & Dalla Vecchia 2012).

These new features are necessary to correctly reproduce the blast waves of a single supernova explosion within high-resolution simulations (see Sect. 7.2).

3.3. Artificial viscosity

The artificial viscosity is an important ingredient in SPH methods designed to correctly capture shocks. Instead of using the default Gadget-2 artificial viscosity based on Monaghan (1997), we now use the original formulation of Monaghan & Gingold (1983) with the Balsara switch f_ij (Balsara 1995) which alleviates the spurious transport of angular momentum in the presence of shear flows. The viscosity between two particles i and j, Π_ij is therefore written as ${\mathrm{\Pi }}_{\mathit{ij}}\mathrm{=}\mathrm{-}\frac{\mathrm{(}\mathit{\alpha }{\mathit{c}}_{\mathit{ij}}\mathrm{+}\mathit{\beta }\hspace{0.17em}\hspace{0.17em}{\mathit{\mu }}_{\mathit{ij}}\mathrm{)}{\mathit{\mu }}_{\mathit{ij}}}{{\mathit{\rho }}_{\mathit{ij}}}\hspace{0.17em}{\mathit{f}}_{\mathit{ij}}\mathit{,}$(1)where α and β are two dimensionless parameters determining the viscosity strength, c_ij their averaged sound speed and μ_ij is given by ${\mathit{\mu }}_{\mathit{ij}}\mathrm{=}\frac{{{v}}_{\mathit{ij}}\mathrm{·}{{r}}_{\mathit{ij}}\hspace{0.17em}{\mathit{h}}_{\mathit{ij}}}{{\mathrm{\left|}{{r}}_{\mathit{ij}}\mathrm{\right|}}^{\mathrm{2}}\mathrm{+}\mathit{ϵ}\hspace{0.17em}{\mathit{h}}_{\mathit{ij}}^{\mathrm{2}}}\mathit{,}$(2)where r_ij and v_ij are the differences in position and velocity of the two particles, respectively, h_ij is the mean of their SPH smoothing lengths, and ϵ = 0.01 is a small parameter to avoid numerical divergence when the particles are close. This formulation is similar to that implemented in Gadget-2 if we impose the parameter β = 3 / 2 α and replace μ_ij by ω_ij = v_ij·r_ij/ | r_ij | as discussed in (Springel 2005). We show in Section 7.2 that this modification improves the reproduction of the blast waves in Sedov-Taylor experiments.

To prevent spurious dissipation even far away from shocks, Morris & Monaghan (1997) and later Rosswog et al. (2000) introduced particle-based time dependent viscosity coefficients α_i(t). The value of α_i(t) depends on a decay timescale and a source term, making it maximal in the presence of a shock and minimal away from it. This is further improved with the approach proposed by Cullen & Dehnen (2010), which avoids some problems related to the previous methods, such as the damping of sound waves or the delay between the peak of the α_i(t) coefficients and the shock front. However, this last improvement occurs at unacceptable increase of CPU time for our purposes.

In the present version of GEAR, we use individual time dependent viscosity coefficients¹α_i(t) as proposed by Rosswog et al. (2000) and fix β_i(t) = 2 α_i(t).

4. Initial mass function sampling

In most of the chemo-dynamical codes, each stellar particle is considered a single stellar population (SSP) and its stellar mass distribution at birth follows the shape of a given IMF. As time passes, the energy and elements released by dying stars must be spread over the ISM. There are various ways of modelling this.

The simplest approach, which is still widely used in N-body simulations (see for example Rosdahl et al. 2015), is to calculate the total energy and metals that are released over time by a full stellar particle, given its mass and IMF. The energy and metals are then ejected in the ISM at once. At the opposite end, pure chemical models (e.g., Matteucci et al. 2009) precisely compute the equations of chemical evolution (Tinsley 1980). In this case ejecta are gradually released over a longer period of time, the longest of which corresponds to the lifetime of the least massive star that ends its life as a supernova (or AGB when winds are considered).

In GEAR, each newly formed stellar particle is treated along the pure chemical evolution model approach. Energy and metals released by the SNeII and SNeIa are spread over the nearest neighbouring particles (see Sect. 5 for more details). This procedure requires accurate computation of, for each stellar particle i and at every timestep, the number of exploding SNeIa (N_SNIa,i(t)) and SNeII (N_SNII,i(t)). These numbers are directly dependent on the choice of IMF and its implementation.

Hereafter, we explore three different techniques to numerically sample the IMF. We study their impacts on the chemical evolution of a galaxy, in particular on the scatter in abundance ratios which are discussed in Sects. 8 and 9.

4.1. Continuous IMF sampling (CIMFS)

The continuous IMF sampling (CIMFS) method (e.g., Kobayashi et al. 2000; Revaz & Jablonka 2012) involves discretising the Tinsley’s equations to obtain an analytic form of the number of Type II and Type Ia supernovae, N_SNIa,i(t,Δt) and N_SNII,i(t,Δt) (Poirier 2004; Revaz et al. 2009), which explode within a timestep Δt. In this formulation, both values are ‘real’ numbers that depend continuously on Δt as well as on the stellar particle mass considered, the smaller mass producing less supernovae. This method works well for large stellar particle masses. Though with the advent of high-resolution simulations, it can lead to odd situations where N_SNII,i or N_SNIa,i are much smaller than 1, i.e. fractions of a supernova explode during a timestep. This is equivalent to the dilution of a single supernova energy and elements over a large time interval. Precisely, the time interval over which one supernova fully explodes is the one for which N_SNIa,i(t,Δt) = 1 or N_SNII,i(t,Δt) = 1. Diluting these explosions over time has an obvious impact on the ISM both in terms of dynamics and chemical enrichment.

For timesteps imposed by the Courant-Friedrichs-Levy condition, the situation described above comes to a head for Type Ia SNe at medium resolution, where the mass of a stellar particle becomes smaller than 10⁵M_⊙. As an example, for a gas density of 0.1 atom / cm³ and a temperature of 10⁶ K, the timesteps are constrained to be of about 0.1 Myr. Assuming a Kroupa IMF (Kroupa 2001) and the Kobayashi et al. (2000) SNIa model, the typical number of type Ia supernovae exploding during this timestep is 0.002. This corresponds to a dilution of the explosion over 500 Myr. This problem is worse for smaller masses or timesteps. Because Type II SNe explode at a higher rate, they are less prone to this problem. However, for a stellar mass resolution of 10³M_⊙, the dilution of the Type II supernova ejecta can reach 50 Myr, a period larger than the longest SNII lifetime.

4.2. Random discrete IMF sampling (RIMFS)

In principle, one could randomly sample the IMF at particle creation to determine when and which stars would explode as supernovae. However, such an approach would require unpractical amounts of memory for any galaxy. An alternative is to follow a stochastic approach that reproduces the discretisation of the IMF without the requirement of storing information for each star.

At every timestep, one calculates, for each stellar particle i, the number N_SNx,i(t,Δt) of potentially exploding supernovae (x standing for Ia and II). The integer part of N_SNx,i(t,Δt) is assumed to explode as supernova. The remaining fractional value is compared to a random variable χ, taken from the uniform distribution over [ 0,1 ]. If $\mathit{\chi }\mathit{<}{\mathit{N}}_{\mathrm{SNx}\mathit{,i}}^{\mathrm{frac}}\mathrm{\left(}\mathit{t,}\mathrm{\Delta }\mathit{t}\mathrm{\right)}$, the final number of dying stars is ⌊ N_SNx,i(t,Δt) ⌋ + 1 (where ⌊·⌋ corresponds to the largest previous integer), otherwise it is ⌊ N_SNx,i(t,Δt) ⌋. In the following, we refer to this method as the random IMF sampling (RIMFS) method.

The drawback of this technique is that when the stellar particle mass is small, typically smaller than 10⁴M_⊙, the number of supernovae generated from the IMF is also small, resulting in Poissonian noise. This is particularly pronounced at the high-mass end of the IMF as massive stars are comparatively rare. This is illustrated in the top panel of Fig. 2, which shows the cumulative number of SNII and the number of exploding SNII per timestep issued from a 2048 M_⊙ stellar particle. As stellar masses are directly related to the chemistry via metal ejection, an artificial scatter in masses induces a scatter in abundance ratios.

4.3. Optimal discrete IMF sampling (OIMFS)

The second choice for the discretisation of the IMF is the so-called optimal IMF sampling (OIMFS), described by Kroupa et al. (2013) and originally based on the work of Weidner & Kroupa (2004). It aims to decrease the noise introduced by a random discrete sampling of the IMF, i.e. it avoids gaps in the initial stellar mass distribution.

Two important masses are introduced. First, ${\mathit{m}}_{\max }^{\mathit{\star }}$, the absolute upper mass limit for stars above which no stars may form. It defines the absolute maximum stellar mass in the IMF. The second mass is m_max, the maximal star mass of the actual sampled IMF, i.e. the maximal star mass in the SSP. The OIMFS formalism imposes that there is only one star of mass m_max, i.e. ${\mathrm{\int }}_{{\mathit{m}}_{\max }}^{{\mathit{m}}_{\max }^{\mathit{\star }}}\frac{\mathrm{\Phi }\mathrm{\left(}\mathit{m}\mathrm{\right)}}{\mathit{m}}\hspace{0.17em}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{1}\mathit{,}$(3)where Φ(m) /m dm is the number of stars in the mass interval [ m,m + dm ]. The lower stellar masses are determined iteratively by ensuring the interval between [ m_{i + 1},m_i ] contains only one star with a mass we set to be m_{i + 1}. The total mass of an SSP for a given IMF and a maximum stellar mass m_max is therefore ${\mathit{M}}_{\mathrm{SSP}}\mathrm{\left(}{\mathit{m}}_{\max }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{m}}_{\mathit{L}}}^{{\mathit{m}}_{\max }}\mathrm{\Phi }\mathrm{\left(}\mathit{m}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{m}\mathrm{+}{\mathit{m}}_{\max }\mathit{,}$(4)where m_L is the minimum stellar mass.

Pflamm-Altenburg & Kroupa (2006) developed an algorithm to sample an IMF following this scheme, releasing a set of C routines². This algorithm has been implemented in GEAR fixing the parameter ${\mathit{m}}_{\max }^{\mathit{\star }}$ to 50 M_⊙, which is the upper mass limit of our SNe nucleosynthesis tables. The mass M_SSP is given by the stellar particle mass and is directly related to the resolution of the simulation.

Figure 2 illustrates the reduction of noise in the IMF obtained with the OIMFS (b) compared to the RIMFS (a). The number of exploding SNeII per unit time is shown for a 2048 M_⊙ stellar particle.

In the RIMFS approach, the individual stellar masses are uncorrelated, allowing multiple supernovae within short periods of time or a long gap between two explosions. This results in a momentary change in the cumulative number of SNeII compared to the CIMFS scheme. With the OIMFS scheme (panel b), these deviations are erased, and the fit to the cumulative number of SNeII in the continuous IMF sampling is clearly improved.

The OIMFS scheme is however limited in its application by the direct relation given by Eq. (4) between the most massive star in the SSP m_max as well as the SSP mass M_SSP. This relation shown in Fig. 3. Clearly, if the resolution is too high, for example in the extreme case where M_SSP = 100 M_⊙, the maximal stellar mass considered is only 7 M_⊙, lower than the minimal mass of a supernova. A more reasonable choice of stellar particle mass is thus around M_SSP = 10⁴M_⊙. Further discussions are conducted in Sect. 8.3.

4.4. An implicit limit in the mass resolution?

From the above description, one sees that each method presents an implicit limitation in mass resolution: the CIMFS is no longer reliable below a mass resolution of about 10⁵M_⊙. For higher resolutions, the ability to discretise the IMF is lost and energy and metals are diluted over unrealistic time intervals. The OIMFS breaks down below 10⁴M_⊙ owing to its intrinsic formulation that sets up a maximal star mass lower than the maximal IMF mass. The RIMFS induces noise when the stellar particle mass resolution is below about 10⁴M_⊙.

One can argue that a stellar particle does not have to sample a full IMF by itself, but instead a complete IMF should result from different contributions of several stellar particles. In this sense, the RIMFS method correctly samples the IMF as long as enough mass/particles are considered and are close to each other. However, as gravity quickly decorrelates stellar particles in numerical simulations, this assumption is only correct as long as 10⁴M_⊙ of stellar particles occupy a volume corresponding to the zone where one particle injects energy and elements. This corresponds to the volume defined by the SPH radius (r_SPH = h_i) of the particle. Computing the ratio between r_SPH and the size r₄ of a region containing 10⁴M_⊙ is not trivial. To accomplish this, we have extracted those values directly from simulations taken from Revaz & Jablonka (2012), where a large range of mass resolutions were explored. r_SPH is taken as the minimum SPH radius among the stellar particles and r₄ is computed in the central regions of the galaxies. The ratio r_SPH/r₄ is displayed on Fig. 4 as well as the ratio r_SPH/r₅ where r₅ similarly defines a region containing 10⁵M_⊙.

Despite the scatter, it appears that bellow a stellar mass of 10³M_⊙, the SPH volume of a stellar particle contains fewer stars than that needed to correctly sample a full IMF. While this allow us to extend the resolution over one order of magnitude in mass when considering the RIMFS method, we still face a mass resolution limit.

We further discuss the impact of this mass limit onto the chemical evolution of galactic systems in Sect. 8.

5. Spreading of elements

For a given stellar particle, when a supernova explodes both the energy and the synthesised elements must be ejected into the surrounding gas. In practice, one wants to (i) choose N_ngb neighbouring gas particles j around a stellar particle i, and (ii) send to each of them a fraction α_ij of the total elements released by the particle i. The fraction α_ij can be written in a general form as ${\mathit{\alpha }}_{\mathit{ij}}\mathrm{=}\frac{{\mathit{X}}_{\mathit{ij}}}{{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{ngb}}}{\mathit{X}}_{\mathit{ik}}}\mathrm{·}$(5)Possible values for X_ij will be discussed below.

5.1. Setting the proper weighting

The most natural way of setting α_ij is to use the SPH formulation ${\mathit{X}}_{\mathit{ij}}\mathrm{=}\frac{{\mathit{m}}_{\mathit{j}}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}}{{\mathit{\rho }}_{\mathit{j}}}\mathit{,}$(6)where, in addition to the kernel function W, which depends on h_i (the SPH smoothing length), and r_ij (the distance between particles i and j), the amount of each element distributed to the particle j is weighted according to its volume V_j = m_j/ρ_j. A simpler alternative is to weight the kernel according to the particles mass${\mathit{X}}_{\mathit{ij}}\mathrm{=}{\mathit{m}}_{\mathit{j}}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(7)In both cases, the final metallicity inevitably traces the shape of the kernel. As shown in Sect. 7, using a simple isolated supernova explosion, this radial dependency of the ejecta generates an artificial abundance scatter in the gas. Its impact on the final stellar abundances is discussed in Sect. 8 in the context of dSph galaxies.

By default, the kernel function is a cubic spline (Monaghan & Lattanzio 1985). A simpler alternative consists in using a step function $\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathrm{\Gamma }}_{\mathit{ij}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1}& \mathrm{if} {\mathit{r}}_{\mathit{ij}}\mathrm{\le }{\mathit{h}}_{\mathit{i}}\mathit{,}\\ \mathrm{0}& \mathrm{otherwise,}\end{array}$(8)assuming that each particle within the SPH radius receive the feedback independent of its distance to the source. Under this Eqs. (6) and (7) then become ${\mathit{X}}_{\mathit{ij}}\mathrm{=}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}} \mathrm{and} {\mathit{X}}_{\mathit{ij}}\mathrm{=}{\mathit{m}}_{\mathit{j}}\mathit{,} \mathrm{respectively}\mathrm{.}$(9)Based on different simulations, we observe no statistically distinguishable difference in terms of star formation or final chemical properties, when using either the particle volume or the particle mass-weighting recipe. This agrees with the results from Tornatore et al. (2007) obtained in the context of galaxy clusters, where only a few tenths of dex difference in the final [Fe/H] of the intra cluster medium exists between the volume or mass-weighting schemes. Hereafter, we use the mass-weighting. However, the choice of the kernel is important in this case, as illustrated both by simple supernovae explosion experiments (Sect. 7) and more complex dSphs simulations (Sect. 8).

5.2. Choosing neighbouring particles

A natural way of choosing the number of neighbours affected by a supernova is to follow the SPH weighting scheme. This is our default choice where we have set N_ngb = 50. However, in the rigorous Lagrangian SPH formulation, such as that used in Gadget-2, the number of particles is defined as (Springel & Hernquist 2002; Hopkins 2013) $\mathit{N̅}{\mathrm{ngb}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{h}}_{\mathit{i}}^{\mathrm{3}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{{\mathrm{ngb}}_{\mathit{i}}}}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{orequivalently}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{M̅}\mathit{i}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{h}}_{\mathit{i}}^{\mathrm{3}}{\mathit{\rho }}_{\mathit{i}}\mathit{,}$(10)ensuring that a fixed mass is contained within a given volume. Here, N_ngb,i is the effective number of neighbouring particles (an integer number) and is not strictly equal to (a real number). It is important to emphasise that h_i as well as N_ngb,i are affected by the mass distribution around the particle via their kernel dependency. This is particularly important if the mass is distributed in a ring-like structure, like it is, for example, after a Sedov-Taylor explosion (see Sect. 7.2). In this case, the effective number of neighbouring particles may be much larger than expected. For larger values of N_ngb,i, each neighbouring gas particle receives a lower mass of elements and induces a bias.

An alternative to this scheme is to fix the size of the neighbouring region instead of setting a constant number of particles, where the ejecta is released into a given volume instead of a given mass. The radius of the volume can be fixed as in Kobayashi & Nakasato (2011) or set according to the propagation of the supernova blast wave which depends on the local properties of the gas. In the latter case, we have used the blast radius R_E (see Eq. (9) of Stinson et al. 2006) computed by Chevalier (1974), McKee & Ostriker (1977), corresponding to the size of a supernova bubble has reached when its internal pressure drops to that of the ambient interstellar pressure.

An advantage of this method is that it is independent of the resolution as the amount of mass affected by ejecta is approximatively constant. The drawback is a drastic increase in the CPU time, as, according to our tests, the number of neighbouring particles is found to be on average larger.

We discuss the effect of these methods further in Sects. 7 and 8.

6. Mixing of metals

In common SPH methods, star formation feedback transfers metals amongst the nearest neighbours where the metals remain for all time. Such a scheme ignores the natural mixing that occurs as a reslut of small-scale turbulence or other unresolved processes. In most cases the metal content of particles is assumed to be constant until enriched by the next supernova. It also results in an artificial scatter as particles at different distances from the source receive varying amount of metals. Furthermore, as the number of neighbours is typically resolution dependent, this scheme results in a strong resolution dependence with the scatter becoming worse at higher resolutions.

These problems have been extensively discussed in the literature (Carraro et al. 1998; Okamoto et al. 2005; Tornatore et al. 2007; Martínez-Serrano et al. 2008; Greif et al. 2009; Shen et al. 2010), and two solutions have emerged: (1) the introduction of a diffusion equation into the SPH scheme and (2) considering the metallicity to be a smoothly varying function.

6.1. Metal diffusion

Including a classical diffusion equation in SPH has been proposed by different authors (e.g. Carraro et al. 1998; Martínez-Serrano et al. 2008; Greif et al. 2009; Wiersma et al. 2009; Shen et al. 2010). Hereafter, we use the approximate solution of Greif et al. (2009). More precisely, we implemented its Eq. (5). which predicts the time evolution of any scalar quantity associated with a particle. This equation involves a diffusion coefficient D_i, which reads ${\mathit{D}}_{\mathit{i}}\mathrm{=}\mathrm{2}\hspace{0.17em}\mathit{d}\hspace{0.17em}{\mathit{\rho }}_{\mathit{i}}\hspace{0.17em}{\mathit{h}}_{\mathit{i}}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{i}}\end{array}\mathit{,}$(11)where h_i is the SPH smoothing length of the particle i, ρ_i its density, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{i}}\end{array}$ is an estimation of the RMS velocity at the position i (Klessen & Lin 2003), and d is a free parameter used to calibrate the diffusion. Diverging slightly from Greif et al. (2009) we use the following estimation of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{i}}\end{array}$ described as $\stackrel{}{{{\mathit{v}}_{\mathit{i}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathit{i}}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{ngb}}}{\mathit{m}}_{\mathit{j}}\hspace{0.17em}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\hspace{0.17em}{\left|{{v}}_{\mathit{i}}\mathrm{-}{{v}}_{\mathit{j}}\right|}^{\mathrm{2}}\mathit{,}$(12)where v_i is the velocity of particle i and W(r,h) is the SPH kernel function. Different values of the diffusion parameter d are discussed in the context of isolated SN explosions (Sect. 7) and dSph modellings (Sect. 8).

6.2. Smooth metallicity

In particle based simulations, any elemental abundance can be computed as the ratio between the mass ${\mathit{m}}_{\mathit{i}}^{\mathrm{X}}$ of the element X and the total hydrogen mass ${\mathit{m}}_{\mathit{i}}^{\mathrm{H}}$ for each particle i, namely, $\begin{array}{ccc}\mathrm{[}\mathit{X}\mathit{/}\mathrm{H}{\mathrm{]}}_{\mathit{i}}\mathrm{=}{\log }_{\mathrm{10}}\left[\frac{\mathrm{(}{\mathit{m}}_{\mathit{i}}^{\mathit{X}}\mathit{/}{{\mathit{m}}_{\mathit{i}}^{\mathrm{H}}}^{\mathrm{)}}}{{\mathrm{(}{\mathit{m}}_{\mathit{i}}^{\mathit{X}}\mathit{/}{{\mathit{m}}_{\mathit{i}}^{\mathrm{H}}}^{\mathrm{)}}}_{\mathrm{\odot }}}\right]\mathit{,}& & \end{array}$(13)where $\mathrm{(}{\mathit{m}}_{\mathit{i}}^{\mathit{X}}\mathit{/}{\mathit{m}}_{\mathit{i}}^{\mathrm{H}}{\mathrm{)}}_{\mathrm{\odot }}$ is the solar abundance ratio of an element X.

With an SPH scheme however, extensive physical quantities linked to one particle are defined through a convolution running over the nearest neighbours. Following Wiersma et al. (2009)³, Eq. (13) may be reformulated with the density ${\mathit{\rho }}_{\mathit{i}}^{\mathit{X}}$ and ${\mathit{\rho }}_{\mathit{i}}^{\mathrm{H}}$ to become $\begin{array}{ccc}\mathrm{[}\mathit{X}\mathit{/}\mathrm{H}{\mathrm{]}}_{\mathit{i}}^{\mathit{s}}\mathrm{=}{\log }_{\mathrm{10}}\left[\frac{\mathrm{(}{\mathit{\rho }}_{\mathit{i}}^{\mathit{X}}\mathit{/}{{\mathit{\rho }}_{\mathit{i}}^{\mathrm{H}}}^{\mathrm{)}}}{{\mathrm{(}{\mathit{m}}_{\mathit{i}}^{\mathit{X}}\mathit{/}{{\mathit{m}}_{\mathit{i}}^{\mathrm{H}}}^{\mathrm{)}}}_{\mathrm{\odot }}}\right]\mathit{,}& & \end{array}$(14)with ${\mathit{\rho }}_{\mathit{i}}^{\mathit{X}}$ defined by $\begin{array}{ccc}{\mathit{\rho }}_{\mathit{i}}^{\mathit{X}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{ngb}}}{\mathit{\rho }}_{\mathit{j}}^{\mathit{X}}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}}\hspace{0.17em}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{ngb}}}{\mathit{m}}_{\mathit{j}}^{\mathit{X}}\hspace{0.17em}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathit{.}& & \end{array}$(15)This ratio $\mathrm{[}\mathit{X}\mathit{/}\mathrm{H}{\mathrm{]}}_{\mathit{i}}^{\mathrm{s}}$ is referred to as the smooth metallicity because it results from the weighted contribution of neighbouring particles. This scheme smooths the metallicity gradient between close particles, as would be expected by a turbulent ISM. Owing to its definition, this method naturally takes the resolution of the simulation into account.

6.3. Fundamental differences between the two schemes

Contrary to the diffusion scheme, the smooth metallicity scheme does not explicitly redistribute metals among particles. Its effect is local, extending the enriched region according to the size of the SPH smoothing length, instantaneous, as it is performed at every timestep; and static, as the enriched region does not change if particles are not moving sensitively. In contrast, the diffusion redistributes metals among particles and continuously extends the enriched region with time as long as the velocity dispersion is non-zero, which is always the case in practice. In this sense, it is a dynamical process.

7. Isolated supernovae explosions

Simple tests provide the first steps to understanding the effect of feedback and elemental dispersion inside numerical simulations. Towards this purpose we perform a series of simulations of isolated supernovae explosions inside a homogeonous gaseous medium.

We look in particular at the impact of the resolution, the adaptive timesteps, the artificial viscosity, and the method of the metal deposition, including the smooth metallicity scheme and/or the metal diffusion. We also take care to ensure the conservation of energy.

7.1. Initial conditions

In these scenarios a single supernova is detonated within a periodic box of side 1 kpc filled with a pristine gas of uniform density. The initial density is set to 0.1 atom / cm³ and the temperature to T = 10⁴ K, which corresponds roughly to the ISM of the first supernovae explosions in our dSph models (Revaz & Jablonka 2012). The metallicity is set to a primordial abundance of (by mass) 76% hydrogen and 24% helium. In order to avoid undesired noise, the gas particles in the box are first relaxed to form a glass-like structure.

To examine the impact of the resolution, these simulations were run with five different number of particles between N_gas = 16³ and 256³ (from 600 M_⊙ down to 0.15 M_⊙, in terms of mass resolution). The smallest resolution is of order the dSph simulations of Revaz & Jablonka (2012) and the highest is for convergence tests alone.

For each resolution, we explored the impact of the artificial viscosity improvement (Sect. 3.3), the adaptive timesteps improvement (Sect. 3.2), and the diffusion or smooth metallicity schemes (Sect. 6). The parameters used for each model are summarised in Table 1.

7.2. Sedov-Taylor solution

We confirm that in our highest resolution simulation (top panel of Fig. 5) the Monaghan (1997) artificial viscosity (see Sect. 3.3) formulation in combination with the adaptive timestep improvement proposed by Durier & Dalla Vecchia (2012, see Sect. 3.2), reproduces the Sedov-Taylor solution. During the supernovas expansion, a non-negligible amount of particles remain at low densities inside the expanding bubble. The shell of this bubble is sufficiently thin and dense, agreeing with the analytic density profile. This is also the case for other variables such as the temperature and radial velocity. The total energy is conserved at a 0.3% level.

At the lowest resolution, where gas particles have a mass of 500 M_⊙ or higher, typically used for galaxy formation simulations, the blastwave is not well recovered. As the SPH kernel has approximately a fixed number of particles (Eq. (10)), the total mass inside the kernel rises at low resolution. This naturally results in the dilution of the feedback energy per unit mass and consequently the increase of internal energy/temperature per particle is weakened (Dalla Vecchia & Schaye 2012). The impact of this can be up to two dex between the N = 16³ and the N = 256³ models. However the transfer of the internal energy to kinetic energy only depends weakly on the resolution (at a 2% level). As a consequence, the net impact outside the kernel remains similar, independent of the resolution.

When radiative cooling is considered in these Sedov-Taylor experiments, the resolution dependence of temperature and density becomes important since the cooling function is directly related to these quantities. At low resolution, the central particles do not reach the low density and high temperature regime where the cooling is significantly reduced (Dalla Vecchia & Schaye 2012). Following Stinson et al. (2006, and as in Revaz & Jablonka 2012 we use an adiabatic time t_ad = 5 Myr, during which the cooling of any particle that received energy from the supernova is switched off. This fix avoids an instantaneous cooling of the hot gas and improves the convergence.

7.2.1. Adaptive timesteps

Without the improved timestep scheme (model “c”) proposed by Durier & Dalla Vecchia (2012), gas particles in the ISM surrounding the explosion fail to be activated early on and arbitrarily large errors can then occur.

7.2.2. Artificial viscosity

Figure 5 compares the two artificial viscosity schemes. The top panel corresponds to the viscosity scheme adopted by GEAR (model “a”) while the bottom panel represents the standard Gadget-2 implementation (model “b”). At high resolutions and using the standard Gadget-2 artificial viscosity, some gas particles can penetrate into the ISM through the supernova shell (see the gas surface density) whereas this is not the case in the Monaghan & Gingold (1983) formulation.

The difference between the two methods is explained by the presence of the additional factor h_ij/r_ij in the definition of μ_ij compared with the standard Gadget-2 viscosity (see Eq. (2)). This factor increases as particles come closer and boosts the effect of the bulk viscosity⁴. At low resolution this difference is negligible, however, as the shock is not strong enough for this to become important. Finally, we stress that the total energy conservation remains very similar between the two schemes.

7.3. Metal distribution and the effect of smooth metallicity or diffusion

7.3.1. Radial distribution of ejecta

The [Fe/H] ratio of gas particles as a function of time is shown in Fig. 6 for model “a” with N = 16³ and N = 128³ particles. Without the smooth metallicity scheme or a diffusion term, only the particles that are directly touched by feedback contain metals. Furthermore, the [Fe/H] profile directly traces the kernel at the moment of the supernova explosion because the metal is distributed according to the SPH kernel, and higher resolution models distribute metals over a smaller physical volume. This reduced volume (and consequently reduced mass of gas that receives ejecta) means that higher resolution models end up with much higher metallicities.

The width of the distribution however remains unchanged. Owing to the fact that the energy transfer into a kinetic form is resolution independent, the velocity of particles is on average smaller in the low resolutions and the growth of the polluted region is subsequently three to four times slower in low-resolution cases. This reduction in speed helps to reduce the size discrepancy of the final polluted region between different resolution models.

7.3.2. Effects of smooth metallicity and metal diffusion schemes

In contrast to the standard enrichment schemes, both the smooth metallicity scheme (model “d”) and metal diffusion scheme (models “e” to “h”) allows particles initially outside the feedback region to contain metals. Figure 7 shows the metallicity 5 Myr after a supernova explosion for the standard SPH scheme, the smooth metallicity scheme, and diffusions schemes (covering a diffusion coefficient d of 0.003, 0.001, 0.0003 and, 0.0001).

Both methods reduce the strong metallicity gradient that exists in the standard SPH scheme. Particles outside the region initially touched by feedback are now able to contain metals produced from the supernova, and the level dependent is on the model (and diffusion coefficient) used. Consequently, the metallicity within the inner region is reduced in both schemes either because metals diffused outwards (diffusion scheme) or the presence of low-metallicity neighbours (smooth metallicity) in the kernel of the central particles.

The major differences between the two methods occurs because of their fundamental nature, i.e. whether they are static and instantaneous or a dynamical process. At low resolutions (N = 16³, h ≅ 0.15 kpc), the smooth metallicity scheme affects a region approximately twice the size of the region initially enriched; this region is much larger than that affected by diffusion regardless of the diffusion coefficient. At high resolutions (N = 128³, h ≅ 0.035 kpc), the smooth metallicity scheme only slightly broadens the initial polluted region. Meanwhile, the diffusion scheme remains weakly affected by the resolution and at high resolutions becomes the most efficient mechanism of metal mixing.

In Fig. 7, the metallicity distribution function (MDF) for each model considered is shown in the upper right corner. Independent of the resolution, the difference of the two methods is clear. With the diffusion scheme, a lot of particles receive a tiny fraction of elements, thus dominating the low-metallicity region of the MDF. Hence, the peak of the MDF is washed out. On the contrary, the smooth metallicity scheme only slightly shifts the MDF towards lower metallicity while keeping the peak of the distribution at an intermediate metallicity.

We finally point out that the effect of both methods is, by construction, resolution dependent, and the mixing is always strongest at the lowest resolution.

7.4. Two supernovae and abundance ratio scatter

In the previous experiments where only a single supernova exploded in a box of pristine gas, the final abundance ratios of all gaseous particles were strictly identical since they all received the same ratio of elements. Studying the scatter in elemental abundances requires simulations with multiple supernovae of different yields.

Here, we study the impact of various processes on the scatter in a dual supernovae scenario. First, an α-rich ([ Mg / Fe ] ≅ 1) 40 M_⊙ supernova is exploded and then 5 Myr later, a second, α-depleted ([ Mg / Fe ] ≅ −0.35) 15 M_⊙ supernova explodes, ensuring a large scatter in abundances ratio.

We performed 24 simulations, each a combination of (i) a mixing scheme: the smooth metallicity, a strong diffusion (d = 0.003), a weak diffusion (d = 0.0001), or no mixing; and (ii) an ejection scheme: the normal SPH kernel, no kernel dependence (w_ij = Γ_ij), a fixed radius (R_ej = 0.125 kpc), or a blast radius (R_ej = R_E), as described in Sect. 5. For each of these simulations, we explored two resolutions: N = 16³ and N = 64³. In the fixed radius scheme, the value R_ej = 0.125 kpc corresponds to the blast radius R_E obtained for a medium with a density of 0.1 atom / cm³ and a temperature of 10⁴ K. As we found only a small difference between the fixed radius and blast radius, we only show results for the former to avoid clutter.

Figure 8 shows the [Fe/H] of the gas as a function of the distance to the centre of the explosion, and [Mg/Fe] of the gas as a function of [Fe/H] for a subset of the most relevant models with a resolution of N = 64³. At the top right of each panel, we also show the final metallicity distribution. For each simulation, the abundance of the gas at four different times is indicated: t₁ = 0.6 Myr (black), just after the first explosion; t₂ = 5.2 Myr (blue), just before the second explosion; t₃ = 5.3 Myr (green), just after the second explosion; t₄ = 9.4 Myr (red), at the end of the run.

7.4.1. Normal SPH kernel, no mixing

With the normal SPH kernel and without mixing, 49 particles are enriched by the first explosion, which is in agreement with the choice of 50 neighbouring particles. As a result of the absence of mixing, the metallicity distribution traces the SPH kernel with the central regions more heavily enriched than the other regions. As each metal species are distributed identically, the touched particles have a constant [ α/ Fe ] value, forming the blue horizontal line. Between the two explosions, the particles enriched by the first supernovae are overpressurised and travel outwards forming the supernovae blast wave. Interestingly, the second explosion impacts 84 particles. As evoked in Sect. 5.2, such a large deviation compared to the expected ± 50 particles is a direct consequence of the inhomogeneous medium generated by the blastwave. In addition to the 49 particles enriched by the first supernova, 35 particles comprised of pristine gas receive metals and these all share the same α-depleted yields, that of the 15 M_⊙ supernova (forming the red horizontal line at the bottom of the [Mg/Fe] vs. [Fe/H] plot). In the central regions of the blast wave, the initially enriched particles are mixed with α-depleted materials and end up with a [ Mg / Fe ] of around 0.5. The resulting scatter is large due to the presence of extreme α-rich and α-depleted regions. As the number of particles touched by the second supernova increases, the ejecta are also more diluted.

7.4.2. Step function kernel (w_ij = Γ_ij), no mixing

In this situation, all particles touched by one explosion receive the same amount of ejecta. The radial [Fe/H] distribution is thus perfectly flat until the edge of the kernel. The second explosion simply shifts the metallicity upwards and all particles are grouped in a small cluster on the [Mg/Fe] vs. [Fe/H] plot. The scatter in abundance ratio is thus extremely low. As opposed to the previous case, the particles touched by the feedback receive equal amounts of energy and thus travel the same length and maintain a similar density. Consequently, as the density is kept constant and homogeneous at the time of the second explosion, only the same 49 particles are touched by the ejecta of the second supernova.

7.4.3. Fixed radius (r = 0.125 kpc), no mixing

When the ejection radius is fixed, the bias related to the adaptive SPH radius is staved off. The second supernova only concerns the particles that are still lying in 0.125 kpc and, thus, because of the kernel, are the most enriched. This double enrichment results in the radial profile bump traced by the green line. As a result of the α-depleted nature of the second supernova, these particles see their [Mg/Fe] decreased. Further away (where particles have low [ Fe / H ]), the α abundances remain unchanged. The low [Mg/Fe] floor present in the normal SPH simulation without mixing is completely absent.

7.4.4. Normal SPH kernel, smooth metallicity

The smoothed metallicity dilutes the ejecta across a larger number of particles, shifting the mean to a lower metallicity. Interestingly, this mixing scheme does not avoid the large scatter in [Mg/Fe] but only reduces it slightly by averaging some extreme values.

7.4.5. Normal SPH kernel, diffusion

The addition of a strong diffusion coefficient (d = 0.003) allows the transfer of metals from the inner metal-rich regions to the outer metal-poor regions. When the second supernova explodes, no pristine gas is affected by the α-depleted ejecta. Particles far from the second supernova, which would be pristine if not for diffusion, receive only a small fraction of ejecta relative to that received via diffusion and consequently maintain a [Mg/Fe] ratio around 1. On the contrary, in the central regions that have lost elements through diffusion, the impact of the second supernova is stronger and makes the [Mg/Fe] ratio gradually decrease down to around 0.

When the diffusion coefficient is decreased (d = 0.0001), the result is intermediate between the normal SPH case and the case with the strong diffusion. In this case the high-metallicity regions follow the standard case, while the low-metallicity are similar to the strong diffusion case.

8. Simulations of dwarf spheroidal galaxies

We simulate dwarf spheroidal galaxies (dSph) to study the effect of the numerical schemes on observable properties. Here, we systematically alter IMF modelisation (Sect. 4), how elements are distributed after supernovae explosions (Sect. 5), and mixing of metals via the smooth metallicity approach (Sect. 6.2) or diffusion (Sect. 6.1).

We focus on star formation rates as well as on stellar abundances, both of which are well constrained by recent observations (see Sect. 1). As dSphs form nearly all stars in situ with only minimal accretions, they are more relevant for a direct comparison than more massive galaxies. Furthermore, as dwarf galaxies are small systems, their modelling is less CPU intensive, which allows us to explore a large parameter space.

8.1. Initial conditions

The dSph models are based upon updated versions of the models described in (Revaz & Jablonka 2012). The initial conditions are based on a two-slope density spherical profile $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}}^{\mathrm{)}}}^{\mathit{a}}{\mathrm{(}\mathrm{1}\mathrm{+}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}}^{\mathrm{)}}}^{\mathit{b}\mathrm{-}\mathit{a}}}\mathit{,}$(16)where a is the slope of the inner profile, fixed to 0 to correspond to a core favoured by observations of normal low brightness and dwarf galaxies (Blais-Ouellette et al. 2001; de Blok & Bosma 2002; Swaters et al. 2003; Gentile et al. 2004, 2005; Spekkens et al. 2005; de Blok 2005; de Blok et al. 2008; Spano et al. 2008; Walker & Peñarrubia 2011; Oh et al. 2011), b is the slope of the outer profile, fixed to 3, in agreement with the outer NFW profile (Navarro et al. 1996, 1997) from ΛCDM simulations. r_s, the scale radius is set to 2 kpc. The models are truncated at a maximum radius r_max. Initially, the models only contain dark matter and pristine gas both sharing the same profile with a respective mass ratio of 0.15 in agreement with the standard cosmological baryonic fraction.

The initial equilibrium of the halo is obtained by calculating the velocity dispersions derived from the Jean’s equations, assuming a spherical symmetry. For the gas, these velocities are transformed in thermal energy by multiplying by 0.5 such that the gas sphere smoothly condenses when the simulation is started.

8.2. Dwarf spheroidal models

Two different sets of initial conditions have been considered, as much as possible both reproducing the observed properties of the Local Group dSphs Sextans and Fornax, which are two galaxies representative of the variety among dSphs. Fornax is expected to be rather massive owning to it high luminosity (15.5 × 10⁶L_⊙) and metal-rich stellar population (⟨ [ Fe / H ] ⟩ = −1.07). It is characterised by an extended star formation rate (Coleman & de Jong 2008). At the opposite end, Sextans is substantially less luminous (L = 0.53 × 10⁶L_⊙) and characterised by a metal-poor stellar population (⟨ [ Fe / H ] ⟩ = −2.02) dominated by old stars (Lee et al. 2003). The parameters used for each one of these two models are given in Table 2.

Each of these two models is run using different resolutions spanning nearly three orders of magnitude in mass. The resolution is defined by an integer r. The corresponding stellar mass m_⋆,r and spatial resolution (gravitational softening) ϵ_r, are computed using the relations ${\mathit{m}}_{\mathit{\star }\mathit{,}r}\mathrm{=}{\mathit{m}}_{\mathrm{0}}\mathit{/}{\mathrm{2}}^r \mathrm{and} {\mathit{ϵ}}_r\mathrm{=}{\mathit{ϵ}}_{\mathrm{0}}\mathrm{·}{\mathrm{2}}^{\mathrm{-}r\mathit{/}\mathrm{3}}\mathit{,}$(17)where m₀ = 65536 M_⊙ and ϵ₀ = 0.24 kpc correspond to the poorest resolution considered (r= 0). With this definition, incrementing the resolution by 3 corresponds to dividing the mass by a factor eight and the spatial resolution by two.

The complete set of simulations including their parameters is given in Table 3, representing more than 80 simulations. A parameter set to “SPH” means that the default SPH setting is used for both w_ij and R_ej. We also test additional parameters, such as the adiabatic time (Sect. 7.2) applied to either SNIa or SNII. We defined a fiducial set of parameters shown on the top of the table to ease the comparison between the models and increase readability. When the value of a parameter is not given, the one of the fiducial model is used. For all models, the star formation parameter c_⋆ is fixed to 0.01 and the supernova efficiency ϵ_SN is set to 0.3.

8.3. The effect of the IMF sampling

In Sect. 4 we presented different approaches to sample the IMF. Briefly, these can be summarised as: continuous IMF sampling (CIMFS, see Sect. 4.1), where each stellar mass contributes in ratio to its proportion of the IMF (including fractional contributions, if this is less than one); random sampling (RIMFS, Sect. 4.2), where the IMF is treated as a probability distribution and supernovae are generated stochastically; and “optimal” IMF sampling (OIMF, Sect. 4.3), where the IMF is split into mass bins such that each bin contains one supernova at the lower end of the mass bin. The choice of IMF sampling method has a direct impact on both the feedback energy and the elements injected as supernovae go off inside the galaxy. Our aim here is to address the impact of the IMF choice on the global chemical enrichment of a galaxy.

In Fig. 9, we delineate the star formation history and evolution of the stellar mass for the three different schemes as a function of the resolution from r= 4 (m_⋆ = 4096 M_⊙) to r= 6 (m_⋆ = 1024 M_⊙). We analyse and compare these simulations in detail below.

8.3.1. CIMFS and the SNIa adiabatic time

The first and foremost feature observed in the CIMFS method is the quenching of star formation after about 3 Gyr, independent of the resolution or the dSph model. Between 2 and 10 times fewer stars are generated in the case of the CIMFS which drastically impacts the final metallicity and abundances.

This quenching is a direct consequence of the heating of IMF by the SNIa and may be understood as follows. For the r= 4 resolution runs (m_⋆ = 4096 M_⊙), the typical timestep is 0.1 Myr. According to the Kroupa IMF (Kroupa 2001) and the Kobayashi et al. (2000) SNIa model, the number of SNIa exploding during this timestep is about 10^-4. This means that small fractions of SNIa continuously explode over a long time interval, forcing the adiabatic time switch to be always active. Consequently, the gas cooling is artificially quenched along with the star formation rate. The problem worsens with increasing resolution.

In contrast, only full SNe explode in the RIMFS and CIMFS approaches, respecting the physical delay in time between all of them (see Fig. 2). As t_ad = 5 Myr is small compared to the mean time between two SNe, the cooling is still effective in maintaining a continuous star formation rate.

In a fourth model (CIMFS + t_ad,SNIa = 0), we tested the suppression of the adiabatic period only after the SNIa explosions, since it is sometimes used in galaxy formation simulations. Setting the adiabatic time for the SNIa to zero in the CIMFS generates a star formation rate akin to the RIMFS and OIMFS. However, we still have the problem that SNIa dilute the ejection of metals in time affecting the chemical evolution of the system. As in this case, at every timestep all stellar particles experimenting SNe explosions expel at least a fraction of elements, and the mixing is artificially boosted. However, it is clear that CIMFS with t_ad,SNIa ≠ 0 does not reproduce any feature accurately and consequently, we will no longer consider those models in the future.

8.3.2. The effect of the resolution

The second clear feature on Fig. 9 is the important effect of the resolution on the star formation rate. As long as the resolution is low (r≤ 4, m_⋆ ≥ 4096 M_⊙) all three sampling methods (RIMFS, OIMFS and CIMFS + t_ad,SNIa = 0) nicely converge.

As resolution increases the cumulative stellar mass begins to diverge, at the r= 5 resolution (m_⋆ = 2048 M_⊙), the RIMFS cumulative stellar mass deviates of about 10% from the OIMFS and CIMFS + t_ad,SNIa = 0. For Sextans, with r= 6 (m_⋆ = 1024 M_⊙), the RIMFS produces 40% more stars than the CIMFS + t_ad,SNIa = 0 and 30% more that the OIMFS. That the OIMFS is closer to the CIMFS because it represents a better fit to the continuous solution, as seen in Fig 2.

The deviation of the three methods with increasing resolution is a direct consequence of the method chosen for IMF sampling, a major problem for high-resolution simulations. As we see in what follows, all of these methods bias the predictions of the chemical abundances once a given resolution is reached.

8.3.3. Impact on the chemical abundances

In addition to its effect on the global star formation rate, we also want to address the impact of the IMF sampling scheme on the final abundances.

A reliable comparison between the different methods is not straightforward, as they all generate slightly different star formation histories that alter the final abundances. We firstly compare the models at 1 Gyr when the amount of stars formed in the Sextans models are all similar, according to the bottom of Fig. 9. In Fig. 10, we show the magnesium abundance of the gas ([Mg/Fe]), our tracer of the α-elements, as a function of the metallicity of the gas ([Fe/H]). In Fig. 9, we indicate the three different IMF sampling methods and the three resolutions. In each panel, a point corresponds to one gaseous particle where the ratios [Mg/Fe] and [Fe/H] is computed with the smooth metallicity scheme (Eq. (14)). On top of that, the continuous and dashed black lines indicate the 1σ dispersion of each sample. The values of the 1σ dispersions are also reported on the bottom of the plot for a direct comparison between the models.

At low resolution (r= 4), all models display equivalent distributions with very low dispersions (1σ< 0.1). This arises because, at such low resolution, the neighbouring particles of a stellar particle cover a large volume, subsequently, ejecta are efficiently spread over the entire galaxy resulting in an excellent mixing. As expected, if the yields of all type II supernovae are mixed altogether, a mean value of [ Mg / Fe ] = 0.5 is obtained, forming a clear plateau.

With increasing resolution, the dispersion increases and the differences between RIMFS, OIMFS and CIMFS + t_ad,SNIa = 0 become clear. Firstly, it is clear that in all cases, the dispersion decreases with increasing [Fe/H]. This is the direct consequence of the internal mixing that each particle experiments when it receives metals from successive supernovae. Secondly, the plots clearly show that the RIMFS scheme always generates dispersions about a factor of two larger than both the OIMFS and the CIMFS + t_ad,SNIa = 0 models. The mixing is obviously stronger for the CIMFS + t_ad,SNIa = 0 because, in this particular scheme, at every timestep, each stellar particle that is older than the shortest star lifetime spreads ejecta. A neighbouring gaseous particle is therefore more likely to receive ejecta from the whole range of stellar masses considered in the IMF. This is in contrast with the RIMFS scheme where extreme cases, such as only a single 50 M_⊙ SNe exploding for a stellar particle, may occur. In such a case, neighbouring gaseous particles are biased towards very high [Mg/Fe] ratios. Similar cases such as low-mass supernovae exploding would bias the neighbouring particles towards low [Mg/Fe] ratios. The final dispersion is then larger for RIMFS than for the other schemes.

The OIMFS scheme, despite its discrete nature, aims to sample the IMF in a more continuous way, ensuring that a wide range of supernova masses are represented. The dispersion obtained is indeed very similar to the CIMFS + t_ad,SNIa = 0 dispersion. However, this method suffers from an important bias. As described in Sects. 4.3 and 4.4, for a given resolution (m_⋆,r = M_SSP), an upper local stellar mass m_max exists, and stars with higher masses are ignored. For the models considered here, r= 4,5,6, the corresponding maximal masses are respectively 40.6, 34.5, and 26.9 M_⊙. Hence, with increasing resolution the OIMFS scheme misses the ejecta of the most massive stars, which are those with a high α-element ratio. As clearly seen for the r= 6 model, the low-metallicity plateau is offset towards lower [Mg/Fe], with respect to the CIMFS + t_ad,SNIa = 0 and RIMF schemes.

As observations only view the stars as they are today, we show in Fig. 11 the [Mg/Fe] ratio for the Fornax models, at t = 14 Gyr. Despite the fact that the quantity of stars formed in these models is different, owing to different star formation rates, these plots allow us to appreciate the effect of the different IMF sampling schemes on the final abundances. For comparison, the observed dispersion in [Mg/Fe] extracted from all stars of Fig. 1 is shown as a yellow curve on the bottom of each plot. While the trends discussed in the gas abundances are less pronounced in the stars, they all remain true: (i) at low resolutions, the metallicity dispersion is low (even lower than that observed) and the three methods converge; (ii) at higher resolutions, the dispersion in the RIMFS scheme is always larger than the CIMFS + t_ad,SNIa = 0 and OIMFS schemes; and (iii) the [Mg/Fe] plateau is offset owing to the lack of massive supernovae in the OIMFS case.

8.3.4. Convergence

We briefly discuss here some convergence issues related to the resolution. Figure 12 shows the star formation and cumulative stellar mass for the Fornax model with resolutions increasing from r= 0 to 9, i.e, covering two dex in terms of mass resolution. Here, no convergence is obtained with the RIMFS scheme⁵, which primarily occurs for two reasons.

Firstly, for the high-resolution model (r= 9), the mass of one stellar particle is only 128 M_⊙ inevitably causing bias in the randomly generated IMF as described in Sect. 4.2. We found that the RIMFS simulation with a resolution lower that r= 6 generates precisely 87 × 10^-4 SNeII per solar mass formed, in agreement with the choice of the IMF. For the high-resolution model, this value can drop to 75 × 10^-4 owing to incomplete sampling. Consequently, the feedback associated with high-resolution simulations is smaller, fostering a higher star formation rate. The divergence between models is considerably smaller for the CIMFS + t_ad,SNIa = 0 case. There, according to the definition of the CIMFS, the number of SNeII per solar mass formed is always precisely the same, independent of the resolution. In addition, energy is released at every timestep, preventing the runaway star formation.

Secondly, because the SPH radius is smaller for higher resolutions, energy is deposited into a smaller volume of gas, preventing the (temporary) shutdown in star formation over large sections of the galaxy seen in low-resolution models.

8.4. Reproducing the scatter in [α/Fe]

A fundamental point when simulating the chemical evolution of a galaxy is to correctly reproduce the scatter in the element abundance ratios. In this section we study the impact on the dSph chemical evolution of the elements spreading scheme discussed in Sect. 5 along with the two mixing methods introduced in Sect. 6, namely the smooth metallicity (Wiersma et al. 2009) and diffusion process based on velocity dispersions (Greif et al. 2009).

8.4.1. The scatter without mixing

We first discuss our models in the absence of any additional mixing, by focusing on the final stellar [Mg/Fe] vs. [Fe/H] distributions for Fornax and Sextans. All three IMF sampling methods suffer from different issues at the resolution considered here (r= 6, m_⋆ = 1024 M_⊙). The CIMFS artificially boosts the mixing of elements, the OIMFS induces a spurious offset in the α-elements abundances ratio, and the RIMFS is subject to Poissonian noise. We however found the RIMFS scheme to be better suited to study the influence of different numerical techniques on the scatter of elements abundances ratio. Thus for the rest of the discussion, we only use the RIMFS scheme. The results are indicated in dark blue on panels (a) and (e) of Fig. 13, with observed values overplotted in yellow (from Fig. 1). We compare only the low-metallicity end of the stellar distribution ([ Fe / H ] < −2.5), where both Milky Way and dSph stars sit on the same high [ α/ Fe ] plateau, before they diverge because of differing star formation histories. This choice allows us to compare the scatter with simulations without running into complications arising from the small statistics in most dwarfs.

At low metallicity in the simulations, the scatter is above 0.4 dex, about 2–3 times higher than that observed in metal-poor stars. In particular, stars are found with extreme [Mg/Fe] values, above 0.85 and below −0.5. As seen by Fig. 1 very few metal-poor stars have been found with subsolar [Mg/Fe] and only two have values below −0.5 (Aoki et al. 2014; Jablonka et al. 2015). At high metallicity, the scatter decreases as each star has had a higher number of supernovae enriching it, however, many stars are still found with very low [ α/ Fe ] ratios. By examining the extreme stars we find that very low-metallictiy stars are enriched either purely from low-mass M< 20 M_⊙ supernovae or, in rare cases, by a single SNIa. Meanwhile α-rich stars all contain at least the ejecta of one SNII more massive than 25 M_⊙ and few if any below 20 M_⊙. At higher metallicity, the evolution is more complex but extremely low α abundances at higher metallicity ([ Fe / H ] > −2) only occur because of the presence of SNIa ejecta.

In all cases, stars with extreme α abundances result from the pollution of an incomplete IMF. This is supported by the mean number of SNII supernovae responsible for the abundances in the Sextans model stars: only 13 for low-metallicity stars and 48 for the more metal-rich but extremely α-depleted stars, while the mean values for all stars is about 70. The frequency of these stars in simulations and the excessive scatter point to limitations in the metal ejection scheme and possibly insufficient mixing in the standard SPH technique. We examine some proposed improvements of the method below.

8.4.2. Metal spreading schemes

Simple experiments involving two supernovae (Sect. 7.4) have shown how the scatter in element abundance ratios in the gas may be influenced by the method used to distribute metals throughout the nearby gas. Here, we apply the same recipes to our dSph models and check their impact on the final abundances of stars.

With the first scheme explored (the step function distribution Γ), all neighbouring particles receive the same amount of material, independent of their distance, avoiding any gradient related to the kernel (w_ij = Γ_ij, Fig. 8b). The effect on our Fornax simulation is shown in panel (b) of Fig. 13. In contrast to what would be expected from the double supernovae experiments, the scatter is not decreased. This lack of decrease suggests the scatter does not originate in the kernel, but by stochastic supernovae explosions having different masses and yields, as discussed in Sect. 8.4.1. At high metallicity ([ Fe / H ] > −2.5), where the stochastic effect is less important, the scatter is slightly lowered.

Panels (c) and (d) show the impact of fixing the ejection radius to either a constant value (R_ej = 0.125 kpc) or to the blast radius (R_ej = R_E). None of these two recipes is able to decrease the scatter, which is in agreement that the scatter does not arise because of the choice of the kernel. In the blast radius case, the effect is even more dramatic, and the scatter is increased at all metallicities. The constant ejected radius allows a better mixing compared to the blast radius. The resulting scatter is similar to the default case, except above [ Fe / H ] > −2.0, where it is slightly higher. In both cases, the mean value of the [Mg/Fe] low-metallicity plateau is too high, lying well above the mean observed value of around 0.5.

One can understand this through an examination of the supernova test run with a fixed radius (r = 0.125 kpc) and no mixing (Fig. 8). Using a fixed ejecta radius, the second supernova is only able to affect particles that still lie in R_ej. The most distant particles touched by the first massive supernova have moved outwards and are only polluted by α-enhanced ejecta. Consequently, these particles form an α-enhanced plateau at low metallicity.

8.4.3. Smooth metallicity

The smooth metallicity technique is a simple and natural way of mixing elements inside the SPH smoothing kernel. Its effect, in the context of dSphs, is shown in Fig. 13a and e for the Fornax and Sextans models, respectively, where the [Mg/Fe] vs. [Fe/H] distribution of the stars (in light blue) is compared to the case where no mixing is applied (in dark blue). In these simulations, the star formation histories are similar, and the final stellar masses are within 1% in the Fornax simulations and to within 3% in the Sextans simulations. The reason for this small impact is because only the local cooling function of the gas is impacted and its impact on the global dynamics and star formation history is minimal. This similar star formation history allows us to reliably compare the simulations. As seen in Fig. 13 the smooth metallicity scheme, strongly decreases the scatter by a factor of 2 to 3, becoming consistent with the scatter observed in stars (shown in yellow). In particular, most of the stars with extreme abundances are forced towards the mean, reaching realistic values.

The final metallicity distribution function (MDF) is only slightly modified by the smooth metallicity scheme. With the most metal-rich particles losing metals and the metal poorest gaining, and this results in a metal-rich cut-off in the MDF, an effect more pronounced in Sextans because of the small number of stars formed.

8.4.4. Metal diffusion

Panels (f), (g), and (h) of Fig. 13 show the impact on the stellar [Mg/Fe] abundances of the metal diffusion for diffusion coefficients of d = 0.003, d = 0.001, and d = 0.0003 compared to the smooth metallicity scheme represented in light blue. We also run a simulation with d = 0.0001, however in this case the effect of the diffusion is very low compared to the fiducial case and, hence, is not discussed here.

It is clear that the diffusion plays a similar role to the smooth metallicity, and for the choice of 50 neighbours a nice agreement is obtained for a diffusion coefficient of d = 0.001. Higher values of the diffusion coefficient are able to reduce the scatter even more, while lower values slightly increase the scatter. This similarity is also confirmed by the MDF shown at the bottom of each panel. A similar result is found in the case of the Sextans simulations (Fig. 19).

8.4.5. Discussion

As both the smooth metallicity and diffusion produce similar results it is difficult to have a clear preference for one method over the other. The disadvantage of the diffusion approach is the need for an additional parameter. Its value may be approximately derived with physical arguments, but at the galactic scale numerical errors and resolution make this computation obsolete. This requires the value of this coefficient to be calibrated as presented here. On the other hand, the smooth metallicity approach appears more natural for an SPH technique, where the information stops at the scale of the SPH smoothing radius. It is important to point out that both methods are sensitive to the resolution. This is illustrated in Fig. 14; when the number of particles is multiplied by 4, the scatter increases by about 0.05 dex at low metallicity. Finally, it is also important to emphasise that the CPU time required for the two methods is identical and may not be used as an argument that one has an advantage over the other.

9. Simulations of spiral galaxies

A mixing scheme should work across multiple resolutions and galaxy sizes. As such, we test the effect of the smooth metallicity and metal diffusion in much more massive Milky-Way like spiral galaxies.

9.1. Initial conditions

The initial conditions of our models are inspired from the works of Bullock et al. (2001) and Kaufmann et al. (2006), which consist of an isolated halo containing a rotationally supported gas component that gently flows towards the centre and forms a rotating disk. Initially, both the gas and the dark matter follow an NFW density profile with a concentration parameter c equal to 8. The baryonic fraction is set to 10%. The total angular momentum of the system is given by λ = 0.12 according to the notation of Bullock et al. (2001). Its radial distribution is such that the halo initially starts with a solid body rotation around the z-axis. Both dark matter and gas are truncated at the virial radius. The total mass is M = 9.14 × 10¹¹M_⊙. We run the models using a single resolution where the stellar mass is 10⁵M_⊙.

We have explored three models. The first model is run without any kind of mixing, the second uses the smooth metallicity scheme and the last model uses the metal diffusion scheme with a diffusion parameter of d = 10^-3, corresponding to the optimal value according to the dSph simulations (Sect. 8). All of the models use the RIMFS scheme, which is, equivalent to the OIMFS and the CIMFS schemes at this resolution (see Sect. 8.3). The star formation parameter c_⋆ is fixed to 0.7, the supernova efficiency ϵ_SN is 0.08 and, as in the dSph simulations, an adiabatic time of 5 Myr is used.

9.2. Results

All three models (no mixing scheme, dispersion scheme, and smooth metallicity scheme) present very similar morphological and physical features.

Hence, in Fig. 15, we only show the smooth metallicity scheme to illustrate these properties, where the galaxy is shown after 14 Gyr of evolution. The first panel (a) shows the extended gaseous disk together with its spiral structure. Its corresponding temperature, computed in the plane z = 0 is given in panel (b). The green regions indicate a dominant quasi-isothermal gas component around 10⁴ K. The red hot spots (T> 10⁵ K) trace the regions recently affected by supernovae explosions. The metallicity map for the gas (c) exhibits a radial gradient similar to that observed in spirals. The [Fe/H] ratio decreases from about 1 at the centre, down to –4 at about 50 kpc. Panel (d) shows the exponential stellar disk with a bar clearly present at the centre. The last panel is an edge-on view of the gaseous disk (top panel of (e)) as well as the stellar component (bottom panel of (e)). Superheated low density bubbles are visible within the gaseous disk, a consequence of supernova explosions. In addition to these features, the model presents realistic dynamical properties given in Fig. 16. The total rotation curve shown in black is nearly constant around 200 km s^-1. The mean azimuthal gas velocity (in green) slightly decreases in the central region owing to the larger velocity dispersion (asymmetric drift effect). The gas velocity dispersion is nearly constant between 10 and 20 km s^-1. Figure 17 shows the star formation rate of the three models, all of which are similar. Furthermore, the star formation history is in good agreement with that recently proposed by Snaith et al. (2014), which reproduces the chemical abundances of long-lived stars of the Milky Way; the majority of the stars are formed between 0 and 5 Gyr. The resulting α-abundances for the three models are shown in Fig. 18. The first three panels show the [Mg/Fe] vs. [Fe/H] distribution of the stars at t = 14 Gyr. We selected the stellar particles located 6 to 10 kpc from the centre, in the disk plane, to compare with observations that take place in the solar neighbourhood. We compare the distribution with observational values (shown in yellow) obtained from high-resolution spectroscopy of Milky Way stars only, as in Fig. 1. Overplotted are lines corresponding to the mean [Mg/Fe] value, and the 1σ deviation. As can be clearly seen, the three models result in different levels of scatter. The fourth panel directly compares the scatter between the models and observations.

The base model (run without smooth metallicity or diffusion) as well as the diffusion model are clearly inconsistent with the observations. At very low metallicity [ Fe / H ] < −3, most of the stars are found with high [Mg/Fe]. In the range −3 < [ Fe / H ] < −1, a low [Mg/Fe] tail is present with stars having [Mg/Fe] as low as −1.5. These two features which are not observed are the direct result of a lack of gas mixing. Gas particles initially located in these external regions are not affected by SNe explosions. Eventually, such pristine particles may fall towards the galactic disk and be impacted as the α-rich ejecta of a massive supernova (M> 30 M_⊙). As the star formation rate is high during the first few Gyr, the polluted gas particles may quickly be transformed into stars. The opposite process can happen at later time (and thus slightly higher [Fe/H]), when lower mass supernovae explode, impacting the gas with low [Mg/Fe]. At low metallicity, the resulting scatter is up to four times the observed scatter of about 0.2 dex.

When smoothing the abundance on a scale corresponding to the SPH resolution scale (the smooth metallicity scheme), all these features are washed out and the scatter becomes very similar to the observed scatter. The model still tends to be slightly too metal rich, mainly due to an over-active stellar formation in the last giga-years.

The diffusion scheme, while lowering the scatter compared to the base model, is not sufficient to smooth out these features and results in an exaggerated scatter.

10. Conclusions

Using self-consistent, N-body, chemo-dynamical models of dSphs and Milky Way-like galaxies, we have investigated the impact of different numerical schemes on the star formation history of galaxies and the chemical properties of stars therein; these schemes are the IMF sampling schemes (Sect. 4), elemental spreading schemes (Sect. 5), and metal mixing schemes (Sect. 6). We focused in particular on the scatter in [α/Fe] observed in metal-poor stars both in dSphs and Milky Way halo.

As feedback processes (both element and energy injection) are at the heart of this work, we ensured the reliability of our code GEAR in a simple case, where one supernova explodes in an homogeneous medium. We nicely reproduce the Sedov-Taylor solution at very high resolution when the individual and adaptive timestep scheme improvements of Durier & Dalla Vecchia (2012) are used. Here we found that the artificial viscosity formulation of Monaghan & Gingold (1983) produces superior results to that of Monaghan (1997) (the default Gadget-2 viscosity formulation). Additional simple experiments involving two supernovae led us to understand the effect of different metal spreading schemes, that are essential to the understanding of pollution in complex galactic systems.

10.1. The IMF sampling schemes

We have tested the influence of different IMF sampling methods: a continuous scheme (CIMFS), an “optimal scheme" (OIMFS) and a random scheme (RIMFS).

The most important result is that all schemes impose, at different levels, a limit to the stellar mass resolution of the simulations. From theoretical considerations, we have found that:

• As the CIMFS continuously ejects metals and energy, below a stellar mass resolution of about 10⁵M_⊙ the products of SNeIa are diluted over unrealistic long periods of time, generating an artificial mixing. When coupled with the adiabatic time method (Stinson et al. 2006), the cooling is strongly underestimated, resulting in a severe quenching of star formation. Removing the adiabatic time within the CIMFS scheme prevents this anomalous quenching, but the dilution problem remains.

• The OIMFS is designed to reduce the noise in the IMF compared to the random sampling. Unfortunately, it does not allow the use of a low particle mass (and subsequently high resolution). Below about 10⁴M_⊙ the maximal stellar mass is limited to about 30 M_⊙. Ignoring the yields of the most massive stars induces a bias that results in an offset in the α-element abundance plateau at low metallicity.

• Because of its stochastic nature, the RIMFS induces noise in the IMF when the stellar mass resolution reaches below about 10⁴M_⊙. Assuming that the IMF is distributed among neighbouring particles, the stellar mass resolution may be decreased down to about 10³M_⊙. The IMF may no longer be considered complete for lower masses, with pockets of enriched gas showing unrealistically high scatter.

In practice, based on simulations of Fornax and Sextans dSph, we have found that for a low resolution (where the mass of stellar particles is greater or equal to 4096 M_⊙), the number of stars inside a single stellar particle is large enough to allow for a good sampling of the IMF regardless of the scheme used. This occurs with all three schemes, which results in similar star formations and abundance patterns. However, the different schemes severely diverge in our tests performed at even higher resolution. At a resolution of 1024 M_⊙, the RIMFS is the most reliable sampling method. At higher resolution (128 M_⊙), the RIMFS shows a lack of convergence owing to incomplete IMF sampling. This incomplete sampling fails to inject enough energy across the galaxy to counterbalance the cooling and effectively regulate star formation.

Consequently, at very high resolution (<1024 M_⊙), no method is reliable. This is a severe limitation for the forthcoming simulations of galaxies that approach or exceed this resolution. Indeed, it is now well established that supernova feedback plays a major role during the process of galaxy formation. However, the precise supernova rate is directly dependent on the IMF sampling. This is also true for gas cooling which directly depends on its metal enrichment. Accurately reproducing the feedback and consequently star formation and chemical enrichment at these low masses will require the development of novel methods to sample the IMF.

10.2. The scatter in abundance ratios

Focusing on the RIMFS scheme at the maximal reliable resolution of 1024 M_⊙, we have demonstrated that the scatter in stellar abundance ratios is about three times larger than that observed in the metal-poor stars of dSph galaxies. This scatter is due to the stochastic explosions of supernovae imprinting small pockets of gas with markedly different and sometimes extreme abundances. In a first attempt to reduce this scatter we have tested the impact of the scheme used to spread metals into the ISM. We found that

• no statistically distinguishable differences exist when combining a classical SPH kernel with either a volume or particle mass-weighting;

• replacing the classical SPH kernel with a step function (gas particles receive ejecta independent of their distance to the exploding supernova) does not reduce the scatter, which is still dominated by stochastic explosions of supernovae; and

• using a constant radius (either fixed or equal to the blast radius) greatly worsens the problem.

In a second attempt, we tested the impact of introducing a mixing scheme. We implemented both the smooth metallicity technique (Wiersma et al. 2009) and metal diffusion (Greif et al. 2009). Both are able to reduce the scatter to a realistic value inside dwarf spheroidal galaxies. The smooth metallicity scheme has the advantage that it is independent of any additional parameters and very natural with regard to the SPH method. Metal diffusion needs the introduction of a diffusion coefficient, a free parameter set that reproduces the observations (d = 10^-3).

As a final test, we simulated a Milky Way-like galaxy that agrees with the star formation history deduced by Snaith et al. (2014). The abundance of α-elements is only reproduced when using the smooth metallicity technique. The diffusion scheme with the same parameters calibrated in the context of dSphs is not able to adequately reduce the number of stars with extreme abundances.

We have shown that the current best practices in chemo-dynamical simulations delivers reliable prediction concerning the chemical properties of galaxies in restricted contexts only, where the stellar mass resolution is limited to value below (masses above) 10³M_⊙. We found the smooth metallicity scheme, combined with a random initial mass sampling scheme (RIMFS) to be the best combination to reproduce the dispersion of abundances. Increasing the mass resolution in future N-body chemo-dynamical modelling while avoiding any bias will require a drastic redesign of the classical star formation recipes presently used.

---

Acknowledgments

We thank the anonymous referee for his/her useful remarks which strongly improved the paper. We also thank Fabrice Durier for providing his adaptive timestep routine and Pierre North for precious suggestions and comments. We are grateful to the International Space Science Institute (ISSI) at Bern for their funding of the team “The first stars in dwarf galaxies”, as well as SCITAS, the scientific IT and application support of the Swiss Federal Institute of Technology in Lausanne (EPFL). As for the simulations, data reduction and galaxy maps have been performed using the parallelised Python pNbody package (http://lastro.epfl.ch/projects/pNbody/). This work was supported by the Swiss National Science Foundation.

References

All Tables

All Figures