format: text

text:
# Supernova Enrichment of Dwarf Spheroidal Galaxies

## 1 Introduction

In cosmological models dominated by cold dark matter (CDM), the amplitude of fluctuations increases toward shorter wavelengths. Small galaxies are therefore likely to be the first to form, and large galaxies subsequently form through the merger of smaller systems (e.g. White & Rees, 1978; Blumenthal et al., 1984; Cole et al., 1994; Klypin, Nolthenius, & Primack, 1997; Navarro, Frenk & White, 1997). Small systems observed today, such as dwarf spheroidal galaxies (dSph’s), may thus represent surviving “building blocks” of larger galaxies, and therefore much attention has been paid to understanding the gas evolution and first generation of star formation within them (e.g. Lin & Murray, 1994; Navarro & White, 1994; Mori et al., 1997; Weil, Eke & Efstathiou, 1998; Klypin et al., 1999; Sommer-Larsen, Gelato & Vedel, 1999; Burkert, 2000; Mori & Burkert, 2000; Dong, Murray, & Lin, 2003).

Their small sizes might indicate that dSph’s can also serve as relatively simple laboratories to study the effects of various physical processes upon star formation. In spite of their dynamical simplicity, however, they are found to have a wide range of star-formation and enrichment histories, with few clear trends (Mateo, 1998; Grebel, 2001; Harbeck et al., 2001). The data therefore indicate that multiple factors have affected star formation in dSph’s, making them important and interesting objects for the study of early star formation in galaxies. Factors that may have influenced star formation in dSph’s include their exposure to external ultraviolet radiation, the depth of their gravitational potential, their overall mass, their interactions with other dwarf systems, and their proximity to neighboring massive systems. The work described here is part of a series of papers in which we examine these various physical factors, with the goal of determining their relative importance upon the evolution of the dwarf galaxies.

An important observational fact concerning dSph’s is that almost all have multiple generations of stars, with many showing clear evidence of having undergone bursts of star formation, generally with one dominant, and multiple weaker events (Mateo, 1998; Grebel, 2001). In some systems, multiple generations of star formation may also be inferred from the observed spread in the metallicity of individual stars within them (Lehnert et al., 1992; Suntzeff et al., 1993; Côté, Oke, & Cohen, 2000; Shetrone, Côté, & Sargent, 2001). Multiple generations of stars are also inferred from Hubble Space Telescope observations of stellar populations in dwarf systems that are resolved both spatially as well as temporally (Gallagher et al., 1998; Dohm-Palmer et al., 1998a, b, 2002). Most systems with multiple generations of stars, however, display less self-enrichment than would be expected from their star-formation histories, and the total mass of heavy elements contained within these galaxies is far less than that expected to be produced by the massive stars that would have formed along with the currently observed low-mass stellar population (Smecker-Hane et al., 1996; Grebel, 1997).

Star formation within dSph’s may be strongly inhibited at early times by the ultraviolet background radiation (Efstathiou, 1992; Quinn, Katz, & Efstathiou, 1996; Kepner, Babul, & Spergel, 1997; Navarro & Steinmetz, 1997; Weinberg, Hernquist & Katz, 1997; Barkana & Loeb, 1999; Benson et al., 2002). External triggers may still, however, lead to cooling and star formation (Lin & Murray, 1994; Dong, Murray, & Lin, 2003). Following an initial burst of star formation, the remaining gas within the galaxy may be heavily affected by photionization, stellar winds, and supernova explosions. Such stellar feedback has been studied extensively as a means of driving the evolution of stellar systems. Photoionization feedback has been proposed as a means by which star formation may be self-regulated in dwarf systems (Babul & Rees, 1992; Andersen & Burkert, 2000; Dong, Murray, & Lin, 2003). In systems the size of globular clusters, photoionization-driven winds may lead to extensive mass loss, sufficient even to completely terminate star formation (Peebles & Dicke, 1968; Tenorio-Tagle et al., 1986). In small galactic systems as extended as dSph’s, however, it is ineffective at driving mass loss (Noriega-Crespo et al., 1989).

Energy input from the multiple supernova events following a burst of star formation may also be a means of driving mass loss. Substantial mass loss has been found in simulations of low-mass disk systems subjected to multiple supernovae (Mac Low & Ferrara 1999). Mass loss is expected to be less efficient, however, from more rounded systems, in which heated gas must travel a greater distance through the ambient medium before it is able to leave the system. Studies of the evolution of supernovae in high redshift dwarf systems find that the high densities and resultant high cooling efficiencies and external pressures make it relatively easy for dwarf galaxies to retain mass following supernova explosions. In order for them to experience significant mass loss, high ($`10`$%) star formation efficiencies are required (Dekel & Silk, 1986; Mac Low & McCray, 1988; Madau, Ferrara, & Rees, 2001; Mori, Ferrara, & Madau, 2002). Even at later times, high star formation efficiencies may be required in order for supernovae to drive extensive mass loss, due to radiative losses (Larson, 1974). We note, however, that mass loss through the Lagrange points may also lead to significant rates of mass loss, if the gas within the galaxy is ionized either by external or internal sources (Murray, Dong & Lin, 2003). Extensive star formation is, therefore, not required in order to account for the total removal of gas from dSph’s, as observed today.

The presence of multiple generations of stars within dSph’s implies that their early star formation efficiencies were not so high as to lead to extensive mass loss. The poor self-enrichment of later generations, however, suggests that enriched material from supernovae was preferentially lost, leaving behind gas that was at most slightly enriched. It was from this poorly enriched gas that subsequent generations of stars then formed. It is therefore of interest to examine the evolution of enriched material within dwarf galaxies that have undergone relatively inefficient bursts of star formation. In this paper, we examine the loss of enriched material from dSph’s which are subjected to a number of supernova events. We proceed in § 2 by describing the numerical method and the setup of the models. The results of the simulations are described in § 3 and § 4. The possible consequences for the evolution of dwarf galaxies are discussed in § 5.

## 2 Numerical Method

### 2.1 The Numerical Code

The models discussed below have been computed using Cosmos, a massively parallel, multidimensional, radiation-chemo-hydrodynamics code for both Newtonian and relativistic flows developed at Lawrence Livermore National Laboratory. The relativistic capabilities and tests of Cosmos are discussed in Anninos & Fragile (2003). Tests of the Newtonian hydrodynamics options and of the internal physics relevant to the current work are presented in Anninos, Fragile & Murray (2002), and so we shall not discuss those in detail here.

Because we shall be examining the effects of having multiple supernova explosions at random locations throughout the cores of dwarf galaxies, we cannot make any assumptions regarding the symmetry of the systems. As a consequence, our models are run in three-dimensions on a Cartesian mesh. Most of the simulations are run with resolutions of $`128^3`$ zones.

### 2.2 Radiative Cooling

Chemistry is not followed in these simulations. Instead, local cooling is given by the following cooling function, based in part on an equilibrium (hydrogen recombination and collisional excitation) cooling curve

