format: text

text:
# Deriving the Metallicity Distribution Function of Galactic Systems

## 1 Introduction

The metallicity distributions of stars in different environments yield important information about the age and formation history of stellar systems. Particular attention has been paid to the Metallicity Distribution Function (MDF) of G-dwarfs in the solar neighbourhood, as this provides one of the strongest constraints on Galactic Chemical Evolution models. These models solve the sets of differential equations governing the formation and destruction of the elements through stellar processes (Tinsley 1980). The over-prediction of low-metallicity stars, the so-called “G-dwarf problem”, is typically resolved by allowing the disk of the Milky Way to form gradually from accreting gas on a timescale of about 7 Gyr. Recent models (e.g. Chiappini et al. 1997; Goswami & Prantzos 2000) have successfully reproduced the observed G-dwarf MDF using dual-phase accretion models in which the halo/thick disk component evolves independently of the thin disk and on a rapid timescale.

An implicit assumption of most theoretically predicted G-dwarf metallicity distributions is that all G-dwarfs have lifetimes older than the age of the Milky Way. In fact, the stars in the Rocha-Pinto & Maciel (1996) sample have masses ranging from 0.7 M to 1.1 M, while Wyse & Gilmore’s (1995) dataset spans 0.8 M to 1.2 M. Some of these stars have lifetimes shorter than the age of the disk. Furthermore, the earliest formed G-dwarfs will have even shorter lifetimes owing to their low metallicity, which reduces their opacity and raises their luminosity. For instance, a 1.3 M star leaves the main-sequence after about 4 Gyr at solar metallicity, or after less than 3 Gyr for Z = 0.05 Z (Schaller et al. 1992). In the low-metallicity environment of the early Milky Way, a 1 M star is expected to leave the main-sequence after about 6 Gyr. The preferential loss of the oldest G-dwarfs from empirical datasets tends to bias the MDF towards higher metallicity and suppresses the low-metallicity tail. However, the protracted formation history of the Galactic disk means that stellar lifetime effects alone cannot account for the paucity of metal-poor G-dwarfs with respect to “closed box” models of chemical evolution. We have investigated the impact of these effects on the shape and peak of measured MDFs in different stellar environments, the results of which are presented in Section 4.2. Rocha-Pinto & Maciel (1997) and Bazan & Mathews (1990) should be referred to for additional commentary in this regard.

Since K-dwarfs are less massive than G-dwarfs, their lifetimes are longer than the age of the Galaxy. The K-dwarf metallicity distribution should therefore trace the “true” MDF because surveys of this stellar population should not suffer from incompleteness due to stars evolving off the main-sequence. Because of their faintness, however, accurate metallicities of large numbers of K-dwarfs have only recently been obtained. Kotoneva et al. (2002) constructed a K-dwarf distribution using a sample of 220 nearby K-dwarfs drawn from the Hipparcos catalogue (ESA 1997). They were able to select K-dwarfs by absolute magnitude and restrict the mass range to 0.7 - 0.9 M, for which evolutionary effects are unimportant.

Galaxy Evolution tool (Fenner & Gibson 2003) - GEtool111http://astronomy.swin.edu.au/GEtool/ \- was used to numerically solve the equations governing the rate of change of gas mass and elemental abundances due to star formation, supernovae feedback, and the continual accretion of gas onto the disk. The model presented here is capable of reproducing the main observed properties of the Milky Way, including the present day radial gas distribution, radial abundance gradients, age-metallicity relation, stellar metallicity distribution and elemental abundance ratios. Section 2 outlines the features of the model, while Section 3 describes the adopted stellar nucleosynthetic prescriptions. The predicted K-dwarf metallicity distribution is presented in Section 4.1, where it is compared against the measured K-dwarf MDF in order to infer the possible formation history of the Galaxy. In Section 4.2 are shown results from a study into the effect of finite stellar lifetimes on the measured MDFs of populations in different star forming environments.

## 2 The Model

### 2.1 Infall scheme

Our Galactic evolution model assumes two main formation phases. The formation of halo stars is associated with the first episode of infalling material, while the second phase leads to disk formation. To aid the comparison with other dual-phase Galactic chemical evolution studies the model is restricted to simulate *two* rather than *three* infall phases. The thick disk is assumed to correspond to a merger heated thin disk and is therefore chiefly associated with the early stages of the second formation epoch. The difficulty inherent in using this kind of dual-phase model to describe a three component system is further discussed in Section 4.1. Simulating the evolution of a separate thick disk component will be explored in forthcoming work.

The initial phase of star formation occurs on a timescale of $`0.5`$ Gyr and enriches the initially primordial gas to a metallicity of \[Fe/H\]$``$$``$1.2. The second formation phase is delayed by 1 Gyr with respect to the first phase and has a more protracted star formation history; High-Velocity Clouds (HVCs) may represent the present-day source of this Galactic star formation fuel. Drawing on observations of the chemical composition of these clouds (Wakker et al. 1999; Gibson et al. 2001; Sembach et al. 2002), the second accretion phase assumes the gas is slightly metal-enriched (Z = 0.1 Z) with an enhancement of $`\alpha `$-elements relative to iron. The adopted value of \[$`\alpha `$/Fe\] = $`+`$0.4 is consistent with $`\alpha `$/Fe ratios found in metal-poor stars (e.g. Ryan, Norris & Beers, 1996). The precise dependence of key observational constraints on the assumed chemical composition of the infalling gas will be presented in a forthcoming paper (Fenner & Gibson 2003).