$$\mathrm{\Lambda }(T,\rho )=\left[\underset{i}{}\dot{e}_i(T)(f_I\rho )^2+\dot{e}_M(T)(f_M\rho )^2\right]\times \{\begin{array}{cc}e^{(TT_{min})/\delta T}\hfill & \text{if}TT_{min},\hfill \\ 1\hfill & \text{otherwise},\hfill \end{array}$$
(2-1)
where $`\dot{e}_M`$ is the temperature-dependent cooling rate from metals (including carbon, oxygen, neon, and iron lines, assuming a metallicity $`Z=0.03Z_{}`$), $`\dot{e}_i`$ is the cooling rate from hydrogen and helium lines, $`f_M`$ is the mass fraction of metals taken to be $`6\times 10^4`$, $`f_I`$ is an estimate of the ionization fraction, defined as $`\text{min}(1,\text{max}(0,(T_{eV}T_c)/3))`$ with $`T_c=1eV`$ to match roughly the upper and lower bounds in the equilibrium model described in Anninos, Fragile & Murray (2002), $`T=(\mathrm{\Gamma }1)e\mu m_p/(k_B\rho (1+f_I))`$ is the gas temperature in Kelvin, $`\rho `$ is the gas density, $`e`$ is the internal energy density of the gas, $`\mu `$ is the mean molecular weight, taken to be unity, $`m_p`$ is the proton mass, and $`k_B`$ is Boltzmann’s constant. The exponential is introduced to suppress cooling at low temperature (here $`T_{min}=10^4`$ K, with width $`\delta T=500`$ K). The initial temperature of the gas in the systems considered here is 8200 K, representing a state which is partially ionized by the background. The rapid decrease of cooling below $`10^4`$ K simulates the effects of heat sources and the exponential decay of electron density below about 1 eV in equilibrium models. This prevents the gas from radiatively cooling significantly below its initial temperature.

Our choice of metallicity is motivated by the fact that we are examining the evolution of young dSph’s, which have not had a chance to undergo significant self-enrichment. The value $`Z=0.03Z_{}`$, is appropriate to the older generations of stars observed within dSph’s (Grebel, 2001).

### 2.3 Gravity

The galaxies are modeled as non-self gravitating gas within a fixed dark-matter potential. The omission of self-gravity is reasonable, given that the gravitational potential of dwarf galaxies with $`M_g10^9M_{}`$ is dominated by the dark-matter halo (Persic, Salucci, & Stel, 1996), where $`M_g`$ is the mass of the visible (gas plus stellar) matter. The shape of the potential is given by a modified Hubble profile (Binney & Tremaine 1987)

$$\varphi (r<R_t)=\frac{G\stackrel{~}{M}}{R_c}\left\{1\frac{\text{ln}[x+(1+x^2)^{1/2}]}{x}\right\},$$
(2-2)
where $`x=r/R_c`$, $`R_c`$ is the core radius, and

$$\stackrel{~}{M}=M_d\left\{\text{ln}[x_t+(1+x_t^2)^{1/2}]x_t(1+x_t^2)^{1/2}\right\}^1,$$
(2-3)
where $`M_d`$ is the dark-matter mass, $`x_t=R_t/R_c`$, and $`R_t`$ is the tidal radius. This potential has a simple form, while retaining the essential behavior for the shape of potentials which are generally observed in dSph’s (Burkert 1995). The gas is initialized to be isothermal and in hydrostatic equilibrium within the potential, such that the density

$$\rho =\rho _0e^{\varphi /c_s^2},$$
(2-4)
where $`c_s`$ is the isothermal sound speed. The value of $`\rho _0`$ is scaled to give the desired gas mass $`M_g`$ within $`R_t`$. The effects of a tidal radius are included by flattening the potential beyond $`R_t`$. Beyond $`R_t`$, the gas density is also decreased by a factor of 0.01, and the temperature is correspondingly increased in order to maintain hydrostatic equilibrium. This simplified treatment of the galaxy’s surroundings allows us to focus on the physics within the dwarf galaxy itself.

### 2.4 Boundary Conditions

To minimize the effects of the computational boundaries on the galaxy evolution, we choose the size of the grid, $`l`$, such that $`2R_t=0.8l`$. Our boundary zones remain static throughout these simulations, maintaining their initial density and temperature. This choice of boundary conditions resulted in fewer disruptive reflections than flat (vanishing first derivative) boundaries and less spurious mass loss than outflow boundaries. We have verified that, in the absence of supernovae, the models remain static for several dynamical times.

### 2.5 Galaxy Parameters

Two different model galaxies have been examined. Their parameters are given in Table 1, which lists the dark-matter mass $`M_d`$, gas mass $`M_g`$, core radius $`R_c`$, the central gas density $`\rho _0`$, the dynamical time of the galaxy $`\tau _d`$ (defined below), the total gravitational binding energy of the gas, $`\mathrm{\Phi }`$, and the total thermal energy of the gas, $`E_{th}`$, for each model. The parameters of galaxies 1 and 2 are chosen to be close to those of Leo II and Fornax, respectively, objects near the ends of the distribution for dSph’s (Mateo 1998). For both model galaxies, the ratio of the tidal to the core radii is chosen to be $`R_t/R_c=3`$. We note that the central densities of the two models are similar, in spite of their having different masses, due to the larger core radius of Model 2. In Figure 1, we plot the square root of the absolute value of the normalized gravitational potential ($`|\varphi \varphi _t|^{1/2}`$) as a function of radius for the two galaxy models. This function gives a good estimate of the local escape velocity.

It may also be noted from Table 1 that the total binding energy of the gas in the two models is of the order of the energy input of 1-10 supernovae, which is precisely the range covered in our models (see below). It might therefore be expected that the energy input from supernovae would lead to extensive mass loss from the systems. As shall be seen, however, this is only the case when radiative cooling is neglected. As has been noted in earlier studies (Larson, 1974), radiative losses cause the deposition of energy from supernovae to be relatively inefficient. We discuss this point further and examine the effects of radiative losses in § 3.2.

### 2.6 Treatment of Supernovae

In this study, supernovae are simulated by adding 1051 erg of internal energy to the gas over finite regions. The internal energy is injected over an approximately spherical volume with a “top-hat” profile of 3-4 zones radius. Computing time constraints limit our resolutions, and we are therefore unable to maintain zoning sufficiently fine as to resolve the supernova shocks and the cooling regions immediately behind them. The results are therefore vulnerable to numerical over-cooling of the shocked gas, which would lead to an underestimate of the effects of the supernovae. This potential problem is exacerbated by having the energy initially deposited into the internal energy of the gas, after which it evolves to a Sedov-Taylor solution. In addition, the artificial viscosity used to capture shocks spreads out shock heating over a few zones. Gas may, therefore, spend multiple timesteps at elevated temperatures, possibly near the peak of the cooling curve. This can cause further over-cooling of the gas, significantly affecting the hydrodynamic evolution. In order to prevent over-cooling of the gas we simultaneously employ two techniques. First, we do not allow cooling within a supernova until is has roughly doubled in radius. Second, we follow the common procedure of preventing cooling of the gas during shock passage (Shapiro & Struck-Marcell, 1985; Blondin, Fryxell, & Konigl, 1990; Anninos & Norman, 1994). To accomplish this in a three-dimensional problem, where shocks may occur in any direction, a simple scalar criterion is needed. We find, in practice, that the ratio of the artificial viscous heating term to gas pressure is an adequate filter. Thus, cooling is turned off whenever $`Q/P>0.1`$. The effects of these techniques will be explored in § 3.2.

In order to follow the evolution of the supernova-enriched material, we add a tracer fluid to the zones in which the supernova energy is initially deposited. The tracer is a passive quantity that is advected in the same manner as density.

The simulations discussed below include either one or ten supernovae. In the cases with only a single supernova, it is set off either at the center of the system or at the core radius. In these cases we are able to take advantage of symmetries in the problem to reduce the required computational resources. For a single supernova set off at the center of the galaxy, we employ $`x`$, $`y`$, and $`z`$ symmetries and run a single octant with reflective boundary conditions at the $`x=0`$, $`y=0`$, and $`z=0`$ planes. For a single supernova at $`R_c`$, we can invoke symmetry in two dimensions and run a quadrant with reflective boundaries along two planes. If ten supernovae are included, the first is set off at $`R_c`$ at the start of the simulation. The rest are chosen to occur at random orientations to the center of the system. The radii at which they occur are chosen at random, weighted by the enclosed mass of gas, and constrained to be within $`R_c`$. The latter constraint is due to the expectation that bursts of star formation are expected to be located predominantly in galactic cores, where gas densities are highest (Dong, Murray, & Lin, 2003). The timescale over which supernovae occur is determined by the timescale of the burst of star formation that led to the formation of the progenitor stars, and by the spread in masses, and therefore lifetimes, of the progenitor stars. To examine the effects of reasonable variations in these factors, the supernovae in the models are set off at times chosen randomly within total spans of either 10 or 100 Myr.

### 2.7 The Form of the Interstellar Medium

In most of our models, the gas density varies smoothly, following Equation 2-4. Following a burst of star formation in a dwarf galaxy, however, such a smooth density distribution is unlikely to actually be found. Instead, ionizing photons and winds from massive stars, and turbulent motions shall all lead to a very irregular gas distribution for the interstellar medium (ISM). The low density regions in such a distribution may enhance the loss of material from supernovae, and so it might be expected that the form of the gas distribution could play a role in the ability of a galaxy to retain enriched gas. In order to explore the role of the gas distribution, two of the models include 100 “bubbles,” each five zones in radius, whose density is decreased by a factor of five relative to the surrounding gas. The physical interpretation of these bubbles is that they are cleared by ionizing photons and winds from massive stars. The temperature within the bubbles is increased, such that they remain in pressure equilibrium with the surrounding gas. The bubbles are selected to have random orientations around the center of the galaxy and they are not, in general, co-located with the supernovae, which is justified by the fact that during their main-sequence life span, massive stars may move up to 100 pc from their initial position. The distances of the bubbles from the galaxy center are selected at random, weighted by enclosed mass and constrained to be no more than 0.6$`R_t`$ from the center of the galaxy. Once the bubbles are placed, the density of the entire system is increased by a common factor to preserve the original mass of the galaxy. To simulate the presence and evolution of the hot stars causing the bubbles, cooling is not allowed within the bubbles for 0.5 dynamical times after the supernovae end.

In physical units, the bubble radii correspond to approximately 40 pc for Model 1 and 100 pc for Model 2, comparable to the radii of Strömgren spheres. The total filling factor for the bubbles within $`0.6R_t`$ is approximately 30%, i.e. the ISM is highly irregular in form. For typical stellar initial mass functions, the number of supernovae and bubbles in our models represent overall star formation efficiencies for the galaxies of less than 1%.

### 2.8 Other Physics

The models discussed below do not include the effects of conduction or of magnetic fields. The effects of thermal conduction may be estimated from the conductive flux

$$F_c=\kappa T^{5/2}T,$$
(2-5)
where $`\kappa =6\times 10^7`$ erg cm-2 s-1 K-7/2 (Spitzer, 1956). The total energy flow across the interface between an expanding bubble of heated gas may be compared with the radiative energy losses from within the volume of the bubble. Radiative losses rapidly lead to a relatively shallow temperature gradient, greatly reducing the conductive flux, with the result that conductive losses are found to be much less important than radiative emission in our models, justifying the omission of conduction.

The presence of even weak magnetic fields would greatly hamper conduction, further reducing its importance (Rosner & Tucker, 1989). Stronger magnetic fields could potentially affect the hydrodynamics. A strong, ordered magnetic field would tend to direct the flow of gas to be one-dimensional, whereas a strong field that was tangled on small scales would act as an additional pressure term, possibly helping to contain gas within the galaxies. The magnetic field would be expected to noticeably affect the flow if the magnetic pressure were comparable to the ambient gas pressure, which, for our models, requires $`B5\mu `$G, greater than or comparable to values measured in the local ISM (Heiles, 1987). Whether dwarf systems could generate fields with strengths comparable to those in our own Galaxy is uncertain, as are the actual field strengths present today in dwarf galaxies. Magnetic fields therefore remain an uncertain quantity. If present in young dSph’s at magnitudes comparable to the local ISM, then their affect would be expected to be comparable to that of the ambient gas pressure, having a quantitative affect upon our results, but not affecting the qualitative conclusion that enriched material can be preferentially lost from dwarf systems.

## 3 Model Results

The parameters of the models that we have computed are shown in Table 2. Listed in the table for each model are the galaxy model used, the number of supernovae included $`N_{SN}`$, the timescale over which the supernovae occur $`\tau _{SN}`$, and the number of bubbles included in the galaxy $`N_{\mathrm{bub}}`$. Variations of many of these base models are also examined (see below). All of the runs are summarized in Table 3, where the variations are indicated by numbers following the letters of the base models. The tabulated results include the percentage of tracer lost from the core as well as that lost beyond $`R_t`$, and also the percentage of gas mass lost from the galaxies. The core tracer loss is measured after at least 2.5 dynamical times for the galaxies, while the tracer and mass losses beyond $`R_t`$ are estimated from the global minima (see below). In this section, we discuss the results of the models, and we follow with a discussion of the reasons for the differences among them in § 4.

### 3.1 Single Supernova Events

In Models A1 and B1, single supernovae are set off in the centers of the two model galaxies. For both models, all of the tracer remains within the core of the galaxy. The lack of tracer loss is expected, given the core masses of the low and high mass models, which are approximately $`1.6\times 10^5`$ and $`3\times 10^6`$ M, respectively. Such large masses of gas have no trouble containing the ejecta from single supernova events.

We have also run versions of Models A and B, in which a single supernova was triggered at $`R_c`$ (Models A2 and B2, respectively). Although this obviously resulted in significant loss of tracer from the core, only the low mass galaxy (Model 1) showed any loss of tracer from the system (30%). We conclude from this that multiple supernovae appear to be required in order for large fractions of enriched material to be lost from within rounded galaxies, even as small as those considered here.

### 3.2 Multiple Supernovae in Undisturbed Systems

#### 3.2.1 Low Mass Galaxy Models

In Models C and D, 10 supernovae are set off at random locations within the core of the initially undisturbed galaxy Model 1. The supernovae are triggered at random times over a span of 10 Myr for Model C and 100 Myr for Model D.

We further utilize Model C to test the effects of resolution upon our results and to explore the impact of different choices for radiative cooling and metallicity. Model C1 is our control model, as it utilizes our standard cooling options, resolution, and metallicity. Hence, this model is most comparable to Models A, B, D, E, F, G, and H. Model C2 utilizes the same cooling options, but has a spatial resolution reduced by a factor of two. Model C3 tests the effect of not shutting off cooling at shock fronts by ignoring the $`Q/P>0.1`$ criterion described in § 2.6. Similarly, Model C4 tests the effect of not shutting off cooling in the initial supernova region by ignoring that requirement. Model C5 ignores both the $`Q/P`$ condition and the requirement to shut off cooling in the initial supernova region, and thus represents results with cooling turned on “everywhere,” and so likely suffers from significant over-cooling. Model C6 tests the effect of our choice of metallicity, by using a metallicity $`Z=0.3Z_{}`$, an order of magnitude higher than that of Model C1. Finally Model C7 represents the extreme case in which radiative cooling is turned off everywhere, and the gas may cool only by adiabatic expansion. We proceed by first comparing the results of Models C1 and D. The effects of the variations of Model C are discussed in § 3.2.3.

In Figure 2 are shown snapshots of two-dimensional slices of the density evolution through the central galactic plane for Model C1. The snapshots cover a span of approximately 100 Myr. The evolution of one of the tracer fluids from Model C1 is shown in Figure 3, at times corresponding to the frames in Figure 2. This tracer is associated with the initial supernova (triggered at $`R_c`$), visible in the second frame of Figure 2. Figure 4 shows the time evolution of the amount of tracer contained within $`R_t`$ (solid curves) and $`R_c`$ (dashed curves) for Models C1 and D. The enriched material from each supernova is labeled with a different tracer, and the evolution is shown individually for each. On the curves, time is presented in code units, where a unit of one represents approximately the dynamical time of the galaxy, defined as

$$\tau _d\frac{R_t}{c_s},$$
(3-1)
where $`c_s`$ is the sound speed of the gas. The time unit for galaxy Model 1 is 48.2 Myr.

Significant ($`80`$%) loss of tracer is seen from within $`R_c`$ for both models, although Model C1 loses a greater fraction than does Model D. In Model C1, the loss of tracer is accompanied by an initial loss of approximately 50% of the gas from the core. The core mass returns to its initial value, however, within $`\tau _d`$. This recovery of the initial core mass without an accompanying recovery of tracer implies that the tracer does not mix well with the gas in the core. In Model D, where the supernovae occur over a longer timescale, the evolution of the core mass is much gentler, yet the same conclusion holds.

Additionally, in Model C1, nearly one-half (47%) of the tracer moves beyond $`R_t`$, while almost 20% does so in Model D. This material also does not mix well with the surrounding gas, as shown in Figure 3, where the tracer can be seen to evolve as a nearly cohesive unit. As can be seen in Figure 4 (and subsequent figures), some of the tracer that moves beyond $`R_t`$ later returns to the galaxy. This behavior is an artifact of our boundary conditions. At the boundaries, the gas properties are held constant in time, which can artificially prevent material from leaving the grid, and drive it back to the galaxy. In actuality, material that moves beyond the tidal radius is unlikely to return to the dwarf galaxy on the time scales shown here, but shall instead follow a separate orbit around the parent galaxy. The gas which escapes from the galaxy in Model C1 does so with velocities of only 5-7 km s-1, much less than the orbital speed of the dSph satellites around the Milky Way. During an orbit, therefore, the gas would not move substantially away from the dSph, and would possibly re-encounter it on the next apogalacticon passage (Dong, Lin, & Oh 2002).

The loss of tracer from the galaxy is accompanied by overall mass loss, as can be seen in Figure 5, which shows the mass contained within $`R_t`$ (solid curves) and $`R_c`$ (dashed curves) for Models C1 and D. The inflow due to the outer boundary conditions is apparent in both figures. As discussed above, in the absence of a confining external medium, gas is actually unlikely to return to the system after it moves to beyond $`R_t`$. The total mass loss to be expected in the models may therefore be represented by the first minima shown in Figure 5, which represent 36% and 21% for Models C1 and D, respectively.

#### 3.2.2 High Mass Galaxy Models

Models E and F were similar to C and D, except that the higher mass galaxy model (Model 2) was used. The evolution of the tracer fluids for these models is shown in Figure 6, displayed as in Figure 4. The time unit for galaxy Model 2 is 130 Myr. A large fraction of the tracer is lost from the core, though less than seen in Models C1 and D, and almost no tracer is lost beyond $`R_t`$. Such differences are expected given the larger mass of this galaxy model. Little mass loss is seen from the system in Models E and F, as expected from the tracer evolution.

#### 3.2.3 The Effects of Numerical Techniques

As discussed above, we tested the effects of resolution upon our results by reducing the spatial resolution of Model C2 by a factor of two relative to Model C1. In order to keep other factors as similar as possible, the physical radii over which the supernova energies were added was kept the same, requiring that the energy be added over many fewer zones. The results were similar to Model C1, with almost identical amounts of tracer loss from the core and overall mass loss. Somewhat less tracer was lost, however, beyond $`R_t`$ in the lower resolution model (33% as compared to 47%). Such a difference would be expected if the lower resolution allowed more cooling in the supernova-heated gas, reducing its buoyancy.

In Model C3, cooling was allowed at shock fronts, but was still prevented during the initial expansion of the hot gas. As compared to Model C1, Model C3 lost noticeably less tracer and mass beyond $`R_t`$, with at least 25% reductions relative to Model C1. These results confirm our suspicion that allowing cooling in the under-resolved shock fronts results in a noticeable increase in radiative losses, though the effect upon the evolution is not overly large.

The results of Model C4, in contrast, suggest that the technique of turning off the cooling in the initial supernova region until it doubles in radius has only a small effect upon the final results. This result implies that the gas radiates away relatively little energy before the Sedov-Taylor solution is reached.

In Model C5, radiative cooling is allowed everywhere and at all times. As would be expected from the above results, the results of Model C5 are very similar to those of Model C3.

Model C6 is similar to Model C1, except that we assume a metallicity for the gas of $`Z=0.3Z_{}`$, ten times larger than that of Model C1. As can be seen in the table, this has no noticeable impact upon the results. From this, we conclude that, by the time the post-shock gas temperature drops to below $`10^4`$ K, where metal cooling dominates, the evolution of the gas is determined by the momentum given to it by the shock, and is no longer affected by shock-heating.

In Model C7, we examine the effects of turning off radiative cooling completely, i.e. the gas evolves adiabatically from its initial conditions. As would be expected from earlier work (Larson, 1974), the differences with all of the other models are dramatic. In Model C7, most of the original gas in the galaxy ($`80`$%) is expelled due to the energy input of the 10 supernovae, consistent with the estimate of the binding energy of galaxy Model 1. Similarly, almost all of the tracer material is lost.

From these results, we conclude, consistent with earlier studies, that radiative losses are crucial to the ability of dwarf galaxies to retain gas following supernova explosions. Numerical techniques employed to improve the accuracy of simulations, by turning off cooling in regions where overcooling is likely to occur, have quantitative effects upon the results, but do not affect qualitative conclusions that significant amounts of both tracer and gas are lost to the systems. The metallicity of the gas is found not to be a significant factor.

### 3.3 Multiple Supernovae in Systems with an Irregular ISM

Models G and H use, respectively, the low- and high- mass galaxy models. As in Models C and E, 10 supernovae are triggered over a timescale of 10 Myr. They differ from those earlier models in containing 100 low-density bubbles, which are intended to simulate an irregular ISM, as described in § 2.7. In Figure 7 are shown snapshots of density evolution in the central galactic plane for Model G, covering the same time span as those in Figure 2. The frothy nature of the initial gas distribution is clearly illustrated in the first few frames. As described in § 2, the bubbles are allowed to cool after $`\tau _{SN}+0.5\tau _d`$ (34.1 Myr for galaxy Model 1 and 75 Myr for galaxy Model 2). The cooling evolution is apparent in the eleventh frame of Figure 7. In Figure 8 is shown the evolution of the same tracer fluid as in Figure 3, now for Model G.

The evolution of the tracer and the mass for the two models are shown in Figures 9 and 10. The reduction of the average density in the cores of the models, and the presence of low density channels lead to dramatic differences from the earlier models. Nearly all of the enriched material is lost from the cores of both models. Approximately 70% is lost beyond $`R_t`$ as well for Model G, while approximately 50% is lost beyond $`R_t`$ for Model H. The systems also show mass losses of 40-50%.

As in Model C, the loss of tracer from the cores of Models G and H is accompanied by mass loss. Unlike the tracer loss, the mass loss later reverses, again implying poor mixing of the enriched gas within the galactic cores. Unlike the earlier model, the core masses eventually overshoot their initial values. The overshoot is due to the cooling of the bubbles, which causes some inflow into the central regions. As before, the enriched material which is lost beyond $`R_t`$ is poorly mixed, as can be seen in Figure 8.

## 4 Differences Between the Models

Models which included only a single supernova generally showed little loss of tracer material. In the case of supernovae centered in the galaxies, no tracer was lost even from the core of the galaxy. In the case of off-center supernovae, only the low mass system lost some of the enriched gas beyond the tidal radius. Such a lack of any effect may be surprising, given the weak gravitational binding energy of the gas, $`\mathrm{\Phi }`$, in the two model galaxies, as shown in Table 1. In Model 1, for example, $`\mathrm{\Phi }=9.8\times 10^{50}`$ ergs, comparable to the energy of a single supernova. As discussed below, however, the relatively high gas densities of the systems allow them to radiate away the deposited energy relatively efficiently, with the result that a single supernova has little overall effect upon the systems.

The occurrence of multiple supernovae enhances the loss of enriched material in a few different ways. Early supernovae disturb the gas in the core of the galaxies, reducing the average density. The cooling efficiency of the gas is therefore reduced, as is the mass of gas available to physically contain the ejecta from later events. Later events may also help to give a “boost” to the enriched material from earlier events, as can be seen in Figure 3, being especially apparent at the times 3.8 and 7.7 Myr. Note that the apparent acceleration in the evolution of the tracer between 24.1 and 33.7 Myr is an illusion caused by the increasing time interval between frames.

The differences among the models with multiple supernovae are due to their different masses, and their relative abilities to cool within a dynamical timescale. As shown in Table 1, Models 1 and 2 have very similar central densities. Their different core radii therefore lead directly to different core masses, as discussed in § 3.1, with Model 2 having a greater amount of mass available to contain supernova ejecta than Model 1.

In the absence of cooling, the hot, enriched material from the supernovae would still eventually leave the galaxies due to their buoyancy. Radiative cooling prevents this, and the different size scales for the two galaxy models lead to very different ratios of cooling to dynamical times. The dynamical timescales from Equation (3-1) are $`\tau _d=`$ 47 and 128 Myr for galaxy Models 1 and 2, respectively. The cooling timescale is defined as (Field 1965)

$$\tau _c\frac{\frac{3}{2}k_BT}{n_HL(T,Z)},$$
(4-1)
where $`L=\mathrm{\Lambda }/n_H^2`$ is the cooling efficiency in erg cm3 s-1, and $`n_H`$ is the number density of hydrogen. Close to their initial temperatures, and using our assumed metallicities, $`\tau _c`$ 36 and 39 Myr, for Models 1 and 2 respectively. The initial cooling of the supernova-heated gas shall be much more efficient, with a correspondingly shorter timescale. For the gas to cool to its initial temperature, however, shall require timescales on the order of those calculated above. For Models 1 and 2, we find, respectively, $`\tau _c/\tau _d=`$0.77 and 0.30. The values for $`\tau _d`$ would be smaller, and the timescale ratios larger by a factor of three had we used $`R_c`$ instead of $`R_t`$ in Equation (3-1).

The higher mass galaxy model is therefore able to cool much more on a dynamical timescale than is the lower mass model. Following the conversion of the deposited energy into kinetic energy, the high mass model can thus radiate away much of the energy before the enriched material can evolve significantly in space, further enhancing the ability of the more massive system to retain the enriched gas. This result is consistent with the observed correlation between metallicity and total mass in dSph’s (Mateo, 1998).

The above factors work together to explain the differences seen between the pair of Models C1 and D, which use the low-mass galaxy model, and the pair of Models E and F, which use the high-mass galaxy model. In the first pair, almost all of the enriched gas is lost from within $`R_c`$, and both lose a significant amount beyond $`R_t`$. Both models also show a substantial amount of mass loss. In Models E and F, however, less enriched material is lost from within $`R_c`$, almost none is lost from $`R_t`$, and there is much less mass loss.

The ability of the galaxy to cool between supernova events can also explain the differences between Models C1 and D, seen in Figure 4. In Model C1, the supernovae occur over a total time of 10 Myr, and the galaxy is unable to cool significantly between events, whereas the longer 100 Myr supernova timescale of Model D allows the system to radiate a significant amount of energy between supernovae. As compared to Model D, then, the supernovae have a much stronger cumulative effect in Model C1, leading to more rapid loss of enriched material from the core, as well as loss beyond $`R_t`$. The implication of these models is that star-burst events do not necessarily lead to greater enrichment of residual gas than does continuous star formation.

The differences between Models E and F are less striking than those between C1 and D, due primarily to the fact that the more massive system is able to retain enriched material as a whole more efficiently than the less massive system. As can be seen in Figure 6, there is less loss from the core of enriched material from a few of the supernovae in Model F as compared to Model E, which primarily accounts for the quantitative differences between the two models. Those supernovae happened to occur at late times, as the system was “recovering” from the earlier supernova events.

Models G and H were identical, respectively, to Models C1 and E, except for the presence of bubbles. The large filling factors of the bubbles in the models dramatically altered the evolution of the tracer. The bubbles reduce the amount of mass in the core of the galaxy, and increase the average cooling time of the galaxy. They also provide many channels for enriched gas to exit the core. As a result, most of the tracer is lost beyond $`R_c`$ in both models, and a significant fraction is lost beyond $`R_t`$. Both models also lose 40-50% of the ambient gas. This result can be used to infer an upper limit on the rate of massive star production in order for there to be sufficient retained gas to fuel subsequent generations of star formation.

## 5 Consequences for the Evolution of Dwarf Spheroidal Galaxies

We have performed three-dimensional, numerical simulations of the evolution of enriched material from supernovae in dwarf spheroidal galaxies. The simulations have included either one or ten supernovae, and in the case of multiple supernovae, they have occurred over timescales of either 10 or 100 Myr. Two different galaxy models, 1 and 2, have been considered, chosen to approximate the parameters for the Leo II and Fornax dSph’s, respectively. Within these galaxy models, both smooth and irregular ISM distributions have been considered. The results of the models, discussed in § 3 and 4, are summarized in Table 3. We now discuss the possible consequences of our results for the evolution of dSph’s.

Our results indicate that, for extremely low star formation efficiencies, such as would lead to individual supernovae widely separated in time, it is difficult to expel enriched material completely from dwarf spheroidal galaxies (dSph’s) having an undisturbed ISM. In lower-mass systems, however, most enriched gas is lost from the cores of the galaxies following multiple supernovae, and a significant fraction is lost to the galaxy as a whole. In a system that has undergone a recent burst of star formation, however, the ISM is likely to be highly disturbed. We find that the presence of an irregular ISM greatly enhances the loss of enriched material both from the core as well as from the galaxy as a whole, in both low and high mass models. In all cases, the material is lost quickly, leaving the core within $`0.2\tau _d`$ after the supernovae. Enriched material that is lost to the galaxy is typically lost within $`2\tau _d`$. As can be seen in the tracer plots, enriched gas does not mix well prior to being lost to the galaxy, leading to highly enriched “blobs” of material escaping beyond $`R_t`$.

The preferential loss of enriched gas from the cores of dSph’s might potentially lead to an “inverse” metallicity gradient, the consequences of which for subsequent star formation depend upon where later generations of stars form. If subsequent star formation is triggered by relatively quiescent processes, such as tidal compression, or runaway cooling resulting from self-shielding from external radiation, then star formation is expected to occur primarily in regions where the gas density is highest, i.e. in the core. In such cases, star formation shall occur in gas which is not heavily enriched from the previous generation of stars, with the result that multiple generations of stars could be present in the galaxy without strong metallicity variations. Given that the ISM is likely to be highly disturbed by the presence of hot stars, our models indicate that, in this case, poor self-enrichment efficiency is likely to be the norm in dSph’s.

Star formation triggered by localized compressions may, in contrast, occur at larger radii, where the gas may have been previously enriched. This might occur, for example, due to collisions between dwarf systems, encounters with a jet from a nearby massive galaxy, or due to the infall of material onto the dwarf galaxy. The latter case might well be the result of infall of material that was originally expelled from the dwarf itself. When material in our simulations is lost from the galaxies, it exits at below the sound speed ($`10`$ km s-1), much less than the orbital speed of dwarf galaxies around systems such as the Milky Way. The ejected gas shall, therefore, follow a similar orbit around the parent galaxy, and may collide with gas remaining within the dwarf system due to orbit crossing at apogalacticon (Dong, Lin, & Oh 2002). Such an encounter might lead to star formation at large radii in the dwarf, and a subsequent generation of stars showing significant metallicity enhancement.

Enriched gas which accumulates in the outer regions of dwarf galaxies may eventually be lost to the galaxy through ram-pressure stripping, tidal effects, or due to subsequent generations of star formation in the core. Thus, as much as 90% of the supernova ejecta could be lost, and may enrich the halo of a nearby massive galaxy. The halo of the Milky Way, for example, could have been enriched to its average metallicity by a population of dwarf spheroidals ten times larger than the number currently observed.

This work was performed under the auspices of the U.S. Department of Energy by University of California, Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. This work is partially supported by NASA through an astrophysical theory grant NAG5-12151.


text:
# Detectibility of Mixed Unburnt C+O in Type Ia Supernova Spectra

## 1 Introduction

Interest in Type Ia supernovae has been high over the last twenty years, since it was recognized that they closely meet the astronomer’s criteria for a “standard candle”. This interest has grown with the realization that one can use the observed light curve to make a luminosity correction (similar to the Cepheid P-L relation, Phillips, 1993; Riess et al., 1995, 1996; Phillips et al., 1999). SNe Ia played an important role in current determinations of the Hubble constant with the *Hubble Space Telescope* (Sandage et al., 1996; Branch, 1998) With the discovery of the “dark energy” (Riess et al., 1998; Perlmutter et al., 1999; Goldhaber et al., 2001) enthusiasm for studies of Type Ia supernovae has grown. Heretofore, all calibrations of the luminosity of SNe Ia have been empirical in nature; relying largely on the sample obtained by the Calan-Tololo survey (Hamuy et al., 1996). While significant efforts are underway to improve the local sample of SNe Ia, theoretical understanding of the explosion mechanism and progenitor system are required in order to produce confidence in the validity of the empirically based results.

In spite of the fact that the progenitor system for SNe Ia is still uncertain (Branch et al., 1995), there is little doubt that it requires the thermonuclear explosion of a C+O white dwarf. Given this fact, progress is being made in the understanding of the explosion mechanism itself and calculations have progressed from 1-D and 2-D, to full 3-D (Khokhlov, 1993, 2001; Hillebrandt & Niemeyer, 2000; Reinecke et al., 1999, 2002c, 2002a, 2002b; Gamezo et al., 2003)

The calculation of the explosion mechanism is a formidable numerical problem. There is good agreement that the explosion begins as a deflagration with a laminar flame, which then develops into a fully turbulent flame via various instabilities, dominated by Rayleigh-Taylor instabilities. Even in 3-D, the flame cannot be fully resolved and a variety of subgrid models have been developed. Recent work has focused on calculations of the deflagration of C+O Chandrasekhar mass white dwarfs. While some calculations have had trouble obtaining the observed energy of the explosion, there seems to be general agreement that a fully resolved 3-D calculation can produce significant kinetic energy and enough nickel to reproduce the general features of a Type Ia explosion. Recently, Gamezo et al. (2003) performed a series of 3-D deflagration model calculations in which they obtained a total kinetic energy of about $`10^{51}`$ ergs. The key feature of the simulations is the highly convoluted turbulent flame surface that allows extensive interpenetration of burnt and unburnt materials. Similar results were found by Reinecke et al. (2002a, b). However, the strong interpenetration of burnt and unburnt material leaves a significant amount of unburnt C+O in the central regions of the ejecta. This result led Gamezo et al. (2003) to conclude that deflagration models would be in conflict with observed supernova spectra because there would be strong carbon or oxygen lines produced by the model, which are not seen in observed SNe Ia spectra. Based on this inference, Gamezo et al. (2003) rule out the deflagration model completely as a viable model for the explosion mechanism of SNe Ia. In drawing this conclusion Gamezo et al. (2003) rely on the work of Jeffery et al. (1992), Fisher et al. (1997), Mazzali (2001), and Kirshner et al. (1993) which all use a highly parameterized direct spectral analysis technique for producing synthetic spectra. In particular, all these analyses make the Schuster-Schwarzschild approximation and thus cannot address the formation of lines below the imposed photosphere. These analyses are extremely useful for line identifications in observed spectra and for determining the velocity interval at which certain ions *must* exist. However, these analyses cannot (and do not claim to) determine whether a given species must be absent. That is, from these analyses one can conclude that a certain ion must be present in a certain velocity interval, but one cannot conclude that any particular ion must be absent.

In order to examine the true viability of a full 3-D hydrodynamical calculation, a full 3-D radiative transfer calculation should be performed. However sophisticated 3-D radiative transfer codes do not yet exist that can tackle the problem of detailed radiative transfer in SNe Ia. In order to examine the viability of deflagration models we present the results from PHOENIX calculations using a parameterized version of the tuned 1-D deflagration model W7.

## 2 Calculations

The W7 model (Nomoto et al., 1984; Thielemann et al., 1986) is a good match to normal SNe Ia spectra (Harkness, 1991a, b; Nugent et al., 1995a, b; Wheeler et al., 1995; Lentz et al., 2001), especially the outer layers visible during the pre-maximum phase of the supernova. The outer layers of W7 are unburnt C+O enriched with solar abundance metals, expanding at velocities, $`v>15000`$ km s-1, with a total mass of $`0.07`$ $`\text{M}_{}`$, followed by intermediate mass elements, with an interior composed mostly of 56Ni but also including some other iron group elements. The model is a 1-D deflagration model, where the flame speed was adjusted to produce intermediate mass elements at high velocity as observed in SNe Ia spectra.

We have altered the inner composition of the W7 to produce two alternate models with C+O mixed into the center. The models are prepared by removing a portion of the original mass in the W7 model and replacing it with an equal mass of C+O, with the C/O ratio = 1 by number, keeping the density structure the same. The C+O is used for the unburnt fuel mixed to the core found in 3-D deflagration models. In the first model, labeled “25% C+O core”, we have replaced 25% of the inner core with C+O. The mixed material extends from the center to 3000 km s-1, where the mixed fraction is then reduced linearly in velocity space to 0% at 4000 km s-1. In the second model, designed to mock-up the Gamezo et al. (2003) model and labeled “50% C+O core”, 50% of the mass is replaced with C+O from the center to 4000 km s-1. The mixed fraction is reduced linearly in velocity space to 10% at 5000 km s-1. The 10% C+O region extends to 15000 km s-1 where the W7 model becomes entirely unburnt material. Examining Figure 3 of Gamezo et al. (2003) it is diffiult to determine the exact composition structure, since they only distinguish between unburnt (C+O) and “burnt” material, but the exact composition of the burnt material is unspecified because it is too computationally intensive to follow the reaction rates in detail.

The spectra and $`\gamma `$-ray deposition are calculated using the general purpose radiative transfer code PHOENIX, version 11.9 (Hauschildt & Baron, 1999). There is some change in the total gamma-ray deposition with the composition change, but the change is smaller than one might expect, since in this velocity range most of the material in W7 is neutron-rich iron peak elements and not 56Ni. For the 25% C+O core model the change is 2%, and for the 50% C+O core model the change is 20% at both epochs. PHOENIX includes all of the effects of special relativity in the steady-state solution of the transport equation and energy balance. For these exploratory calculations we have treated all species in LTE. Studies by Baron et al. (1996) and Wheeler et al. (1998) have shown the LTE is a reasonable approximation in SNe Ia and since we are looking only qualitatively for lines that are clearly absent in observed spectra the assumption of LTE should not affect our conclusions. The use of PHOENIX with W7 models is described in detail in Lentz et al. (2001).

## 3 Results

Figure 1 displays the observed spectrum of SN 1994D on 1994 April 14, roughly 15 days past maximum light (35 days past explosion, Patat et al., 1996) and the results of three calculations: 1) the original W7 model (denoted base), 2) the “25% C+O core” model, and 3) the “50% C+O core” model.