After Chiappini et al. (1997), we assume that the rate at which material is accreted during these phases declines exponentially. The evolution of total surface mass density $`\sigma _{tot}(r,t)`$ is given by

$$\frac{d\sigma _{tot}(r,t)}{dt}=A(r)e^{t/\tau _H(r)}+B(r)e^{(tt_{delay})/\tau _D(r)}$$
(1)
where the infall rate coefficients $`A(r)`$ and $`B(r)`$ are chosen in order to reproduce the present-day surface mass density of the halo and disk components, which we take to be 10 and 45 M pc-2, respectively. These values are comparable with 17 and 54 M pc-2 used by Chiappini et al. (2001). The adopted timescales for the infall phases are $`\tau _H`$ = 0.5 Gyr and $`\tau _D`$ = 7.0 Gyr at the solar radius $`r_{}`$ = 8 kpc. The “inside-out” functional form for $`\tau _D(r)`$ (Romano et al. 2000) is adopted, with the Milky Way age taken to be 13 Gyr.

### 2.2 Star formation formalism

We adopt a star formation prescription which is based upon the hypothesis that star formation is triggered by the compression of interstellar material from spiral arm motion (e.g. Prantzos & Silk 1998). This implies that the star formation rate (SFR) is proportional to the angular frequency of the spiral pattern and therefore inversely proportional to radius. The star formation law incorporated in this model is given by

$$\psi (r,t)=\nu \sigma _{gas}^2(r,t)\left(\frac{r_{}}{r}\right)$$
(2)
where the value of the efficiency factor $`\nu `$ is constrained by the present day gas fraction.

### 2.3 Initial mass function

The shape of the stellar initial mass function (IMF) influences the quantity of Galactic material locked up in stars of different masses, which in turn determines the rate at which different elements are released into the interstellar medium. The models presented in Section 4.1 use the Kroupa et al. (1993) three-component IMF. With respect to the Salpeter (1955) and Scalo (1986) IMFs, this function has a steeper slope for high masses, resulting in less material being processed by massive stars. We impose an upper mass limit of 60 M on stellar formation in order to recover the observed trend of \[O/Fe\] at low metallicity (e.g. Carretta et al. 2000, Melendez et al. 2001) while avoiding the overproduction of oxygen.

### 2.4 Properties of the halo, thick disk, and thin disk components

The ratio of the local number density of thick disk to thin disk stars lies between 1% and 6% (Gilmore et al. 1995). Both disk components are well-fitted in vertical height, $`z`$, with exponential functions. Taking the thick and thin disk scale heights to be 1 kpc and 330 pc, respectively (Gilmore et al. 1995), the ratio of the thick to thin disk surface densities lies between 2% and 20%. The local stellar density of halo stars relative to thin disk stars is only $`0.1`$% with a mean metallicity of \[Fe/H\] = $``$1.7 (Norris & Ryan 1991), while thick disk stars have typical metallicities in the interval $``$0.4 $`<`$ \[Fe/H\] $`<`$ $``$1.2 with a mean of $``$0.7 (Gilmore et al. 1995). In comparison, the MDF of thin disk stars peaks near \[Fe/H\] = $``$0.2 (McWilliam 1997). In addition to metallicity, kinematic behaviour is used to distinguish between these three components. Thick disk kinematics lie somewhere between those of the halo and the thin disk, and describe a distinct population.

## 3 Yields

The observed behaviour of elemental ratios such as \[O/Fe\] with \[Fe/H\] are important observational constraints on chemical evolution models. Oxygen is chiefly produced in massive stars, whereas a significant fraction of iron is supplied by Type Ia supernovae involving binary systems of lower mass stars. The characteristic timescale on which oxygen is ejected into the ISM is therefore much shorter than that of iron, owing to the mass dependence of main-sequence lifetimes.

In order to recover the observed behaviour of elemental abundance ratios, the instantaneous recycling approximation has been relaxed. Mass and metallicity-dependent lifetimes were taken from Schaller et al. (1992).

### 3.1 Massive stars

The present work incorporates an updated set of stellar yields from Limongi et al. (2000; 2002), supplemented with a finer mass coverage grid (kindly provided by Marco Limongi), covering a range of stellar initial masses and metallicities (13$`<`$$`m/`$M$`{}_{}{}^{}<`$80 and Z/Z=0, 10-3, 1). Stellar yields are one of the most important ingredients in Galactic chemical evolution models, yet questions remain concerning the composition of ejected material, owing to the uncertain role played by processes such as mass loss, rotation, fall-back, and the location of the mass cut in supernovae, which separates the remnant from the ejected material.

The amount of iron released in supernova explosions of massive stars is particularly sensitive to fall-back and the location of the mass cut. There are few supernova observations with which to directly infer the amount of iron ejected. We have therefore taken the opposite approach of using the observed trend of abundance ratios to indirectly constrain iron yields. The Limongi et al. (2000; 2002) yields are presented for a range of mass cuts, allowing a self-consistent treatment of all yields (c.f. Goswami & Prantzos 2000; Argast et al. 2002). We have found that the trend of \[O/Fe\] with \[Fe/H\] is best recovered by suppressing the ejection of iron in stars more massive than about 30 M. This leads to a smaller contribution to the interstellar iron abundance from Type II supernovae than is obtained from the canonical Woosley & Weaver (1995) model. A preliminary relationship between the mass cut and iron yield was presented in Fenner, Gibson & Limongi (2002); a full accounting will be provided in Fenner & Gibson (2003).

### 3.2 Low- and intermediate-mass stars