Figure 2 displays the composition structure of the three models at day 35. While C II lines become more prominent they do not necessarily rule out the existence of a substantial amount of C+O. Fig. 1 shows that as the C+O abundance is increased, the C II $`\lambda 4619`$ line becomes stronger, as does the C II $`\lambda 6580`$ line, and the C II $`\lambda 7234`$ line. In the UV there is no evidence for, e.g., C IV. In the “25 % C+O core” model, C+O mixing has not reached its peak until $`\tau _{\text{std}}=0.76`$ (the total continuum optical depth at 5000 Å).

To verify that we are indeed sampling the line forming region with our mixed models we calculated a model 44 days after explosion (roughly 24 days after maximum light) which is shown in Figure 3. In this case, $`\tau _{\text{std}}=1`$ at 1100, 1450, 2000 km s-1 in the base model, the 25% C+O model, and the 50% C+O core model respectively, since O serves as a source of free electrons. At this epoch the same features become too strong, although since our models are pure LTE and they use only the Kurucz (1993) line list which does not include forbidden lines. Therefore, lines of \[Fe III\] and \[Co III\] which are generally assumed to produce the features near 5500 Å and 6500 Å (Bowers et al., 1997) are not included in our calculations. It is interesting to note that the velocity of the absorption minimum of the C II $`\lambda 6580`$ line in our model spectra is at about $`6000`$ km s-1, which is where C II is both prevalent, in our 50% C+O core model, and is the top of the canonical line forming region $`\tau _{\text{std}}=0.1`$. In the 25% C+O core model, the minimum moves somewhat redward to about $`5400`$ km s-1. However in this model there is no carbon between 4000 and 15000 km s-1 so the minimum is likely produced by a complicated blend of lines as is almost always the case in detailed models. The peak is almost certainly produced by C II $`\lambda 6580`$.

The models of Gamezo et al. (2003) have both C+O mixing to low velocity, and also mixing of burnt material (possibly iron) to high velocity and we have not attempted to simulate that here. In addition the results of Gamezo et al. (2003) end at 1.9 s, before the explosion has become homologous, so it is difficult to translate their plot of compositions as a function of radial coordinate into a model in homologous expansion where both velocity and radius are equally good radial coordinates.

## 4 Discussion

Our results clearly demonstrate that one cannot simply examine the radial composition structure of a model and determine whether a hydrodynamical model is viable or not. On the other hand, we do claim to have shown that the model of Gamezo et al. (2003) *is* a viable model for SNe Ia. In particular, it is not possible to determine the exact composition structure from their paper; future work should attempt to make some educated guesses about the composition structure and calculate synthetic spectra directly from their (angle averaged) model. Sorokina & Blinnikov (2003) calculated multi-band light curves from the model of Reinecke et al. (2002b) and found that having burnt products at high velocity was somewhat helpful particularly in $`U`$ and $`B`$. This is an indication that mixed models may be viable for SNe Ia, whether produced by deflagration or some other mechanism. Nevertheless, a full synthetic spectrum calculation is required before firm conclusions can be drawn about the viability or lack thereof of any particular hydrodynamical model. What we have shown is that unburnt C+O in the central regions are not necessarily ruled out at the epochs we study.