The metallicity-dependent yields of Renzini & Voli (1981) were used for stars in the 1 $``$ m/M$``$8 range. Yields for stars with masses between 8 M and 13 M were estimated by interpolating between the highest mass in Renzini & Voli and the lowest mass in Limongi et al (2000; 2002). It should be noted, however, that it is not yet clear whether stars in the range 8 $`<`$ m/M$`<`$ 13 contribute to interstellar enrichment, or whether the processed material becomes trapped in the remnant following core-collapse.

### 3.3 Type Ia supernovae

The yields for Type Ia SNe were taken from the W7 model of Thielemann et al. (1993). The fraction of binaries resulting in Type Ia supernovae was fixed by requiring that the model reproduce the local age-metallicity relation; a fraction of 3% was adopted, leading to model predictions consistent with observation. Type Ia supernovae supply about one half of the interstellar iron abundance in this model.

## 4 Results

### 4.1 The solar neighbourhood

The predicted star formation rate in the solar vicinity (*solid line*) is shown as a function of time in the upper panel of Figure 1 along with the range of measured values (*vertical bar*). The star formation rate in the halo peaks at almost 3 M pc-2 Gyr-1 after 500 Myr. The predicted present day rate of star formation in the thin disk of 2.7 M pc-2 Gyr-1 agrees with observations (Chiappini et al. 2001) and is about 2.5 times lower than the maximum rate of star formation which occurred $``$ 10 Gyr ago, 2-3 Gyr after the onset of the second infall phase.

The dotted line in the upper panel of Figure 1 depicts the star formation rate given by the Chiappini et al. (2001) dual-phase Galactic chemical evolution model. Like the Milky Way model presented in this paper, the Chiappini et al. model is consistent with many observational constraints, including the Kotoneva et al. (2002) empirical MDF. Despite both models satisfying key empirical constraints such as the K-dwarf metallicity distribution, they do predict very different properties for the halo, thick disk, and thin disk components of the Milky Way. This is largely a consequence of different star formation rates in the earliest epoch of Galactic evolution. The predicted number of stars born during the first billion years differs by a factor of 2.5 between the two models, whereas the *total* number of stars formed differs by less than 20%. The higher initial rate of star formation in Chiappini et al. leads to a more rapid enrichment of the interstellar medium, when compared against the model presented here. The implications of these differences in terms of the Milky Way’s stellar components are quantified further below.

First, the model we have employed here is successful in reproducing the new Kotoneva et al. (2002) K-dwarf MDF, as shown in the lower panel of Figure 1. The predicted MDF has been convolved with a $`\sigma `$ = 0.1 dex Gaussian in \[Fe/H\], consistent with the observational uncertainties in the data. Plotted against the model results (*solid curve*) is the K-dwarf MDF from Kotoneva et al. (*squares with error bars*) and the G-dwarf MDF from Hou et al. (1998) (*dotted line histogram*). The model MDF provides a satisfactory agreement with the Kotoneva et al. dataset, although the model peaks $``$ 0.2 dex higher in \[Fe/H\]. The model also follows closely the Hou et al. G-dwarf MDF, which exhibits a longer metal-poor tail than the K-dwarf distribution. While most authors measure the peak of the local G-dwarf metallicity distribution to lie between -0.2 and -0.3 dex (e.g. Wyse & Gilmore 1995; Rocha-Pinto & Maciel 1996), there is still not complete agreement. Haywood (2001) found that the MDF of a color-selected sample was centered on solar metallicity, however Kotoneva et al. showed that a bias toward metal-rich stars might arise due to a systematic trend between metallicity distribution and stellar color.

The complete absence of stars with \[Fe/H\] $`<`$$``$1 in the Kotoneva et al. (2002) sample is not reproduced by our model (*solid curve*). In Chiappini et al. (1997; 2001) and Kotoneva et al. (2002), any low-metallicity tail in the theoretical MDF is suppressed by neglecting either: 1) stars with \[Fe/H\] $`<`$$``$1.2 (e.g. Chiappini et al. 1997; Kotoneva et al. 2002); or 2) stars born before t = 1 Gyr (which effectively ignores stars with \[Fe/H\] $`<`$$``$0.6) (e.g. Chiappini et al. 2001). In comparison, we have included all long-lived stars in the metallicity distribution function. Stars with \[Fe/H\] $`<`$$``$1.2 represent a very small fraction of present-day stars in the Chiappini et al. model due to their strong initial burst of star formation. Efficient, early, star formation leads to a metallicity \[Fe/H\] $`=`$$``$1.2 being reached on a very short timescale ($`<`$ 0.1 Gyr). Such rapid initial enrichment is one of the main differences between the present model and that of Chiappini et al. We assume a more moderate initial phase of star formation, leading to a metal-poor tail (i.e. \[Fe/H\] $`<`$$``$1.2) in the MDF which contains $``$ 4% of the stars. The age associated with this metal-poor tail is 12.5 Gyr, while the median age of our predicted thin disk dwarfs is 7.5 Gyr. This is in excellent agreement with Hansen et al. (2002) who derive ages of 7.3 $`\pm `$ 1.5 Gyr and 12.7 $`\pm `$ 0.7 Gyr for Galactic thin disk white dwarfs and for the halo globular cluster M4, respectively.