One possible problem with the model of Gamezo et al. (2003) is that late time nebular spectra would show strong forbidden oxygen lines, which have not typically been observed in SNe Ia, although Minkowski (1939) did identify weak narrow lines in SN 1937C of \[O I\] $`\lambda 6300`$ and $`\lambda 6394`$ beginning 184 d after maximum. The lines were absent by 275 d after maximum. Again, until detailed late-time synthetic spectra have been calculated, it is impossible to draw a firm conclusion. It may be that since the material in the center of the Gamezo et al. (2003) model is dense and the radioactive 56Ni is well mixed with the C+O that the oxygen stays ionized for a long time, but it does seem likely that \[O I\] lines should appear at some point.

In addition, the models of Gamezo et al. (2003) mix burnt material into the outer layers as well as mixing unburnt material into the inner layers. We have not accounted for that effect in these simple exploratory calculations. However, we know that SNe Ia show a range of “temperatures”, forming a “temperature sequence” that correlates with the light curve shape sequence (Nugent et al., 1995b, S. Bongard et al., in preparation). Exactly how the sequence is produced is not understood, although it is clear that models that produce more nickel will be both hotter and have higher opacity than those with less nickel (Khokhlov et al., 1993). “Hot” SNe Ia show lines of Fe III (Fisher et al., 1999; Hatano et al., 2002) and Si III (D. Branch et al., in preparation), but in parameterized models these lines are produced via extra thermalization (Nugent et al., 1995b) and not by direct excitation due to non-thermal electrons produced by 56Ni decay. If very inhomogeneous mixing of burnt and unburnt material in the outer layers is a general feature of deflagration models (Khokhlov, 2001), then the observed uniformity of the Si II 6150Å feature probably rules these models out (Thomas et al., 2002).

From our simple parameterized results, we conclude that the viability or non-viability of any hydrodynamical model cannot simply be determined by the presence of a particular element in a certain velocity interval. The converse is not true, the absence of a particular element in a velocity interval where it is observed is fatal to any particular hydrodynamical model. It will be helpful when future hydrodynamical simulations can extend to the regime of homologous expansion and when the full compositions structure of the models are available.

We thank David Branch for very useful discussions and significantly improving the presentation of this work. This work was supported in part by NASA grant NAG5-12127, NSF grant AST-0204771, and an IBM SUR grant to the University of Oklahoma; and by NSF grant AST-9720704, NASA ATP grant NAG 5-8425 and LTSA grant NAG 5-3619 to the University of Georgia. PHH was supported in part by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. Some of the calculations presented in this paper were performed at the San Diego Supercomputer Center (SDSC), supported by the NSF, on the IBM pSeries 690 of the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN), and at the National Energy Research Supercomputer Center (NERSC), supported by the U.S. DOE. We thank all of these institutions for a generous allocation of computer time.


text:
# Biology helps you to win a game

## I Introduction

Is the “survival of the fittest” principle limited to biology only? Perhaps not and there could be other spheres of life in which this principle is applicable. Competition plays a key role and in order to compete and thus survive in any environment or situation, one primarily needs to adapt in order to succeed. Then what is adaptation and evolution? Adaptation is an alteration or adjustment in structure or habits, often hereditary, by which a species or individual improves its condition in relationship to its environment. Evolution is the change in the genetic composition of a population during successive generations, as a result of natural selection acting on the genetic variation among individuals, and resulting in the development of a new species. Here, we show that in the behaviour of various complex systems found in natural and social environments parisi ; huberman ; nowak ; lux ; arthur , that can be characterized by the competition among interacting agents for scarce resources, adaptation to the environment plays a very important role.