The age-metallicity relation (AMR) for our dual-infall model is shown in Figure 2 (*solid curve*). It takes $``$1 Gyr to enrich the interstellar medium in the solar neighborhood to a metallicity of \[Fe/H\] = $``$1.25. In contrast, the Chiappini et al. (2001) model (*dotted curve*) reaches the much greater metallicity of \[Fe/H\] $``$$``$0.5 over the same time interval. The AMR from Chiappini et al. (2001) differs from that in the 1997 version of their two-infall model, owing to different timescales for the first episode of infall. Their halo/thick disk accumulates on a timescale of 2 Gyr in the 1997 model, with a 2 Gyr delay before the thin disk starts to form. In the later model, these timescales are shortened to $``$1 Gyr, giving rise to an even steeper AMR than in the 1997 paper. The kink seen in the Chiappini et al. AMR is due to the metallicity of the ISM becoming diluted when the second phase of infall commences at 1 Gyr. It is not until the thin disk gas surface density reaches a threshold of 7 M pc-2 that Chiappini et al. allow thin disk stars to form. No such star formation threshold is incorporated in our model.

We wish to stress that caution must be employed when drawing conclusions regarding the nature of three stellar population components of the Milky Way (thin disk, thick disk, halo), when the models are inherently limited by their dual-infall nature. That said, if one adopts a chemical criterion of $``$1.0 $`<`$ \[Fe/H\] $`<`$$``$0.5 for the definition of “thick disk”, then the age of the thick disk predicted under the Chiappini et al. formalism would be centered upon 12.5 Gyr, with a vanishingly small age-spread ($`<`$ 0.2 Gyr), and possess kinematics closely aligned with the first of the infall phases (i.e. the halo). Conversely, the thick disk predicted under our formalism would be 2-3 Gyr younger, with an age-spread an order-of-magnitude larger ($``$ 2 Gyr), and possess kinematics more closely resembling those of the second infall phase (i.e. the disk). The Kotoneva et al. (2002) sample of solar neighborhood K-dwarfs identified only 2 of 431 stars as having halo kinematics. Thus halo stars are rare, whereas the Chiappini model would associate a large fraction of stellar mass with the halo component.

In the Chiappini et al. model of Galaxy formation, the first infall phase is associated with the formation of the halo and thick disk. Thus these components are coeval and form prior to, and independently of the thin disk. In such a scenario, the thick disk would have little dispersion in age and the stars would exhibit primarily halo kinematics. The model presented here associates the first infall phase with the build up of the halo. In contrast with Chiappini et al., the formation of the thick disk overlaps with the early evolution of the thin disk (during the second infall phase). Our thick disk, as defined by metallicity, would have predominantly disk-like kinematics and a broader range of ages.

While both formalisms clearly have problems in describing a three-component system, we favour the scenario presented in this paper, in which the intermediate metallicity population (thick disk) has a typical age several Gyrs younger than the halo (as inferred by observation - e.g. Chayboyer et al. 2000; Ibukiyama & Arimoto 2002) and kinematics akin to that of disk-like population, and not a rapid-collapse halo population. In that sense, our model is more consistent with the favored model of Galactic thick disk formation whereby an existing stellar thin disk “puffs up” during the process of a merger (e.g. Wyse 2001).

### 4.2 Lifetime effects and MDFs in external galaxies: A note of caution

Recent progress has been made in measuring the metallicity distribution in the stellar halos of external galaxies (e.g. Durrell et al. 2001; Harris & Harris 2000, 2002). The observed metallicity distribution of these systems, which presumably formed early and on a short timescale, is likely to be biased by the effect of finite stellar lifetimes. In an attempt to quantify the discrepancy between the observed and “true” MDF due to lifetime effects, we constructed a “toy” model of the NGC 5128 halo using the chemical evolution code described in this paper. This NGC 5128 halo model has some key differences with respect to the Milky Way model described in this paper. The main differences are: a single phase of star formation; a shorter formation timescale ($`1`$ Gyr); and initial composition of gas is primordial.

In stellar systems such as halos, a significant fraction of stars classified as G-dwarfs may have evolved off the main-sequence. This has the effect of preferentially removing older and therefore more metal-poor stars from the observed present-day G-dwarf population. In a similar way, the metal-poor tail in the measured G-dwarf MDF is further suppressed by the shorter lifetimes of low mass, metal-poor stars (e.g. Schaller et al 1992). On the other hand, the disk of our own Galaxy has a more extended star formation history, with the lengthier timescale and the higher proportion of young to old stars acting to mitigate the “loss” of G-dwarfs from the present-day MDF.

In Figure 3, the star formation history of the Milky Way disk is compared with a more burst-dominated star formation history, which we here assume to be representative of ellipticals (or pure halo systems). The bulk of the stars in this burst model are born during the first 1.5 Gyr, after which the gas supply is depleted and the star formation rate declines rapidly. Conversely, it takes $`>`$ 5 Gyr for the majority of stars to form in the disk model.

Figure 4 shows the predicted MDF of truly long-lived stars (*black lines*) and G-dwarfs (defined as having masses in the range 0.8 $`<`$ m/M$`<`$ 1.4) (*gray lines*) for the halo (*dashed lines*) and Milky Way (*solid lines*) models. The observed MDF of the NGC 5128 outer halo (Harris & Harris 2002) is plotted as a histogram, where we have used the relation \[Fe/H\] = \[m/H\] $``$ 0.25 to transform the published values of \[m/H\] into units of \[Fe/H\]. The black lines show the MDF that is expected if all G-dwarfs had lifetimes longer than the age of the galaxy (i.e., masses m$`<`$0.8 M). These curves have been normalised to contain the same area as the histogram. The gray curves show the MDF that would be measured if the mass range of stars in the sample is 0.8 $`<`$ m/M$`<`$ 1.4. These have been normalised to reflect the loss of stars in each metallicity bin owing to evolution off the main-sequence. The dotted line also shows the finite lifetime halo model MDF, scaled to illustrate the fit to the data. The loss of stars from a present-day sample of G-dwarfs can be seen to narrow the MDF and shift the distribution toward higher metallicities. While evident in both models, this behaviour is more pronounced in environments such as halos, where the bulk of the stars presumably formed $`>`$10 Gyr ago.

By ignoring lifetime effects, one can overestimate the true mean metallicity of long-lived stars. To gauge the magnitude of this effect, we plot the cumulative metallicity distribution in Figure 5. For this particular halo model, mass-dependent stellar lifetimes serve to shift the median metallicity of the population by $`\mathrm{\Delta }`$\[Fe/H\]=$`+`$0.3.

## 5 Conclusions

We have reproduced the metallicity distribution of nearby K-dwarfs using a model in which the Milky Way formed during two accretion episodes. Several key differences between this and other dual-phase infall models are (i) the assumption that the thin disk formed from metal-enriched $`\alpha `$-enhanced extragalactic material, (ii) the adoption of a new comprehensive set of stellar yields from Limongi et al. (2000; 2002) for calculating the contribution to the enrichment of the ISM from massive stars, and (iii) a more gentle initial phase of (halo) star formation, resulting in a shallower (but still observationally consistent) age-metallicity relation.

Taking the stellar halo as having \[Fe/H\] $`<`$$``$1.2 we predict that this component has an age of 12-13 Gyr and formed on a timescale of $``$1 Gyr. If stars with metallicity $``$1.0 $`<`$ \[Fe/H\] $`<`$$``$0.5 are assigned to the thick disk, then the median age of this population is predicted to be $``$10 Gyr with a greater scatter in ages ($``$ $`\pm `$2 Gyrs) than for the halo. The age difference of $`5`$ Gyr between the halo and thin disk in our model is consistent with empirical estimates based upon white dwarf luminosity functions (Hansen et al. 2002).

The model presented in this paper was compared with a similar dual-phase infall model from Chiappini et al. (2001) that also reproduces the observed K-dwarf metallicity distribution. We showed that even though differences in the rate of star formation in the early Galaxy might not be apparent in present-day observables such as the MDF, they do lead to different predictions for the ages and kinematics of the halo, thick, and thin disk components. This degeneracy is inherent in dual-infall models when applied to the modeling of three-component systems.

## Acknowledgments

It is with pleasure that we thank the referees Sean Ryan and Eira Kotoneva for their very valuable comments on the manuscript.

We wish to thank Chris Flynn for many invaluable discussions. We are also grateful to Chris, Bill Harris, and Gretchen Harris, for providing their observational metallicity distribution functions prior to publication.