These agents could be diverse in form and in capability, ranging for example, from carcinogenic cells in the human body to multinational firms in the global financial market. In these dynamically evolving complex systems the nature of agents and their behaviour differ a lot but they have a common underlying mechanism. In order to have a deeper understanding of the interactions of the large number of agents, one should first consider the individual capabilities of the agents. Its behaviour may be thought of as a collection of simple rules governing “responses” to numerous “stimuli”. The rules of action serve as the agents’ strategies, and the behaviour of an agent is the rules acting sequentially. Therefore, in order to model any complex dynamically adaptive system, a major concern is the selection and representation of the stimuli and responses, since the behaviour and strategies of the component agents are determined thereby. Then the agent needs to adapt to different situations, where the experience of an agent guides it to change its structure so that as time passes, the agent learns to make better use of the environment for its own benefit. However, the timescales over which the agents adapt vary from one individual to another and also from one system to another.

In complex adaptive systems, many interesting temporal patterns are produced, since a major part of the environment of a particular agent includes other adaptive agents and a considerable amount of agent’s effort goes in adaptation and reaction to the other agents. Thus the situation is considerably different and more complicated than in game theory game  and conventional theories in economics, where the study is of patterns in behavioural equilibrium that induce no further interaction.

In this paper, we study a simple game based on the basic minority game challet1 ; challet2 ; cavagna ; riolo ; lamper , where the agents adapt themselves by modifying their strategies from time to time, depending on their current performances, using genetic crossover mechanisms holland ; goldberg ; lawrence ; Marko1 . The game can be a very simple representation of a complex adaptive system. We make a comparative study of their performances with the various mechanisms and in a “test” situation.

## II Model

In this section we give a brief description of the model. The basic minority game consists of an odd number $`N`$ of agents who can perform at a given time $`t`$, any of the two possible actions denoted here by $`0`$ or $`1`$. The minority game was based on the El Farol bar problem, created by Brian Arthur, in which a population of agents have to decide whether to go to the bar every Thursday night, and so there were two possible actions “to attend” denoted by $`1`$ and “not to attend” denoted by $`0`$, depending on whether the bar was too crowded or not arthur . An agent wins the game if it is one of the members of the minority group. All the agents are assumed to have access to finite amount of “global” information: a common bit-string “memory” of the $`M`$ most recent outcomes. With this there are $`2^M`$ possible “history” bit-strings. Now, a “strategy” consists of two possible responses, which in the binary sense are an action $`0`$ or action $`1`$ to each possible history bit-strings. Thus, there are $`2^{2^M}`$ possible strategies constituting the whole “strategy space”.

Each time the game has been played, time $`t`$ is incremented by unity and one “virtual” point is assigned to the strategies that predicted the correct outcome and the best strategy is the one which has the highest virtual point score. The performance of a player is measured by the number of times the player wins, and the strategy, which the player uses to win, gets a “real” point. The number of agents who have chosen a particular action, say $`1`$ which represents “to attend”, is denoted by $`A_1(t)`$ (also referred as “attendance”) and it varies with time. We have plotted the attendance and performance for the basic minority game in Fig. 1.

Now we define the total utility of the system as the number of persons in the minority group at a given time $`t`$. For convenience, we mathematically define a scaled utility (total utility/maximum utility) as

$$U=[(1\theta (x_tx_M))x_t+\theta (x_tx_M)(Nx_t)]/x_M,$$
(1)
where $`x_M=(N1)/2`$, $`x_t`$ is either equal to $`A_1(t)`$ or $`A_0(t)`$, and $`\mathrm{\Theta }(x)`$ is Heaviside step function:

$$\theta (x)=\{\begin{array}{cc}0\hfill & \text{ when }x0\hfill \\ 1\hfill & \text{ when }x>0.\hfill \end{array}$$

The players examine their performances after every time interval $`\tau `$. If a player finds that he is among the fraction $`n`$ (where $`0<n<1`$) who are the worst performing players, he adapts himself and modifies his strategies. The mechanism by which the player creates new strategies is genetic crossover, whereby he selects the two “parents” from his pool of $`k`$ strategies and creates two new “children” lawrence ; Marko1 , as described in Fig. 2.

If the parents are chosen randomly from the pool of strategies then the mechanism represents a “one-point genetic crossover” and if the parents are the best strategies then the mechanism represents a “hybridized genetic crossover”. The children may replace parents or two worst strategies and accordingly four different interesting cases arise: (a) one-point genetic crossover with parents “killed”, i.e. parents are replaced by the children, (b) one-point genetic crossover with parents “saved”, i.e. the two worst strategies are replaced by the children but the parents are retained, (c) hybridized genetic crossover with parents “killed” and (d) hybridized genetic crossover with parents “saved”.

It should be noted that the mechanism of evolution of strategies is considerably different from earlier attempts challet1 ; li1 ; li2 . This is because in this mechanism the strategies are changed by the agents themselves and even though the strategy space evolves continuously, its size and dimensionality remain the same.

The Hamming distance $`d_H`$ between two bit-strings is defined as the ratio of the number of uncommon bits to the total length of the bit strings. It is a measure of the correlation between two strategies:

$$d_H=\{\begin{array}{ccc}0\hfill & \text{ correlated}\hfill & \\ 0.5\hfill & \text{ uncorrelated}\hfill & \\ 1\hfill & \text{ anti-correlated}\hfill & \end{array}$$

which can be plotted as the game evolves.

## III Results

In order to determine which mechanism is the most efficient, we have made a comparative study of the four cases, mentioned earlier. We plot the attendance as a function of time for the different mechanisms in Fig. 3.

In Fig. 4 we show the total utility of the system in each of the cases (a)-(d), where we have plotted results of the average over 100 runs and each point in the utility curve represents a time average taken over a bin of length 50 time-steps. The simulation time is doubled from those in Fig. 3, in order to expose the asymptotic behaviour better. On the basis of Figs. 3 and 4, we find that the case (d) is the most efficient.

In Fig. 5 (a) one can see the evolution of the average Hamming distance of all the strategies of a player in a game, where the player adapts using one-point genetic crossover and the two worst strategies are replaced by the children and the parents are also saved. It should be noted that the Hamming distance can change only when the worst strategies are replaced by the children and the parents are saved, where the bits in a strategy pool can change over time. Otherwise the bits in the pool of strategies remain the same. We observe that the curves tend to move downwards from around $`0.5`$ towards zero, which means that as the time evolves, the correlation amongst the strategies increases and the strategies in the pool of a particular agent converges towards one strategy. The nature of the curves depend a lot on the parameters of the game. In Fig. 5 (b) one can see the evolution of the average Hamming distance of all the strategies of a player in the game, where the player adapts using hybridized genetic crossover and the two worst strategies are replaced by the children and the parents are also saved. Here too, the strategies in the pool of a particular agent converges towards one strategy, and at a faster rate than with the previous mechanism. We observe that increasing memory $`M`$ does not change dramatically the convergence rate, but as we increase the number of strategies in the pools, the convergence slows down.

In order to investigate what happens in the level of an individual agent, we created a competitive surrounding– “test” situation where after $`T=3120`$ time-steps, six players begin to adapt and modify their strategies such that three are using hybridized genetic crossover mechanism and the other three one point genetic crossover, where children replace the parents. The rest of the players play the basic minority game. In this case it turns out that in the end the best players are those who use the hybridized mechanism, second best are those using the one-point mechanism, and the bad players those who do not adapt at all. In addition it turns out that the competition amongst the players who adapt using the hybridized genetic crossover mechanism is severe.

## IV Conclusion

We can summarize our findings by stating that adaptation improves not only the individual player’s performance but also improves the total utility of the system. The best results are found for the players who adapt and modify their strategies using the hybridized genetic crossover mechanism and the children replace the two worst strategies and the parents are saved. The mechanism of adaptation is very simple and can be used to model different complex adaptive systems. It can also be potentially developed to include other features like mutation. We can thus say that in a way, “biology helps you to win a game”.

###### Acknowledgements.

This research was partially supported by the Academy of Finland, Research Centre for Computational Science and Engineering, project no. 44897 (Finnish Centre of Excellence Programme 2000-2005).