BKG acknowledges the financial support of the Australian Research Council through its Large Research Grant Program (Grant ID #A00105171).

## References

Argast, D., Samland, M., Thielemann, F. -K. & Gerhard, O. E. 2002, A&A, 388, 842
Bazan, G. & Mathews, G.J. 1990, ApJ, 354, 644
Carretta, E., Gratton, R. G. & Sneden, C. 2000, A&A, 356, 238
Chaboyer B., Sarajedini, A. & Armandroff, T.E. 2000, AJ, 120, 3102
Chiappini, C., Matteucci, F. & Gratton, R. 1997, ApJ, 477, 765
Chiappini, C., Matteucci, F. & Romano, D. 2001, ApJ, 554, 1044
Durrell, P.R., Harris, W.E. & Pritchet, C.J. 2001, AJ, 121, 2557
ESA, 1997, The Hipparcos and Tycho Catalogues, ESA, SP-1200
Fenner, Y. & Gibson, B.K. 2003, in preparation
Fenner, Y., Gibson, B.K. & Limongi, M. 2002, Ap&SS, 281, 537
Gibson, B.K., Giroux, M.L., Penton, S.V., Stocke, J.T., Shull, J.M. & Tumlinson, J. 2001, AJ, 547, 3280
Gilmore, G., Wyse, R.F.G. & Jones, J.B. 1995, AJ, 109, 1095
Goswami, A. & Prantzos, N. 2000, A&A, 359, 191
Hansen, B.M.S., Brewer, J., Fahlman, G.G., Gibson, B.K., Ibata, R., Limongi, M., Rich, R.M., Richer, H.B., Shara, M.M. & Stetson, P.B. 2002, ApJ, 574, L155
Ibukiyama, A. & Arimoto, N. 2002, A&A, 394, 927
Harris, G.L.H. & Harris, W.E. 2000, AJ, 120, 2423
Harris, W.E. & Harris, G.L.H. 2002, AJ, 123, 3108
Haywood, M. 2001, MNRAS, 325, 1365
Hou, J.L., Chang, R. & Fu, C. 1998, in Pacific Rim Conference on Stellar Astrophysics, ed. K.L. Chan, K.S. Cheng & H.P. Singh (San Francisco: ASP), p. 143
Kotoneva, E., Flynn, C., Chiappini, C. & Matteucci, F. 2002, MNRAS, 336, 879
Kroupa, P., Tout, C.A. & Gilmore, G. 1993, MNRAS, 262, 545
Limongi, M., Straniero, O. & Chieffi, A., 2000, ApJS, 129, 625
Limongi, M. & Chieffi, A., 2002 PASA, 19, 246
Melendez, J., Barbuy, B. & Spite, F. 2001, ApJ, 556, 858
McWilliam, A. 1997, ARAA, 35, 503
Norris, J.E. & Ryan, S.G., 1991, ApJ, 380, 403
Prantzos, N. & Silk, J. 1998, ApJ, 507, 229
Renzini, A. & Voli, M. 1981, A&A, 94, 175
Rocha-Pinto, H.-J. & Maciel, W.J. 1996, MNRAS, 279, 447
Rocha-Pinto, H.-J. & Maciel, W.J. 1997, A&A, 325, 523
Romano, D., Matteucci, F., Salucci, P. & Chiappini, C. 2000, ApJ, 539, 235
Ryan, S. G. Norris, J. E. & Beers, T. C. 1996, ApJ, 471, 254
Salpeter, E.E. 1955, ApJ, 121, 161
Scalo, J.M. 1986, Fund. Cosm. Phys., 11, 1
Schaller, G., Schaerer, D., Meynet, G. & Maeder, A. 1992, A&AS, 96, 269
Sembach, K.R., Gibson, B.K., Fenner, Y. & Putman, M.E. 2002, ApJ, 572, 178
Thielemann, F.-K., Nomoto, K., Hashimoto, M., 1993, in Prantzos, N., Vangiono-Flam, E., Casse, M., eds, Origin and Evolution of the Elements, Cambridge Univ. Press, Cambridge, p. 297
Tinsley, B.M. 1980, Fund. Cosm. Phys., 5, 287
Wakker, B.P., Howk, J.C., Savage, B.D., van Woerden, H., Tufte, S.L., Schwarz, U.J., Benjamin, R., Reynolds, R. J., Peletier, R. F. & Kalberla, P. M. W. 1999, Nature, 402, 388
Woosley, S. E. & Weaver, T. A. 1995, ApJS, 101, 181
Wyse, R.F.G. 2001, in ASP Conf. Ser. 230, Galaxy Disks and Disk Galaxies, eds. J. G. Funes & E. M. Corsini (San Francisco: ASP), 71
Wyse, R.F.G. & Gilmore, G. 1995, AJ, 110, 2771


text:
# Quantum signatures in laser-driven relativistic multiple-scattering

## Abstract

The dynamics of an electronic Dirac wave packet evolving under the influence of an ultra-intense laser pulse and an ensemble of highly charged ions is investigated numerically. Special emphasis is placed on the evolution of quantum signatures from single to multiple scattering events. We quantify the occurrence of quantum relativistic interference fringes in various situations and stress their significance in multiple-particle systems, even in the relativistic range of laser-matter interaction.

The interplay of the strongest forces in atomic physics via ultra intense laser pulses and highly charged ions is governed rather well by quantum relativistic Dirac dynamics reviewions ; reviewfields ; Maquet and Grobe (2002). On one hand, for single particles quantum effects such as tunneling, spin effects and quantum interferences have shown to be rather crucial even in the regime of ultra-short and highly relativistic dynamics Rathe et al. (1997); Panek et al. (2002). On the other hand, for many particle systems, laser-induced plasma physics was shown to be remarkably well described by classical relativistic dynamics plasmaex ; plasmathe . With the intermediate regime from few particle to cluster physics attracting increasing interest cluster , the question arises for the role of quantum effects in laser-induced relativistic dynamics.

In this letter we investigate the quantum relativistic dynamics of laser-driven multiple-scatterings of an electron being represented by a Volkov wave packet at an ensemble of highly charged ions. With an numerical accuracy, which allows for transitions even to the Dirac sea with negative energies, we quantify the interference fringes at each scattering event and the mutual interplay among those events. Clear quantum behaviour in the Dirac wave packet is identified in the highly relativistic regime after multiple scattering.

The system of interest consists of an electron which is driven by an intense laser pulse with time $`t`$ and space $`\stackrel{}{r}`$ dependent vector potential $`\stackrel{}{A}(t,\stackrel{}{r})`$ and scattered multiply at an ensemble of ions with scalar potential $`A_0(\stackrel{}{r})`$. The electronic wave packet dynamics in such an environment is characterised by the Dirac spinor $`\psi (t,\stackrel{}{r})`$ and is governed by the Dirac equation reading in atomic units as throughout the article:

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}=\left[{\displaystyle \frac{\mathrm{}c}{i}}\alpha ^j{\displaystyle \frac{}{r^j}}+\beta mc^2+q\left(A_0\alpha ^jA_j\right)\right]\psi `$ (1)
with electron charge $`q=1`$a.u., electron mass $`m=1`$ a.u. and $`\alpha ^j`$ ($`j\{1,2,3\}`$) and $`\beta `$ being the Dirac matrices Bjorken and Drell (1964). The three components of $`\stackrel{}{r}`$ and $`\stackrel{}{A}`$ are $`r^j`$ and $`A^j`$, respectively, and $`c=137.036`$ a.u. the speed of light.

Our numerical analysis takes advantage of splitting the linear Dirac Hamiltonian into a position-dependent and a derivative-dependent part. We then make use of the so-called “split-operator” technique Fleck et al. (1976), in which we propagate the wave packet successively by the position- and derivative-dependent parts and employ fast-Fourier-transformations, such that all operations are plain multiplications. With time step $`\mathrm{\Delta }t`$, the numerical error introduced this way is of the order $`\mathrm{\Delta }t^3`$ Fleck et al. (1976). For $`\mathrm{\Delta }t2\times 10^5`$ a.u.$`<\frac{\mathrm{}}{2mc^2}`$, transitions between positive- and negative-energy states are resolved, and this way we obtained convergence of our split-step propagation of $`\psi `$ even at large $`t10`$ a.u.. Further the spacing of the grid in position space needs to be suitable to resolve the maximal momenta employed. In spite of large relativistic velocities, this is not problematic, because it is the canonical rather than the kinetic momentum that has to be represented on the grid, with $`\stackrel{}{p}_{can}\frac{q}{c}\stackrel{}{A}=\stackrel{}{p}_{kin}`$. In the case of a high velocity of the particle being exclusively due to an intense laser field, $`\stackrel{}{p}_{can}`$ is even zero in polarization direction. It is non-zero in propagation direction, but its magnitude is small even for intense fields. Once scatterings with nuclei have occurred, however high canonical momenta appear, which, for the parameters used here, can be represented successfully on a grid with a spacing $`\mathrm{\Delta }x_i=0.118`$ a.u., corresponding to a maximum momentum of $`26.6`$ a.u.. The so-called fermion doubling problem, which occurs at momenta close to the highest grid momenta, is consequently also avoided Müller et al. (1998).

For the sake of reducing computing power, we introduced two advantageous techniques. Firstly, in position space, the calculation is restricted to the area centered around the rapidly moving wave packet, involving a “moving-grid” approach. Secondly, the grid size, too, is dynamically adapted in time: While a freely evolving wave packet spreads with time, a multiply scattered one does considerably more. As our simulation has to cover a substantial part of the whole laser pulse, including times where the packet is still quite small, it is possible to save considerable CPU time, noting that the time consuming two-dimensional fast-Fourier-transformations scale as $`O(N^2\mathrm{log}N)`$, where $`N\times N`$ equals the grid size. This “growing-grid” approach is also our pragmatic solution to the well-known boundary problem Alonso et al. (1997); Maquet and Grobe (2002) in Dirac calculations, at least to the point where damping functions and absorbing boundaries become unavoidable. Finally, the whole code is written to take advantage of multiple CPUs.

In a series of contour plots, we present the time evolution of an initially Gaussian shaped wave packet under the influence of a strong laser pulse, which is subsequently scattered at several highly charged ions. We use a four cycle laser pulse with amplitude $`E_0=50`$ a.u. ($`I=8.75\times 10^{19}\frac{W}{cm^2}`$) and frequency $`\omega =1`$ a.u., which features a 1.5 cycles $`sin^2`$ turn-on and turn-off, and one cycle of constant intensity in between. As amplitude and frequency suggest, we are clearly in the fully relativistic regime. The ions are modeled by static softcore potentials $`Ze^2/\sqrt{((\stackrel{}{r}\stackrel{}{r}_{\text{Ion}})^2+a)}`$, with a “Coulomb-like” small softcore parameter $`a=0.01`$ being just large enough to avoid numerical instabilities at the ions’ origins $`\stackrel{}{r}_{\mathrm{Ion}}`$. We chose their charge as a high multiple ($`Z=50`$) of the elementary charge $`e`$ in order to acquire comparable field strengths for laser and ions (at 1 a.u. distance from the ionic center).

In fig. 1, the top left graph illustrates an overview of the successive quantum relativistic scattering scenario of an electron wave packet at two highly charged ions. The initial Gaussian wave packet (positive energy and spin up only) is centered around the origin and its evolution is depicted by various snap shots along its center of mass motion (solid line) in the laser pulse. At first, after a short motion in the negative polarization and positive propagation directions during the first half cycle of the turn-on phase, the particle is visibly accelerated in the polarization direction, reaching the first upper turning point after the first whole cycle is completed. Further on, continuing with a clear Lorentz-force induced drift in the laser propagation direction, the electron wave packet is accelerated in the negative polarization direction to face its first encounter with an ionic core potential. The motion that we observe during the first unperturbed $`1\frac{1}{4}`$ cycles could be modeled rather accurately and easily with a classical Monte Carlo ensemble and is completely in agreement with known free (without nucleus) quantum wave packet results in Román et al. (2001). This includes the Lorentz-contraction of the wave packet along the direction of present velocity and its apparent rotation, i.e. precisely shearing, because of the phase differences sensed by spatially separated parts of the wave packet.

At the first encounter of the electron wave packet with the nuclei at position $`(16.5,0)`$, the laser electric field is rather small at this particular phase. Therefore relativistic Coulomb scattering dominates, involving the interference of the incoming wave with the scattered wave. The corresponding fringe structure which features a distinct maximum in forward direction followed by a series of side maxima, can be viewed in detail in fig. 1b. Further fringes from a scattering at the second ion at $`(30,0)`$ are also visible. With increasing time the wave packet evolves in the negative polarization direction towards the second lower turning point, which is reached after 1.5 cycles. The separations among the fringes have grown substantially which themselves turn continuously more parallel, with an orientation reflecting the direction of motion at the time when the scattering occurred (fig. 1c).

Finally, the wave packet, which at this stage is split in various sub-wave packets, continues in positive polarization direction. The fringes maintain their orientation when the second scattering at the second ion occurs, and are therefore joined by a set of newly created fringes with a different orientation. When two cycles are completed, the electron has reached its second upper turning point. Fig. 1d is a snapshot taken a little while before that, and one can clearly see the pattern that is generated by two sets of differently oriented interference fringes. ¿From then on, without further scatterings, the whole structure essentially remains, apart from further changes imposed by the laser field such as spreading and shearing with an essentially classical center-of-mass motion.

In figure 2, we show in a similar way the effect on the same initially Gaussian wave packet after it has passed a collection of seven ions two times. A classical estimation confirms that such a collection of heavy highly charged ions moves less than 0.04 a.u. due to Coulomb repulsion and the laser field during maximally employed interaction time of $`\mathrm{\Delta }t12`$ a.u. and may thus be assumed as resting. We note that the clear interference structures in single and double scattering are now less apparent though still visible when one moves to complex structures.

A simple analytical model is finally introduced to qualitatively confirm our numerical result. Adopting existing text-book Strange (1998) theory for three-dimensional Dirac scattering of an electron $`\psi `$ at the time-independent potential of an ion at $`\stackrel{}{r}_{\text{Ion}}=\stackrel{}{0}`$, we obtain

$$\psi (\stackrel{}{r})=\varphi (\stackrel{}{r})d^3r^{}G_0(\stackrel{}{r},\stackrel{}{r}^{},E)V(\stackrel{}{r}^{})\psi (\stackrel{}{r}^{})$$
(2)
with unperturbed Dirac wave $`\varphi (\stackrel{}{r})=w^\rho (\stackrel{}{p})e^{\frac{i}{\mathrm{}}\stackrel{}{p}\stackrel{}{r}}`$ ($`\rho \{1,2\}`$), corresponding eigenvalue $`E`$ and $`w^\rho (\stackrel{}{p})`$ being the free-electron spinor amplitude Bjorken and Drell (1964). On the right hand side of eq. (2), we replace $`\psi `$ by $`\varphi `$ (first Born approximation) and insert the relativistic free-particle Green’s function at energy $`E`$ Strange (1998)

$$G_0(\stackrel{}{r},\stackrel{}{r}^{},E)=\frac{1}{\mathrm{}c}\left[c\stackrel{}{\alpha }\underset{¯}{\overset{}{p}}+\beta mc^2+E\right]\frac{e^{\frac{i}{\mathrm{}}pR}}{4\pi R},$$
(3)
where $`R=|\stackrel{}{r}\stackrel{}{r}^{}|`$, $`p=\sqrt{\frac{E^2}{c^2}(mc)^2}`$, $`\underset{¯}{\overset{}{p}}=i\mathrm{}\stackrel{}{}`$, $`\stackrel{}{p}`$ the initial momentum, $`\stackrel{}{p}^{}=p\frac{\stackrel{}{r}}{r}`$ the final momentum and $`p=\mathrm{}k=|\stackrel{}{p}|=|\stackrel{}{p}^{}|`$ its magnitude. For the case of interest $`rr^{}`$, neglecting contributions of order $`\frac{1}{r^2}`$ and higher, and assuming a short-range potential, one finally obtains the outgoing electronic wavefunction

$`\psi (\stackrel{}{r})`$ $`=`$ $`w^\rho (\stackrel{}{p})e^{\frac{i}{\mathrm{}}\stackrel{}{p}\stackrel{}{r}}{\displaystyle \frac{1}{4\pi (\mathrm{}c)^2}}[\stackrel{}{\alpha }\stackrel{}{p}^{}+\beta mc^2+`$ (4)
$`E(\stackrel{}{p}^2)]{\displaystyle \frac{e^{\frac{i}{\mathrm{}}p|\stackrel{}{r}|}}{|\stackrel{}{r}|}}{\displaystyle }d^3r^{}V(\stackrel{}{r}^{})w^\rho (\stackrel{}{p})e^{\frac{i}{\mathrm{}}(\stackrel{}{p}\stackrel{}{p}^{})\stackrel{}{r}^{}}.`$
We are interested in the maxima of $`|\psi |^2=\psi ^{}\psi `$, or more exactly in the angles $`\vartheta _n`$ that point towards the scattering fringes. Using $`V(\stackrel{}{r}^{})=V_0\delta (\stackrel{}{r}^{})`$ becker  with $`V_0>0`$ as the simplest potential, we finally obtain, up to an additive function $`f(r)`$ and a constant pre-factor, the $`\vartheta `$-dependant part of $`|\psi |^2`$ as

$$|\psi |^2((\gamma ^21)\mathrm{cos}\vartheta +1+\gamma ^2)\mathrm{cos}\left(krkr\mathrm{cos}\vartheta \right)+f(r).$$
(5)
In the nonrelativistic case $`\gamma =\frac{E}{mc^2}1`$, the maxima of the above expression can be simplified further to read

$$\vartheta _n=\pm \mathrm{arccos}\left(1\frac{n\pi }{kr}\right),n.$$
(6)
Then to adapt the dynamics in the laser fields, one may choose for a fixed initial momentum $`\mathrm{}k`$ a distance $`r`$ where, in the *absence* of a laser field, the scattering fringes would be observed and calculate the corresponding angles $`\vartheta _n`$. Then with the laser field and, using a classical formula Salamin and Faisal (1996) and now neglecting the ionic potentials, one may propagate over a period $`t=\frac{mr}{\mathrm{}k}`$ a suitably chosen ensemble of classical particles that initially starts at the position of the scattering center with initial momenta of magnitude $`\mathrm{}k`$ in the direction of the scattering angles $`\vartheta _n`$. This simple model qualitatively confirms our numerical results while it fails to predict the final positions and separations of the fringes by better than a factor of two. In addition to the stressed approximations in the analytical approach, the mimicking of the quantum wave packet in the transition regime from scattering to free dynamics in the laser field is too delicate to compete seriously in accuracy with the up-initio quantum relativistic approach.

Concluding, relativistic quantum dynamics was investigated for a multiple-particle system with clear interference fringes being identified and quantified. While quantum signatures in many-particle systems are likely to be washed out, our examples show that there is an intermediate regime in the number of involved particles with clear quantum effects for relativistic dynamics.

Financial support by the German Science Foundation (Nachwuchsgruppe within SFB276) is gratefully acknowledged.