text:
# Chemical composition of Galactic OB stars

## 1 Introduction

Massive blue stars strongly influence the chemical evolution of galaxies by their interaction with the Interstellar Medium, returning to it the material from their atmospheres and interiors. In recent years, there is growing evidence that their atmospheres might be enriched with CNO processed material even in the early stages of their lives, contrary in general to the predictions made by the classical models of stellar evolution (Schaller et al. Schaller92 (1992), Chiosi Chiosi98 (1998)).

Such models predict that only very massive stars, with masses higher than 60 $`M_{\mathrm{}}`$ on the ZAMS, lose enough mass in the wind to show, in the main sequence phase, the inner CNO enriched material (Schaller et al. Schaller92 (1992), Heger Heger98 (1998)).

For lower mass stars, the exposition of CNO processed material in the stellar atmosphere only occurs in later phases of the evolution, after the first dredge-up (Chiosi Chiosi98 (1998)).

However, recent investigations on the effect of stellar rotation on the structure and evolution of massive stars (Maeder & Meynet Mae&Mey00 (2000), Heger Heger98 (1998), Heger & Langer HegyLang00 (2000)) show that early atmospheric contamination can occur by mixing with the stellar interior, this contamination being more important as both the stellar mass and the initial rotational velocity increase.

The physical mechanism responsible for the mixing processes in massive stars may be other than rotation, but rotation appears as a natural driver for this mixing or, at least, as one of the first mechanisms whose consequences on mixing have to be explored (Meynet & Maeder Mey&Mae00 (2000)), especially considering the fact that massive stars are in general quite fast rotators.

Rotationally induced mixing has very important consequences on the evolution of the star, changing the evolutionary tracks towards higher luminosities, widening the main sequence in the HR diagram and making it last longer, and also changing the surface composition of the star. The effects of the CNO cycle show up in the atmosphere, with He and N enrichments and C and O depletions (Heger & Langer HegyLang00 (2000), Meynet & Maeder Mey&Mae00 (2000)). Like Heger & Langer (HegyLang00 (2000)) show in their fig. 9, He enrichments appear later in the stellar atmosphere, while CNO abundance changes are already present.

There is plenty of evidence in the literature of this CNO contamination in early type stars, one of the most important being the so–called He discrepancy (Herrero et al., Herrero92 (1992)). They find that a considerable fraction of Galactic OB stars are overabundant in He.

With regard to CNO abundances, a wealth of work also finds the trend in these elements abundances produced by the CNO cycle. Gies & Lambert (GyLambert92 (1992)) find N (and He) overabundances in a sample of B supergiants in the Galaxy, while Venn (Venn95 (1995)) finds such trend in the abundances of CNO for A supergiants, but not in the He abundances, that she assumes to have solar values. Other more qualitative works, studying equivalent widths, also find this trend in the CNO line intensities of B supergiant stars (Lennon et al., Lennon93 (1993), McErlean et al., McErlean99 (1999)), but not in the He abundance.

Therefore, no clear correlation between the He overabundance and CNO contamination is found in the literature, maybe due to the fact that, as predicted by rotating evolutionary models, changes in the N/C and N/O ratios can appear before any He enrichment.

As can be seen from the references given, little work has been done on the calculation of CNO abundances of O stars. Only Schönberner et al. (Schon88 (1988)) have studied three stars from O8 to O9 and Pauldrach et al. (Paul94 (1994)) and Taresch et al. (Taresch97 (1997)) two individual Galactic early O stars, $`\zeta `$ Pup (O4 If) and HD 93129A (O3 If) respectively.

In the very numerous works on abundances of Galactic B stars (Kilian et al., Kil94 (1994), Vrancken et al., Vran00 (2000), Smartt et al., 1997b , Rolleston et al., Rolles00 (2000), Cunha & Lambert, CunhayL94 (1994), Daflon et al., Daflon99 (1999), Gummersbach et al., Gumm98 (1998)) there is only a marginal presence of late O stars in the samples studied. Here is where we make our contribution, studying the CNO abundances of O type stars in order to compare them with the He abundance.

Of course, abundances derived from spectral analysis are dependent on the model assumptions. For example, as demonstrated by Villamariz (Villamariz01 (2001)) the He abundances obtained for OB stars from NLTE plane-parallel and hydrostatic model atmospheres are lowered by $``$20 % (in the He abundance by number relative to the total H+He, $`ϵ`$(He)) due to the inclusion of line–blocking in the spectral synthesis. A similar effect is found by consideration of a microturbulence of 15 km s-1 in the synthesis of the H and He lines of supergiant stars (Villamariz & Herrero Villamariz00 (2000)), as found earlier for B supergiants by McErlean et al. (McErlean98 (1998)) and Smith & Howarth (Smith&Howarth98 (1998)). The combined result is that the He discrepancy is considerably reduced and it even disappears for some stars. However, it still remains for an important fraction of Galactic OB stars.

On the other hand, analyses with the unified code of Santolaya-Rey et al. (Santolaya97 (1997)), which accounts for the stellar wind and considers spherical extension of the atmosphere, show no systematic changes in the He abundance (see Herrero et al. Herrero00 (2000)). Indeed, for 2 of the 7 stars analysed in this work, higher He abundances are found (see also Herrero et al. Herrero95 (1995)).

In this work we determine the CNO abundances of four O9 stars with known He abundances, as a first step in the investigation of the possible link between rotation and surface CNO contamination in the whole range of OB stars. In sect. 2 we present the objects selected for this study and the equivalent width measurements. In sect. 3 we present the technical details of our methodology, and in sect. 4 we make a careful study of the uncertainties in our abundance determinations. Sect. 5 contains the absolute abundances of the three He normal stars that we analyse, and sect. 6 their relative values. Sect. 7 is devoted to the analysis of the fast rotator HD 191423 and finally the discussion and conclussions of our work are given in sects. 8 and 9 respectively.

## 2 The programme stars

We have studied four O9 stars, two dwarfs and one supergiant with normal He abundances, and one rapidly rotating giant with an important He overabundance ($`ϵ`$(He)=0.20, Herrero et al. Herrero92 (1992)).

Our intention is to make a line by line differential analysis within the four stars, as for example Monteverde et al. (Ilu00 (2000)) and Smartt et al. (Smartt96 (1996), 1997a ). The determination of stellar chemical compositions is known to be dependent on the details of the methodology used. That is why it must be done differentially within stars of similar physical conditions. In this way, systematic errors in the line abundances, due to the simplifications adopted in the model atoms for example, can be avoided.

The four stars have been observed by our group during different observing runs. However, the observing conditions, spectral ranges, resolutions and S/N values are very similar in all cases. HD 209975, HD 34078 and HD 191423 were observed with the 2.5m Isaac Newton Telescope in 1989 using the IDS spectrograph with the H2400B grating, resulting in a spectral resolution of 0.6 Å (see Herrero et al. Herrero92 (1992)), whereas HD 214680 was observed with the same telescope and configuration in a later run in 1992 (see Herrero et al. Herrero00 (2000)). The S/N ratios range from 200 to 300 in the Hβ region.

The stars were analyzed by Herrero et al. (Herrero92 (1992)) by means of NLTE plane-parallel and hydrostatic H/He model atmospheres. Here however we present new parameters with line-blocking included in the line formation calculations, obtained in Villamariz (Villamariz01 (2001)).

Table 1 lists the programme stars together with their parameters. The parameters of HD 209975 and HD 191423 are slightly changed when the microturbulence found from the CNO analysis is applied to the H/He analysis, see sections 5.3 and 7.1.

### 2.1 Equivalent width measurements

We have measured equivalent widths of the CNO lines listed in tables 2 and 3 for the three stars with normal He abundance. The metallic features of the spectrum of the fast rotator HD 191423 (see sect. 7) are very weak and blended because of the high rotational broadening, and for this star equivalent widths of isolated lines cannot be measured.

When identifying lines in stellar spectra, it is often not easy to know which blends are contributing to the spectrum at any given wavelength. To deal with this difficulty, we have considered all possible contributions to each wavelength listed in the VALD database (Kupka et al. Kupka99 (1999), Ryabchikova et al. Ryab99 (1999)) and also in the classical book of Striganov & Sventitsky S&S68 (1968), within $`\pm `$2 Å each side of the line. In tables 2 and 3 line identifications (in the first column) are given with all possible contributions.

The identification of each spectral line is based on its Doppler shift and on its intensity compared to that of the other members of the multiplet, i.e., we only assume the identification of a line when its Doppler shift is consistent with the rest of the identified lines and its relative intensity compares well to the remaining lines of the multiplet (see tables 2 and 3).

The approximate gaussian shape of the line gives an idea of the possibility of line blending, and we finally decide about the reliability of each line a posteriori, when all the line abundances for each ion are calculated (see fig. 1).

To measure the equivalent widths we use our own software developed in IDL. A least squares profile fitting procedure was used, with gaussian profiles fitting the line and polynomials of degree one or two to fit the local continuum.

Errors in the measurements due to the uncertainty in the position of the local continuum (estimated as $`\pm `$ 1/S/N) are quoted together with the equivalent widths. We have found that this is the most important source of uncertainty in the given values, among those coming from the use of gaussian fits, or from the fitting of the continuum.

We have only considered lines of C III, N III and O II for the abundance determination, because they are the most common ions of each element in our O9 spectra.

## 3 Method of analysis

For the three He normal stars for which equivalent widths of isolated lines are measurable, we make use of the classical method of curve of growth to determine the absolute abundances. For each star and element we proceed as follows. We construct a grid of model spectra for abundances 10 times, 5 times and 2 times higher and lower than solar, and for the solar abundance, each for 6 values of microturbulence: 0, 5, 10, 15, 20 and 25 km s-1.

In this way, we construct the curves of growth for each line that gives the line abundance for each value of microturbulence. With the abundances of all the lines identified in the stellar spectrum we obtain the preliminary absolute abundance of the element together with its particular value of microturbulence (see sect. 3.3). This value minimizes the dependence of the line abundances with the line strength in the abundance–equivalent width diagrams. These diagrams are also a diagnostic tool to check the reliability of the observed lines for the abundance determination (see fig. 1).

The model CNO line formation calculations make use of the model atmosphere and model populations of H and He given by the stellar parameters $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$ and $`ϵ`$(He) of each star (see table 1). They are held fixed while the new spectrum is being calculated, with the codes DETAIL and SURFACE (latest versions, Butler 1998, private communication).

It is interesting to note that we make a consistent synthesis for each star and element, for the particular physical conditions determined by the stellar parameters $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$ and $`ϵ`$(He). Such a detailed analysis is not common in the literature. We have studied the methodologies followed in the works with which we will compare our results: Gies & Lambert (GyLambert92 (1992)), Cunha & Lambert (CunhayL94 (1994)), Schönberner et al. (Schon88 (1988)) and Daflon et al. (Daflon99 (1999)) and only the latter two make a consistent synthesis of the metallic spectra.

Gies & Lambert (GyLambert92 (1992)) and Cunha & Lambert (CunhayL94 (1994)) make use of the NLTE metallic abundances of the Munich group, calculated with the LTE line-blanketed model atmospheres of Gold (Gold84 (1984)), while the H/He analysis was made with the LTE line-blanketed model atmospheres of Kurucz (Kurucz79 (1979)), with a more complete blanketing than that of Gold (Gold84 (1984)).

We feel that the additional effort is worthwhile in terms of consistency, specially for the relative abundances determined, which will hopefully remain stable whatever physical assumptions are changed in forthcoming works. Furthermore, Cunha & Lambert (CunhayL94 (1994)) find that using the same model atmospheres for the H/He and the metallic synthesis consistently, apart from leading to different results, removes the dependence of the C and N abundances on the stellar temperature.

### 3.1 The metallic line opacities

When synthesizing the CNO spectra, we have made extensive use of the line opacities from the Kurucz opacity distribution functions (ODF), at solar abundance and 2 kms-1 microturbulence (Kurucz Kurucz96 (1996)).

This is to account for the stellar flux blocking due to the metallic lines, which produces important effects on the synthesized spectra, as can be seen in fig. 2 for the particular case of O ii in HD 34078. Therefore, the use of the ODFs will lead to more appropriate abundances.

For the H/He spectral synthesis we have made use of a different set of metallic lines, from 228 to 912 Å (see Herrero et al. Herrero00 (2000) for further details) instead of the longer range (9 – 100000 Å) considered for CNO. However, in the former range we find the strongest contribution of line-blocking to the H/He occupation numbers, so the new opacities considered would yield the same H/He spectra.

### 3.2 The model atoms

The model atoms used are basically those described in Becker & Butler (Bec&But88 (1988)) for O ii and in Eber (Eber87 (1987)) for C iii. We use a combined N ii–N iii model atom which was constructed using the N ii model of Becker & Butler (Bec&But89 (1989)) and the N iii model of Butler (Butler84 (1984)).

We have included a few new transitions in the O ii and C iii model atoms (see table 4) that were observed in our spectra and are not present in the quoted models. The N ii and N iii photoionization cross sections have been updated with the Opacity Project atomic data (Cunto et al. Cunto93 (1993)). The effects of the new cross sections on the synthetic equivalent widths can be seen in fig. 3.

To calculate the Voigt profile of these transitions, we assume that the collisional broadening is much more important than the radiative natural broadening and omit the latter. We adopt the formula of Cowley (Cow71 (1971)) for the collisional parameters. At the high temperatures of our models, this is a good assumption because Doppler broadening dominates the profiles.

During the abundance analysis we have excluded the C iii 4673 line from the calculation of the C abundance of the three He normal stars, because it yields considerably higher abundances than the remainder of the lines. This indicates either that the line is not being well synthesized, or perhaps that the $`\mathrm{log}`$ gf value is underestimated. The rest of the new transitions give line abundances consistent with the “old” ones.

### 3.3 Microturbulence

A wealth of work in the literature (see Vrancken Vran97a (1997) and McErlean et al. McErlean99 (1999) for example) show how the value of microturbulence can depend on which element is considered in a given atmosphere. The two previous works in particular show how, for B stars, the value obtained from the O ii lines is systematically higher than that obtained from other ions, such as Siiii.

We take this into account and make independent determinations of the abundance together with the microturbulence value for each element.

In synthesizing the CNO spectra it is necessary to decide how to treat microturbulence, i. e., whether to include it in both the calculation of the populations and in the line formation (in both DETAIL and SURFACE), or only in the line formation (only in SURFACE).

Vrancken (Vran97a (1997)) finds that for the high values of microturbulence of the giants of her sample, the inclusion in both DETAIL and SURFACE can lead to lower abundances by 0.1-0.2 dex. As we determine differential abundances, we do not include microturbulence in DETAIL in the present work.

## 4 Abundance uncertainties

The determination of the abundance uncertainties is as important as the abundances themselves.

Our absolute abundances are determined from the mean value of the individual line abundances (each of them calculated for the mean microturbulence, see sect. 5) of all the well behaved lines in the abundance–equivalent width diagrams (see fig. 1). In order to consider all the possible sources of uncertainty in this measurement we proceed following Israelian et al. (Isra98 (1998)).

The standard deviation of the line abundances, when the number of lines is large enough and they are weak to strong, will give us the uncertainty from three different sources: first from the uncertainty in the position of the continuum, second from the mean value of microturbulence adopted for the individual abundances, and third from the limitations of the model atoms, or in fact of the entire spectral synthesis process.

The other sources of uncertainty come from the uncertainty in the stellar parameters, $`\mathrm{\Delta }T_{\mathrm{eff}}`$, $`\mathrm{\Delta }\mathrm{log}g`$ and $`\mathrm{\Delta }ϵ`$(He) ($`\pm `$1 000 K, $`\pm `$0.1 dex and $`\pm `$25 % respectively, see Villamariz Villamariz01 (2001)).

Villamariz (Villamariz01 (2001)) has shown how H/He model spectra of lower gravities, those suitable for supergiants, are the most sensitive to changes in the model calculations. We expect this also to apply to the metallic spectra, so we estimate the line abundance uncertainties for our supergiant HD 209975 and take them as maximum values of the uncertainty for the whole O9 sample.

The stellar parameters of HD 209975 are: $`T_{\mathrm{eff}}`$=34 000 K, $`\mathrm{log}g`$=3.20 and $`ϵ`$(He)=0.08. For the calculation of the line uncertainties we proceed as follows. We construct a grid of six model atmospheres (and H/He populations), putting each parameter in turn to one of its two extreme values and fixing the other two parameters to the central values. With each model of the grid we perform all the abundance calculations as for the central parameters, to obtain new line abundances (see tables 5 to 7).

In this way we calculate the line uncertainties due to each parameter independently, and combine them quadratically to find the final line uncertainties due to the errors in the stellar parameters ($`\mathrm{\Delta }ϵ_\mathrm{i}^{\mathrm{mod}}`$). In particular, we have taken the contributions towards higher and lower abundances separately:

$$\mathrm{\Delta }ϵ_\mathrm{i}^{\mathrm{mod}\pm }=\sqrt{\mathrm{\Delta }^2ϵ_\mathrm{i}^{T_{\mathrm{eff}}\pm }+\mathrm{\Delta }^2ϵ_\mathrm{i}^{\mathrm{log}\mathrm{g}\pm }+\mathrm{\Delta }^2ϵ_\mathrm{i}^{ϵ(\mathrm{He})\pm }},$$
(1)
and $`\mathrm{\Delta }ϵ_\mathrm{i}^{\mathrm{mod}}`$ is the mean of those two values (see tables 5 to 7).

Again following Israelian et al. (Isra98 (1998)), these individual uncertainties are used to compute the mean uncertainty for each element ($`\mathrm{\Delta }ϵ(\mathrm{X})^{\mathrm{mod}}`$), which is then quadratically combined with the standard deviation to give the final uncertainty in the absolute abundance:

$$\mathrm{\Delta }ϵ(\mathrm{X})=\sqrt{\sigma ^2+\mathrm{\Delta }^2ϵ^{\mathrm{mod}}}.$$
(2)
Remember that while the standard deviation in the line abundances is individually calculated for each star and element, the uncertainty due to the model parameters is the same for the whole sample.

## 5 Absolute abundances of the He normal stars

We begin the abundance analysis by determining preliminary absolute abundances and microturbulence values for each element as explained in sect. 3. The line abundance–equivalent width diagrams give us the set of lines suitable for the abundance determination, according to the dispersion found in each case.

With the mean value of microturbulence given by the values found for the three elements, we compute the final line abundances and absolute abundances, given in table 11. The final absolute abundances are the mean of the individual values, and its uncertainty is the combination of the standard deviation and the error coming from the stellar parameters as explained in sect. 4.

Below we give details on the preliminary abundance determination for each star, and then we give the final absolute abundances for the three stars in table 11.

### 5.1 HD 214680 (10 Lac), O9 V

This star belongs to Lac OB1, an association at 368 pc from the Sun according to the Hipparcos data (de Zeeuw et al. Zee99 (1999)). Kane et al. (Kane80 (1980)) have measured LTE abundances for two B2 stars of the association, and Gies & Lambert (GyLambert92 (1992)) and Schönberner et al. (Schon88 (1988)) for HD 214680.

We only have a few lines for each element in common with Gies & Lambert (GyLambert92 (1992)), and the comparison with our equivalent widths shows that theirs are 5 to 30 % lower than ours. However, this difference is what they find when comparing their values to the literature, and they argue that the reason must be their underestimation of the far wings of the line when computing the equivalent width.

Schönberner et al. (Schon88 (1988)) have a total of nine lines of CNO in common with us, and even though their spectra are photographic and with a lower S/N, there is good agreement between our respective values.

The preliminary absolute abundances of CNO for this star are given in table 8. The microturbulence of 22.5 km s-1 indicated by the O ii lines is not a good value. These lines are too weak to be sensitive to microturbulence (the maximum equivalent width is only 93 mÅ, see table 3), and the scatter in the line abundances is not minimum for this value of microturbulence, so we do not take it as being representative.

Note that O ii in spite of being the most numerous ion in the spectrum has an uncertainty in the preliminary absolute abundance of the order of those for C and N, showing how the O ii line abundances are more scattered than for the other two elements.

The high value of microturbulence also found for C is an uncomfortable result. Kudritzki (Kud92 (1992)) and Lamers & Achmad (LamyAch94 (1994)) justify the high values found for supergiants as a result of the presence of stellar winds, but in the case of this dwarf star, without evidence in the spectrum of a strong wind, we cannot attribute this high value of microturbulence to this.

With the microturbulence indicated by N and C we find a mean microturbulence of 16 $`\pm `$ 9 km s-1, with which we compute the final line and absolute abundances.

### 5.2 HD 34078, O9.5 V

This star belongs to Ori OB1, which is located at about 500 pc from the Sun (de Zeeuw et al. Zee99 (1999)). Gies & Bolton (GiesyB86 (1986)) classify this star as a runaway, with its origin in Ori OB1 and an actual peculiar radial velocity of 50 km s-1.

Ori OB1 is a very well known association, with a wealth of abundance studies such as those of Kilian (Kil92 (1992)), Cunha & Lambert (CunhayL92 (1992), CunhayL94 (1994)) and Cunha et al. (Cunha98 (1998)). They all determine abundances of early B dwarf stars, with which we compare our results.

Gies & Lambert (GyLambert92 (1992)) have determined abundances for this star. As for HD 214680 our equivalent widths are higher than theirs, due to the reasons stated.

The preliminary absolute abundances of CNO for this star are given in table 9. The microturbulence of 17.5 km s-1 indicated by the O ii lines is again disregarded for the same reasons argued for HD 214680. And again, the line abundances for O show more scatter than for C and N.

With the microturbulence indicated by N and C we obtain a mean microturbulence of 14 $`\pm `$ 2 km s-1, with which we compute the final absolute abundances.

### 5.3 HD 209975, O9.5 Ib

This star belongs to Cep OB2, an association at 615 pc from the Sun (de Zeeuw et al. Zee99 (1999)). Daflon et al. (Daflon99 (1999)) have measured abundances for dwarfs from O9 to B2 in this association.

HD 209975 has a higher rotational velocity than the two dwarfs studied, and so its lines are more blended and diluted, leading to higher uncertainties in the equivalent widths and to a lower number of isolated lines suitable for the abundance analysis (see tables 2 and 3).

While the set of lines finally excluded for the abundance determination in the two dwarfs were coincident for C iii and N iii, and only differ in a single line for O ii, in the case of the supergiant HD 209975 we find that the set of lines that must be excluded according to their dispersion in the abundance–equivalent width diagrams, are not completely coincident. To assure consistency in our analysis, we decided to exclude both sets of lines for every element, so that only those lines with suitable scatter and also considered in the two dwarfs contribute to the abundance determination.

This leads to the preliminary abundances in table 10. The scatter in the O ii line abundances is especially severe in this star: compare the $`\sigma /\sqrt{\mathrm{n}}`$ values for the dwarfs in tables 8 and 9 of 0.04 and 0.05, to the 0.12 value for HD 209975. Again the microturbulence value indicated by O is not considered to be representative.

With a mean microturbulence of $`\xi `$=20 $`\pm `$ 10 km s-1 we compute the final absolute abundances for HD 209975 given in table 11 ($`\mathrm{log}ϵ_{\mathrm{i1}}`$). As seen in Villamariz & Herrero (Villamariz00 (2000)), the use of a value of $`\xi `$=15 km s-1 for the H/He spectral synthesis of models suitable for supergiants can lead to changes in the stellar parameters of the order of our standard error box of 1 000 K in $`T_{\mathrm{eff}}`$, 0.1 dex in $`\mathrm{log}g`$ and 25 % in $`ϵ`$(He). The changes in the stellar parameters influence the line abundances determined as explained in sect 4, so we apply the correspondent corrections to these abundances due to the value of microturbulence found for this star.

The changes in the parameters are adopted from those found for the O9 supergiant HD 210809 in Villamariz & Herrero (Villamariz00 (2000)): $`\mathrm{\Delta }T_{\mathrm{eff}}`$=+1 000 K, $`\mathrm{\Delta }\mathrm{log}g`$=+0.05 dex and $`\mathrm{\Delta }ϵ`$(He)=— 20 %. With them we compute the line abundances of the second column in table 11 for HD 209975, which are the final values. As seen in the table these results are more consistent with the two dwarfs HD 214680 and HD 34078.

## 6 Relative abundances of the He normal stars

As explained in sect. 2, computing differential line abundances will give us the best estimates of the differences in the chemical compositions of the stars studied. For stars of similar parameters, and therefore similar physical conditions, the systematic errors in the line abundances will be the same for all stars. The errors will thus be reduced by the use of the differential method.

Among the three He normal stars, the two dwarfs HD 214680 and HD 34078, with projected rotational velocities of 50 and 40 km s-1 respectively, are assumed to provide the reference values for the abundances of CNO in unmixed stars. A comparison of their abundances will give information about the consistency of our method, which should yield similar values for both, and therefore also about the scatter that can be expected between stars of similar chemical compositions. Remember that both stars are in the solar neighbourhood (de Zeeuw et al. Zee99 (1999)), so they are not expected to have had different initial chemical compositions.

The supergiant HD 209975, also from the solar neighbourhood and with v sin $`i`$= 100 km s-1, is assumed to provide the CNO reference abundances for stars of lower gravities. The differential analysis to the dwarf stars will then tell us whether its slightly different physical conditions, basically the more extensive atmosphere with lower surface gravity, introduces any changes in the reference CNO abundances.

The material of this star should not be mixed. Its He abundance is normal and its rotational velocity is lower than the 200 km s-1 considered as the approximate initial minimum value of the surface equatorial rotation for efficient mixing (Meynet & Maeder Mey&Mae00 (2000)). But we must keep in mind that the 100 km s-1 is a projected value, and also that the initial rotational velocities of massive stars decrease due to the loss of angular momentum in the stellar wind (Meynet & Maeder Mey&Mae00 (2000), Heger & Langer HegyLang00 (2000)), so the hypothesis of unmixed material must be checked. In adittion, a normal atmospheric He abundance is compatible with increments in the N/C and N/O ratios in the first stages of the main sequence evolution (Heger & Langer HegyLang00 (2000)).

In table 12 we give the line by line differential analysis of the three He normal stars. The final relative abundances are the simple mean of the individual values, and the standard deviation its corresponding uncertainty.

This is what we find:

The two dwarf stars have CNO abundances indistinguishable within the uncertainties, with a scatter in the CN values of 0.1 dex (with HD 34078 more abundant than HD 214680) and the same O abundance.

The supergiant HD 209975 also has CNO abundances that are indistinguishable from those of the dwarfs. The scatter in CN is also 0.1 dex and O abundance is again the same. It must be noted however, that the dispersion of the CN abundances goes in the direction of the CNO contamination: lower C abundance together with higher N abundance (see also figs. 11 and 12).

## 7 The fast rotator HD 191423

This star belongs to the Cyg OB8 association, with a distance modulus of 11.8mag (Humphreys Hum78 (1978)) that situates it at 2.3 kpc from the Sun. HD 191423 is a very fast rotator, v$`\mathrm{sin}i`$=450 km s-1, with the weak metallic lines completely blended, so the method used for the other three stars cannot be used in this case (see fig. 4). We perform a spectral synthesis and as we analyse the H/He spectra, fitting our model lines to the observed ones to find the stellar parameters $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$ and $`ϵ`$(He), we find the best fits to the observed blends that give us the stellar CNO abundances. This methodology is widely used in stellar abundance determinations when the lines available are not isolated (see for example Pauldrach et al., Paul94 (1994) and Taresch et al., Taresch97 (1997) for early O stars, and García-López et al., Ramon95 (1995) and Israelian et al., Isra98 (1998) for late type stars).

We have identified all the metallic spectral features present in the spectrum of the star. Those that correspond to lines used for the three He normal stars are now selected for the abundance analysis. See fig. 4 for the four metallic spectral features available:

(1): C iii–O ii 4067–82

(2): C iii–N iii–O ii 4379

(3): C iii–N iii 4505–21

(4): C iii–N iii–O ii 4620–60

We have made a careful search of all the possible contributions to these blends, as explained in sect. 3 when identifying lines in the stellar spectra. We find that the dominant contributions to each blend are those quoted, considering the lines that we have found in the other three stars in the same ranges, and also the synthetic spectra (see sect. 7.2). We cannot make a line by line differential abundance analysis as for the other three stars, but for the fast rotator we make use of the lines that were suitable for the abundance determinations in the reference stars.

Blend (4) is in emission in Of stars together with the He ii 4686 line (Walborn Wal71 (1971)). In the case of stars with weaker winds, like HD 191423, it can show some filling by the wind emission even while being in absorption, so we do not expect to obtain good fits for this blend.

In order to synthesize the complete CNO+H/He spectrum we proceed as follows. First we calculate individually the populations of each element, making use of DETAIL and the model atoms quoted in sect. 3.2. Then we use a combined CNO+H/He model atom for the line formation calculations with SURFACE to reproduce the blends observed.

The last step in this process, before the comparison with the stellar spectrum, is the convolution of the synthetic profiles with a standard rotational broadening function (and of course with the instrumental profile, negligible in this case). This procedure is not very suitable for objects with such high rotational velocities, where surface deformation and gravity darkening become of considerable importance. Further work must be done to improve our treatment of rotation.

### 7.1 Microturbulence and new stellar parameters

As the entire CNO+H/He spectrum is calculated simultaneously, we have to make use of one particular value of microturbulence for all the elements. We consider the mean value found for the supergiant HD 209975, $`\xi `$= 20 km s-1, because of its similarity with our fast rotator.

As seen for HD 209975, this value of the microturbulence for a giant star like HD 191423 can lead to changes in the stellar parameters, which have been determined for $`\xi `$= 0 km s-1. Therefore we proceed to determine new stellar parameters $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$ and $`ϵ`$(He) at $`\xi `$= 20 km s-1.

We start with a model spectrum at $`\xi `$=20 km s-1 and the initial parameters of the star: $`T_{\mathrm{eff}}`$=34 000 K, $`\mathrm{log}g`$=3.40 and $`ϵ`$(He)=0.20 (see fig. 5, dotted line) and we compute new model spectra in order to find the best fit for the observed H/He lines (for details on the lines selected for the H/He analysis, see Herrero et al., Herrero92 (1992) and Villamariz, Villamariz01 (2001)).

In fig. 5 we see the final stage of the fitting, yielding the new parameters: $`T_{\mathrm{eff}}`$=35 000 K, $`\mathrm{log}g`$=3.40, $`ϵ`$(He)=0.12 and $`\xi `$= 20 km s-1 (dash-dotted). It is interesting to note in the figure how the lines He i 4387 and He ii 4541 cannot be well reproduced with the H/He model atom, since they are blended with metallic lines (see fig. 6). In the following spectral synthesis we see how the fits improve with the use of the combined CNO+H/He model atom (see figures 7 and 8). The blend of He ii 4541 with Si iii 4542 is not reproduced by this combined model atom.

With the new parameters He i 4471 shows the dilution effect present in the spectra of late O and early B giants and supergiants (see Voels et al. Voels89 (1989), McErlean et al. McErlean98 (1998), Herrero et al. Herrero92 (1992)), and He ii 4686 shows some filling with wind emission.

On the basis of these new parameters we compute all the CNO+H/He synthetic spectra.

### 7.2 CNO abundances in HD 191423, O9 III:n

To start the abundance analysis of the fast rotator HD 191423 we begin with the abundances found for the reference supergiant HD 209975:

$`\mathrm{log}ϵ(\mathrm{C})=3.98\pm 0.24\mathrm{log}\frac{ϵ(\mathrm{C})_{}}{5}=4.19`$

$`\mathrm{log}ϵ(\mathrm{N})=4.01\pm 0.17\mathrm{log}ϵ(\mathrm{N})_{}=4.07`$

$`\mathrm{log}ϵ(\mathrm{O})=3.46\pm 0.48\mathrm{log}\frac{ϵ(\mathrm{O})_{}}{5}=3.87`$111In fact this value is closer to $`\mathrm{log}\frac{ϵ(\mathrm{O})_{}}{2}=3.47`$, but a transcription error made us start with the quoted value. Nevertheless our results are not dependent on this. Solar abundances are taken from Grevesse et al. (Grev96 (1996))

With the CNO abundances of our standard grid of 10, 5 and 2 times higher and lower than solar that are closest to the abundances of HD 209975 we compute the first model spectrum (model M0, see fig. 7, short–dashed line).

The blend labeled as (1) is the only one that gives us information about the abundances of O and C. To find a good fit, starting from M0, we must lower those abundances, which we do in models M1 to M3:

M0: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{5}`$, $`ϵ(\mathrm{N})=ϵ(\mathrm{N})_{}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{5}`$

M1: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{10}`$, $`ϵ(\mathrm{N})=ϵ(\mathrm{N})_{}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{5}`$

M2: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{5}`$, $`ϵ(\mathrm{N})=ϵ(\mathrm{N})_{}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{10}`$

M3: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{10}`$, $`ϵ(\mathrm{N})=ϵ(\mathrm{N})_{}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{10}`$

Blends (2) and (3) have O ii components, but as can be seen in fig. 7, they are dominated by the N iii component, the only line that remains at a fixed abundance in models M0 to M3.

Looking at the fits in fig. 7 we obtain the C and O abundances of HD 191423: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{10}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{5}`$ (model M1, dash–dotted), which are held fixed while we search for the N abundance, with models M4 and M5:

M4: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{10}`$, $`ϵ(\mathrm{N})=2.5ϵ(\mathrm{N})_{}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{5}`$

M5: $`ϵ(\mathrm{C})=\frac{ϵ(\mathrm{C})_{}}{10}`$, $`ϵ(\mathrm{N})=5ϵ(\mathrm{N})_{}`$, $`ϵ(\mathrm{O})=\frac{ϵ(\mathrm{O})_{}}{5}`$

In fig. 8 we see that it is not possible to find a good fit for both blends (2) and (3) with a single value of the N abundance. While blend (2) and the N iii line embedded in He i 4387 are well reproduced for $`ϵ(\mathrm{N})=2.5ϵ(\mathrm{N})_{}`$, blend (3) is better reproduced by a higher abundance, $`ϵ(\mathrm{N})=5ϵ(\mathrm{N})_{}`$. Furthermore, the N iii line embedded in He ii 4541 also gives contradictory information, because it is too strong with the three N abundances considered. We then find that the N abundance of this star may be between 2.5 and 5 $`ϵ(\mathrm{N})_{}`$.

Blend (4) shows evidence of wind emission, as with He ii 4686.

Adopting the abundance uncertainties of HD 209975 we give the absolute CNO abundances of the fast rotator HD 191423:

$`\mathrm{log}ϵ`$(C)=—4.48$`\pm `$0.24

$`\mathrm{log}ϵ`$(N) $``$\[—3.70,—3.40\]$`\pm `$0.17

$`\mathrm{log}ϵ`$(O)=—3.82$`\pm `$0.48

From now on we consider the most conservative value for the N abundance, the lower value, to be the N abundance of HD 191423.

## 8 Discussion

In fig. 9 and in table 13 we compare our results with those of Gies & Lambert (GyLambert92 (1992)) and Schönberner et al. (Schon88 (1988)) for HD 214680, with Gies & Lambert (GyLambert92 (1992)) for HD 34078, and with Cunha & Lambert (CunhayL94 (1994)) and Daflon et al. (Daflon99 (1999)) for the associations Ori OB1 and Cep OB2 respectively. In Cunha & Lambert (CunhayL92 (1992), CunhayL94 (1994)) it is shown how the two youngest subgroups of Ori OB1 (Ic and Id) are enriched in O and Si with respect to the other two older subgroups (Ia and Ib). We compare our results with the mean abundances of the older, non-enriched stars. Also the solar abundances of Grevesse et al. (Grev96 (1996)) are quoted.

Due to the strong dependence of the abundance determinations on the particular methodologies followed (see our discussion on this in sects. 3 and 6), the direct comparison of absolute abundances obtained in different works must be done with care, considering a possible “methodological bias”.

For the two dwarfs we see that our results are more consistent with each other than those of Gies & Lambert (GyLambert92 (1992)), which show considerable differences in the N and O abundances of these two similar stars (see fig. 9, crosses). They do not determine the C abundances. Schönberner et al. (Schon88 (1988)) find CNO abundances of HD 214680 to be higher than ours by 0.3 dex (see fig. 9, plus signs), an effect that can be explained by the ODFs used in our work, which as seen in sect. 3.1, lead to lower abundances.

The comparison with Cunha & Lambert (CunhayL94 (1994)) and Daflon et al. (Daflon99 (1999)) shows good agreement considering the uncertainties, and it is remarkable that in both their works and our dwarf stars, the C and N abundances are more dispersed than the O abundance, which is practically the same in all cases.

The CNO abundances that we obtain for our two reference dwarf stars are lower than the solar values by 0.1–0.4 dex, also in good agreement with the abundances of other unmixed early B dwarfs in the solar neighbourhood (Kilian Kil92 (1992), Kil94 (1994), Vrancken Vran97a (1997)). The results therefore support the hypothesis that the Sun is enriched in these elements relative to its neighbour early type stars.

The CNO abundances of our supergiant HD 209975 are different from those of the two reference dwarfs by $``$0.1 dex in the direction of CNO contamination: lower C and O and higher N abundances. However, these small differences are within our uncertainties (also in the relative values, see tab. 12) , so we believe that the differences can be due to the different physical conditions of the star relative to the dwarfs, as we indicated in sect. 6, and not to real different chemical compositions (remember that it also formed in the solar neighbourhood).

This star gives the reference CNO abundances of unmixed O9 giants and supergiants, but we will have to study more low gravity stars in order to establish these reference abundances reliably.

The CNO absolute abundances of our fast rotator HD 191423 could be showing the effect of possible internal mixing together with the different original composition of the cloud from which it was formed. It belongs to the Cyg OB8 association, that according to Melnick & Efremov (Mel&Efre95 (1995)) forms a group with the associations Cyg OB1 and Cyg OB9. Adopting for our star the galactocentric distance of Cyg OB1 (RGC=7.8 kpc, Massey et al. Mass95 (1995)), we see that it is closer to the galactic center than the Sun, at 8.5 kpc (Gummersbach et al. Gumm98 (1998)).

Considering the Galactic abundance gradients found by several authors such as Smartt & Rolleston (1997b ), Gummersbach et al. (Gumm98 (1998)) and Rolleston et al. (Rolles00 (2000)), we must correct for the different initial compositions before comparing its abundances with those of our reference stars.

We do so with the most recent gradients of Rolleston et al. (Rolles00 (2000)):

$`\mathrm{\Delta }`$C= —0.07 dex/kpc $`\times `$ (7.8–8.5)kpc = + 0.05 dex

$`\mathrm{\Delta }`$N= —0.09 dex/kpc $`\times `$ 0.7 kpc = + 0.06 dex

$`\mathrm{\Delta }`$O= —0.067 dex/kpc $`\times `$ 0.7 kpc = + 0.05 dex

that yield the second set of abundances represented in fig. 9. We see that the corrections are small and tend to reduce the N overabundance and to increase the C and O underabundances. Anyway the trend obtained in the CNO abundances of this fast rotator is in the direction of CNO contamination.

Massive stellar evolution with rotation predicts that the mixing is more efficient for higher mass stars and for higher initial rotational velocities, with the surface CNO contamination growing as the star ages. Another prediction is that the rotation of the star decreases with time, due to the loss of angular momentum in the stellar wind (Meynet & Maeder Mey&Mae00 (2000), Heger & Langer HegyLang00 (2000)).

The actual spectroscopic masses of our sample stars are the following:

HD 214680, 18 $`M_{\mathrm{}}`$

HD 34078, 21 $`M_{\mathrm{}}`$

HD 209975, 20 $`M_{\mathrm{}}`$

HD 191423, 29 $`M_{\mathrm{}}`$

and in fig. 10 we can see in the HR diagram that our two dwarf stars are certainly young, while the supergiant HD 209975 and the fast rotator are more evolved.

Therefore, the surface CNO contamination found in HD 191423 is in agreement with the evolutionary predictions: the star started its main sequence lifetime with a very high initial rotational velocity and it is now showing the corresponding surface CNO contamination.

In order to obtain more quantitative conclusions we compare our N/C and N/O ratios with the evolution of these quantities for certain evolutionary models, in figs. 11 and 12 (Meynet 2001, private communication).

We refer the ratios to their initial values, that for the evolutionary models is the solar composition and for our sample stars are those of the reference He normal stars.

Although we consider the CNO abundances of HD 209975 as the reference unmixed values for low gravity stars, we show in figs. 11 and 12 its ratios refered to the reference dwarf HD 214680.

They are compatible with the evolution of a rotating star with initial mass around 20 $`M_{\mathrm{}}`$, that has evolved until its actual (projected) rotational velocity. Therefore, we leave open the question of the possible CNO contamination of this object.

For HD 191423 we can only say that the star must have started its main sequence lifetime with an extremely high rotational velocity, that leads to the present day projected value and surface CNO contamination.

Therefore we find that it is possible to explain our results in the frame of massive stellar evolution with rotation, but we are still far from obtaining quantitative conclusions, due basically to the problem of the projection of the measurable rotational velocities and to the small sample considered.

## 9 Conclusions

We have studied a sample of four Galactic O9 stars in order to investigate the relationship between stellar rotation and surface CNO contamination.

The reference CNO abundances for unmixed O9 stars are given by the three He normal stars studied, all located in the solar neighbourhood. Their chemical compositions are indistinguishable within our uncertainties. These abundances are in agreement with those found in the solar neighbourhood for unmixed early B stars, with values about 0.2 dex below solar.

The He and CNO abundances of the fast rotator HD 191423 show the trend expected from CNO contamination: He and N overabundance and C and O depletions. For the reference supergiant HD 209975 we question whether it shows also some surface contamination.

The comparison of our observed N/C and N/O ratios with the predictions of the evolutionary models with rotation of Meynet & Maeder (Mey&Mae00 (2000)) show a qualitative agreement.

Nevertheless, we must be aware of the limitations of this result, which represents only a first step in the investigation of the relationship between stellar rotation and surface CNO contamination in massive O stars.

###### Acknowledgements.

We wish to thank N. Przybilla, R.J. García López, G. Meynet and N. Langer for providing us with data necessary for this work and for many valuable suggestions and discussions. AH acknowledges support for this work by the spanish DGES under project PB97-1438-C02-01 and by the Gobierno Autónomo de Canarias under project PI1999/008.


format: text

text:
# P-Cygni Type Ly𝛼 from Starburst Galaxies

## 1 Introduction

Ly$`\alpha `$ emission lines render us a very useful redshift indicator for starburst galaxies at high redshift, because their rest-frame equivalent width was predicted to be $`50200\mathrm{\AA }`$ (Charlot & Fall 1993). Thommes et al. (1998) suggested that the existence of strong Ly$`\alpha `$ emission is a necessary condition for primeval galaxies. However, dust enrichment is believed to be a fast process in the early evolution of galaxies, and therefore Ly$`\alpha `$ should suffer severe extinction by dust. It has been argued that even a small amount of dust in a very thick hydrogen medium can destroy Ly$`\alpha `$ photons because their path lengths before escape increase tremendously due to a huge number of resonance scatterings. Hence, dust absorption is believed to be the main reason for the early null detection of Ly$`\alpha `$ emission from high-$`z`$ galaxies by narrow band imaging surveys.

The Lyman break method has been very successful among a few programs to detect high-$`z`$ galaxies. Steidel et al. (1996) presented the rest-frame UV spectra of primeval galaxies at $`z=3`$. They continued to search for even more distant galaxies, and secured the spectra of galaxies at the redshift $`z=4`$ (Steidel et al. 1999). The UV spectra of those primeval galaxies are characterized by a flat UV continuum and the three types of the Ly$`\alpha `$ feature: symmetric emission, broad absorption in the wings, and P-Cygni type asymmetric emission.

From HST GHRS observations of nearby H II galaxies, Kunth et al. (1998) discovered that P-Cygni type Ly$`\alpha `$ emission lines are always accompanied by the low-ionization interstellar absorption lines blueshifted by $`100300\mathrm{km}\mathrm{s}^1`$ with respect to the systemic velocity, while broad wing absorptions of Ly$`\alpha `$ are accompanied by the much less blueshifted interstellar absorption lines. There are spectral similarities between nearby and high-$`z`$ starburst galaxies. Recently there have been a number of observations to obtain rest-frame optical nebular lines from high-$`z`$ galaxies (Teplitz et al. 2000a; Teplitz et al. 2000b; Pettini et al. 2001).

These observations support the idea that the kinematics of surrounding media plays a crucial role in the survival of Ly$`\alpha `$ photons in primeval galaxies, which is similar to the case in nearby starburst galaxies. When the systemic redshifts were determined by using the nebular lines, evidence for bulk motions of several hundred $`\mathrm{km}\mathrm{s}^1`$ is found (Pettini et al. 2001). Ahn & Lee (1998) and Lee & Ahn (1998) investigated the formation of asymmetric Ly$`\alpha `$ line profiles observed in the spectra of high-z galaxies by introducing the galactic superwind model.

It is generally believed that the Ly$`\alpha `$ emission is not a good indicator for the star-formation rate, because of dust extinction. However, Ly$`\alpha `$ emission often bears a nearly unique information on high-z starburst galaxies. Therefore, we have to understand the line formation mechanism, in order to establish the quantitative measure of star-formation rate.

We have been investigating the Ly$`\alpha `$ line formation in thick and static media of neutral hydrogen, and developed an accurate and efficient Monte Carlo method (Ahn, Lee, & Lee 2000, 2001, 2002). In this paper we apply our method to an expanding spherical supershell. In section 2 we describe our Monte Carlo code. We show the detailed mechanism of Ly$`\alpha `$ line transfer in expanding media in section 3, and the effects of both thickness and expansion velocity of the supershell in section 4. In the last section we summarize the results and discuss on the implication of our works.

## 2 Monte Carlo Method

Our model for starburst galaxies is illustrated in Figure 1. We assume that there is a central star-forming region surrounded by an expanding spherical supershell. We also assume that the Ly$`\alpha `$ source is the central super star cluster, because the Ly$`\alpha `$ luminosity of a typical starburst galaxy is equivalent to ionizing photons from $`10^310^4`$ O, B stars (González Delgado 1998; Ahn 2000) and the star-forming region is compact (Frye et al. 2002).

Expanding supershells are frequently seen in nearby starburst galaxies in a kiloparsec scale around the central star-forming region (Marlowe et al. 1995; Martin 1998). Recent HST observations revealed that the neutral column density of a supershell ranges $`N_{\mathrm{HI}}10^{18}10^{22}\mathrm{cm}^2`$, which is deduced from usual Voigt fitting analyses (Kunth et al. 1998). The Ly$`\alpha `$ line-centre optical depth is related with the H I column density $`N_{\mathrm{HI}}`$ by

$$\tau _01.41\left(\frac{b}{12.9\mathrm{km}\mathrm{s}^1}\right)^1\left[\frac{N_{\mathrm{HI}}}{10^{13}\mathrm{cm}^2}\right],$$
(1)
where $`b`$ is the Doppler width for turbulence motion. Therefore, we should deal with scattering media with the line-centre optical depths $`\tau _0=10^510^9`$. However, Monte Carlo methods are usually inefficient in this regime, because of the large computing time due to the large number of resonance scatterings required before escape. Hence we developed an accelerated Monte Carlo method in the previous paper (Ahn, Lee, & Lee 2002).

In an expanding supershell, some photons are scattered back into the opposite part of the supershell, and so they get redshifted. They experience a series of such back-scatterings before they escape the system. These processes enhance the escape probability of Ly$`\alpha `$ resonance photons in optically thick media than in static medium, and the number of scatterings before escape becomes much less than that in static media. Therefore, even though we want to solve the problem of radiative transfer in extremely thick media, it is sufficient for us to use the original code described in Ahn et al. (2001).

We illustrate the back-scattering process in Figure 1. Suppose that a photon scattered at $`P_1`$ to a direction $`𝐤\stackrel{}{P_1P_3}`$ is scattered at a next scattering site $`P_3`$ if the central region is also uniformly filled with neutral hydrogen atoms. Since no neutral hydrogen is assumed to be present in $`\overline{P_4P_5}`$, the photon traverses to $`P_2`$, where $`|\overline{P_5P_2}|=|\overline{P_4P_3}|`$. We determine the distance $`2d|\overline{P_4P_5}|`$ along the wave vector $`𝐤`$, and add $`|\overline{P_4P_5}|`$ to the path length $`|\overline{P_1P_3}|`$ in the code. Here we define the position vectors of the adjacent scatterings to be $`𝐫^{}`$ and $`𝐫^{\prime \prime }`$, $`r^{}=|𝐫^{}|`$, and the angle $`\theta \mathrm{}OP_1P_4`$. An elementary geometrical calculation gives

$$d=r^{}\mathrm{cos}\theta l,$$
(2)
where

$$\mathrm{cos}\theta =\frac{𝐫^{}𝐤}{|𝐫^{}|},$$
(3)
$`l|\overline{P_1P_4}|=r^{}\mathrm{cos}\theta \left[r^2\mathrm{cos}^2\theta r^2+r_{min}^2\right]^{1/2}`$, and $`|\overline{OP_4}|=|\overline{OP_5}|r_{min}`$.

The velocity shifts due to the bulk expansion of the medium should be taken into account. Referring to the geometry shown in Figure 1, the velocities $`𝐯_\mathrm{𝟏}`$ and $`𝐯_\mathrm{𝟐}`$ are given by

$$𝐯_\mathrm{𝟏}=V_{\mathrm{exp}}\mathrm{cos}\theta 𝐤,$$
(4)
and

$$𝐯_\mathrm{𝟐}=V_{\mathrm{exp}}\frac{sl+d}{r^{^{\prime \prime }}}𝐤,$$
(5)
where $`V_{\mathrm{exp}}`$ is the expansion velocity of the supershell and $`V_{\mathrm{exp}}>0`$ by definition. Here $`s|\overline{P_1P_3}|`$ is the path length corresponding to an optical depth $`\tau =\mathrm{ln}R`$, where $`R[0,1]`$ is a random number. We subtract frequency corresponding to the relative velocity $`𝚫𝐯=𝐯_\mathrm{𝟐}𝐯_\mathrm{𝟏}`$ from the frequency $`x`$ whenever a back-scattering occurs. Here we introduce the dimensionless frequency parameter $`x`$ defined by

$$x\frac{\mathrm{\Delta }\nu }{\mathrm{\Delta }\nu _D}=\frac{(\nu \nu _0)}{\mathrm{\Delta }\nu _D},$$
(6)
which describes the frequency shift from the line-centre frequency $`\nu _0`$ in units of the Doppler width $`\mathrm{\Delta }\nu _D\nu _0(b/c)`$. Here $`b`$ is the Doppler parameter that combines both turbulence and the thermal motions of the scattering medium, and $`c`$ is the speed of light. The detailed description of our code can be found in the previous papers (Ahn, Lee, & Lee 2000, 2001, 2002).

While scatterings occur in the supershell, we should also take into account the Doppler shift caused by the relative bulk motion between the two adjacent scattering positions. For the radial vectors of the adjacent scattering locations, $`𝐫_\mathrm{𝟏}`$ and $`𝐫_\mathrm{𝟐}`$, the expansion velocities projected into the wave vector should be considered. After a simple geometrical consideration, we obtain

$$\mathrm{\Delta }v=V_{\mathrm{exp}}\left(\frac{𝐫_\mathrm{𝟐}}{|𝐫_\mathrm{𝟐}|}\frac{𝐫_\mathrm{𝟏}}{|𝐫_\mathrm{𝟏}|}\right)𝐤.$$
(7)
We subtract this shift from the frequency of a Ly$`\alpha `$ photon whenever scatterings take place in the medium. When a Ly$`\alpha `$ photon escapes the medium by a back-scattering, we consider the redshifts due to the final back-scattering; i.e., only $`𝐯_\mathrm{𝟏}`$ in Eq.(4) is taken into account. We also calculate the path lengths of photons traversing in the supershell in order to estimate the amount of dust extinction. However, dust effect is considered very briefly in this paper, and we will concentrate on the Ly$`\alpha `$ line transfer in the expanding supershell.

## 3 The Formation of P-Cygni Type Ly$`\alpha `$ Profiles

We now discuss on typical profiles of Ly$`\alpha `$ resulting from the scattering processes discussed above. The expanding supershells in starburst galaxies have typical turbulence which is given by the Doppler parameter, $`b40\mathrm{km}\mathrm{s}^1`$ estimated from Voigt fittings (e.g. Kunth et al. 1998). Hence we assume that the Voigt parameter of the expanding supershell to be $`a=1.5\times 10^4(40\mathrm{km}\mathrm{s}^1/b)`$, where the Voigt parameter $`a\mathrm{\Gamma }/4\pi \mathrm{\Delta }\nu _D`$ with $`\mathrm{\Gamma }`$ being damping constant. Here we neglect the thermal motion, because the turbulent motion usually exceeds the thermal motion. The FWHM of the input profiles is set to be $`\mathrm{\Delta }x=2`$, because the optical nebular emission lines indicate that $`\mathrm{\Delta }v=70100\mathrm{km}\mathrm{s}^1`$ (Pettini et al. 2001). We consider a supershell with its expansion velocity being $`V_{\mathrm{exp}}300\mathrm{km}\mathrm{s}^1`$ and its column density $`10^{18}\mathrm{cm}^2<N_{\mathrm{HI}}<10^{21}\mathrm{cm}^2`$, in accordance with observations in nearby starbursts (Kunth et al. 1998) and high-z galaxies (Pettini et al. 2001). These two quantities satisfy the condition $`a\stackrel{~}{\tau }_0>10^3`$, which guarantees the validity of the analytic solution given by Adams (1972), Harrington (1973), and Neufeld (1990). Here we define the line centre-optical depth by $`\stackrel{~}{\tau }_0\sqrt{\pi }\tau _0`$, which is a little bit different from the previous researches.

In Figure 2 we show the emergent profile for the model. The number of photons used in the Monte Carlo calculation is 80,000, and the other physical parameters are assumed to be the typical values mentioned above. Especially we choose the neutral column density $`N_{\mathrm{HI}}=2\times 10^{20}\mathrm{cm}^2`$. The dotted line in the figure represents the average number of scatterings for a photon during its transfer, and the dashed line represents the average number of successive wing scatterings just before its escape. The solid line refers to the emergent Ly$`\alpha `$ profile. The vertical bars with numbers refer to the frequency of emission peaks corresponding to $`N`$-th back-scattering. We can see two kinds of peaks in the profile: one located at the bluer frequency, and the other three at the redder part.

We see in the figure that the weak peak at the blue part is composed of photons which have experienced very large number of scatterings in the wings as well as in the core. According to Adams (1972), the number of scatterings for those photons is $`N\tau _0`$. Hence we can see that these are photons transferred in the supershell by so-called excursions (Adams 1972). Some photons around $`x0`$ are also similar ones, which is evident by the large number of scatterings and wing scatterings.

These weak peaks are formed through the following processes. When photons with the line-centre frequency are incident upon the inner wall of the expanding supershell, its scattering is off-centred by an amount of $`x_{\mathrm{exp}}`$, where we define the dimensionless frequency corresponding to the expansion velocity $`x_{\mathrm{exp}}V_{\mathrm{exp}}/v_{\mathrm{th}}`$. Hence, the fraction of transmitted photons through the supershell is approximately given by Eq.(2.23) in Neufeld (1990), and the profile of transmitted photons is given by Eq. (2.32) in the same paper. Since $`x_{\mathrm{exp}}=7.65`$ in the above simulation, we can see from Eq.(2.23) in Neufeld (1990) that a very small amount of photons can transmit the optically very thick supershell and majority of photons are back-scattered. The transmitted photons are redistributed into both sides of $`x=x_{\mathrm{exp}}`$. The radiative transfer of those photons was studied very well both analytically (Neufeld 1990) and numerically (Ahn et al. 2002).

The three peaks at the red part are composed of back-scattered photons. Since the supershell is receding with respect to the incident photon, scatterings usually occur at wing frequencies if the opacity at the frequencies is sufficiently large. Here we define the effective scattering opacity at a wing frequency $`x`$ by $`\tau _e(x)`$. At the first scattering of Ly$`\alpha `$ photon with the expanding supershell, the effective optical depth should be measured at $`x_0+x_{\mathrm{exp}}`$, i.e $`\tau _e(x_0+x_{\mathrm{exp}})`$. Here we consider an input frequency $`x_0=0`$, for simplicity. Note that $`x_{\mathrm{exp}}<0`$. If $`\tau _e(x_0+x_{\mathrm{exp}})1`$, a back-scattering occurs at the wing frequency. Then, the back-scattered photons are redshifted to the frequency $`x=x_0+x_{\mathrm{exp}}`$, and they will scatter in the other part of the supershell along the back-scattered direction. If the supershell in that direction is optically thin, then those back-scattered photons can escape freely to form a prominent peak. In fact, that part of the supershell is also receding with respect to the back-scattered photons, and so the effective opacity now becomes $`\tau _e(x_0+2x_{\mathrm{exp}})`$. The opacity at this frequency is even smaller than that of previous scattering, and the escape probability increases. Hence some photons escape the system and form the first peak at $`x=x_0+x_{\mathrm{exp}}`$. The photons reflected into slant directions are redshifted by an amount smaller than those to the normal direction, and the optical depths in those direction becomes larger. Hence the escape probability of those photons become smaller than that of photons in the normal direction. These photons contribute to the broadening of each peak. It is also noticeable in the figure that the primary peak is not well matched with $`x_{\mathrm{exp}}`$, because its blue part is eroded by the scattering of very thick supershell whose scattering centre is located at $`x=x_{\mathrm{exp}}`$ in the figure.

Ly$`\alpha `$ photons repeat the above-mentioned back-scatterings until $`\tau _e`$ becomes small enough to permit most of the photons to escape the system. We see that photon’s effective scattering frequency at the $`N`$-th back-scattering is given by $`2(N1)x_{\mathrm{exp}}`$, where $`N>1`$. The frequency of photons escaping associated with $`N`$-th back-scattering is given by $`(2N3)x_{\mathrm{exp}}`$. When $`N=1`$, the situation is special, because the source is at rest. In this case, each frequency is $`x_{\mathrm{exp}}`$ and $`0`$, respectively.

We illustrate in Figure 3 how the effective opacity of successive back-scatterings varies, where we consider photons scattering back and forth along a diameter for simplicity. The effective frequencies of Ly$`\alpha `$ photons successively change, as is drwan by the vertical dotted lines in the figure. The numbers in the box for each dotted line refers to the order of back-scatterings denoted by $`N`$. Let us consider the case for $`\tau _0=10^6`$ which corresponds to the case in Figure 2. We can see that the effective opacity for $`N=1`$ is very large, and so most of Ly$`\alpha `$ photons are back-scattered. This means that there is a little amount of transmitted flux with $`x=0`$. The second or $`n=2`$ back-scattering happens at wing frequencies much further away from the $`N=1`$ frequency, and so the escape probability is enhanced, but it is still not large enough for most photons to escape the system. However, the emission peak forms at $`x=x_{\mathrm{exp}}`$, because the flux of incident Ly$`\alpha `$ photons is very large. In the next turn with $`N=3`$, the second peak at $`x=3x_{\mathrm{exp}}`$ forms due to the process similar to the $`N=2`$ case. When $`N=4`$, we expect that the other peak, but it appears very weak because most Ly$`\alpha `$ photons have already escaped from the system. From then on the effective opacities for the following back-scatterings are sufficiently small, but photons have been already exhausted. So no peaks at the higher frequencies appear in the emergent spectrum, as shown in Figure 2. We can summarize mathematically these explanation as follows.

$$I_1(0)=I_0f_\mathrm{T}(x_{\mathrm{exp}})0,$$
(8)
$`I_2(x_{\mathrm{exp}})`$ $`=`$ $`I_0\left[1f_\mathrm{T}(x_{\mathrm{exp}})\right]f_\mathrm{T}(2x_{\mathrm{exp}})`$
$`=`$ $`I_0f_\mathrm{R}(x_{\mathrm{exp}})f_\mathrm{T}(2x_{\mathrm{exp}}),`$
$$I_3(3x_{\mathrm{exp}})=I_0f_\mathrm{R}(x_{\mathrm{exp}})f_\mathrm{R}(2x_{\mathrm{exp}})f_\mathrm{T}(4x_{\mathrm{exp}}),$$
(10)
$$I_4(5x_{\mathrm{exp}})=I_0f_\mathrm{R}(x_{\mathrm{exp}})f_\mathrm{R}(2x_{\mathrm{exp}})f_\mathrm{R}(4x_{\mathrm{exp}})f_\mathrm{T}(6x_{\mathrm{exp}}),$$
(11)
and so on. Generally

$`I_N\left([2N3]x_{\mathrm{exp}}\right)=I_0f_\mathrm{R}(x_{\mathrm{exp}})`$ (12)
$`{\displaystyle \underset{n=2}{\overset{N1}{}}}f_\mathrm{R}\left([2n2]x_{\mathrm{exp}}\right)f_\mathrm{T}\left([2N2]x_{\mathrm{exp}}\right),`$
where the transmitted fraction ($`f_\mathrm{T}`$) and the reflected fraction ($`f_\mathrm{R}`$) are approximately described in Eq.(2.32) and Eq.(2.33), respectively of Neufeld (1990). Here $`I_0`$ is the input flux, and $`I_N`$ is the flux associated with the $`N`$-th back-scattering.

In order to see how much emergent Ly$`\alpha `$ photons are affected by dust extinction, we calculate the total path length of the emergent photons through the medium. Here we adopt $`N_{\mathrm{HI}}2\times 10^{20}\mathrm{cm}^2`$, $`b=40\mathrm{km}\mathrm{s}^1`$, and $`V_{\mathrm{exp}}300\mathrm{km}\mathrm{s}^1`$. We also assume that dust is uniformly distributed in the supershell.

In Figure 4 we show the spectral distribution of both the total path length of Ly$`\alpha `$ photons in the supershell ($`\tau _\mathrm{d}`$) and the total number of scatterings ($`N_{\mathrm{scat}}`$). The latter represented by solid circles in the inset is depicted in logarithmic scale, and shows clearly two groups. One group is located at the upper domain with $`\mathrm{log}_{10}N_{\mathrm{scat}}6`$, and they are photons having experienced excursions. The rest in the lower part with $`\mathrm{log}_{10}N_{\mathrm{scat}}1`$ are those having escaped by a number of back-scatterings.

The crosses represent, in units of $`\tau _0`$, the total path length along which a Ly$`\alpha `$ photon traversed in the supershell. We can see in the figure that most of emergent photons escape by back-scattering, which can be inferred from the small number of scatterings. It is also noticeable that the total path lengths are very diverse in the spectral distribution and most of photons traverses rather short path. Since there is no frequency-dependence in path lengths except for a weak trend of monotonous increases, the redder peaks in emergent profiles from dusty media can be destroyed at most a little bit more than the bluer peaks. This is puzzling because all the observed Ly$`\alpha `$ profiles shows single peak which must be the $`N=2`$ peak in our calculations. The solution of this puzzle may lie in the spatial distribution of dust in the system. If the inner radius of the dust shell is is smaller than that of neutral hydrogen, the dust extinction of the higher order peak can increase. We are now undertaking this subject, and the results will be published in the forthcoming paper.

## 4 Dependence on the Physical Parameters

The emergent profiles depends on a few parameters of scattering medium in our model. First of all, the turbulence is out of our interest, because they are constrained by observations and the emergent profile is dominantly formed by kinematics of the supershell in our case. Hence here we investigate how the Ly$`\alpha `$ peaks depend on the optical depths and the expansion velocity of the supershell.

### 4.1 Effects of supershell’s opacity

We present in Figure 5 the emergent profiles for a few cases with different line-centre optical depths, where the expansion velocity of the supershell is fixed to be $`x_{\mathrm{exp}}=7.65`$. The FWHM of input profiles is also set to be $`\mathrm{\Delta }x=2`$, and the Doppler parameter of the supershell is set to be $`b=40\mathrm{km}\mathrm{s}^1`$. The supershell is assumed to have no dust, and the number of input photons are set to be equal for each calculation or $`N_{\mathrm{ph}}=80,000`$. The arrows in the figure represent the peaks corresponding to $`x_{\mathrm{exp}}`$, $`3x_{\mathrm{exp}}`$, and $`5x_{\mathrm{exp}}`$. When the opacity of the supershell is not sufficiently large ($`N_{\mathrm{HI}}=2\times 10^{18},2\times 10^{19}`$), the input photons are transmitted without scatterings and the peak at $`x0`$ forms, whose FWHM increases by a small amount compared to that of input profile. We can see in the figure that the back-scattered photons form the primary peak at $`xx_{\mathrm{exp}}`$, the secondary peak at $`x3x_{\mathrm{exp}}`$, and the tertiary peak at $`x5x_{\mathrm{exp}}`$. As the opacity of the supershell decreases, the peaks at higher frequencies disappear accordingly. We have explained the mechanism of those variations in the previous section.

### 4.2 Effects of Kinematics

We also investigate the effects of the expansion velocity of the supershell. The results of Monte Carlo calculations are shown in Figure 6. We can see that the emission peaks gradually merge into a single peak as $`|x_{\mathrm{exp}}|`$ decreases. This is because the supershell behaves like a spherical medium which is static with respect to the Ly$`\alpha `$ source as $`|x_{\mathrm{exp}}|0`$. In an extremum with $`|x_{\mathrm{exp}}|<1`$, the emergent profiles will have symmetric double peaks which is same to those given in Ahn et al. (2002). If $`|x_{\mathrm{exp}}|`$ is much larger than unit, then we can see the emission peaks of back-scattered photons in the emergent profile. However, as $`|x_{\mathrm{exp}}|1`$, the frequency gaps between those peaks become close, and they eventually mixed up to be a single broad peak. We can also expect a heavy extinction for small $`|x_{\mathrm{exp}}|`$ cases.

It is also noticeable that the similar problem for the contraction supershell has a exact mirror-symmetry to the expanding supershell. Hence it is evident that we obtain the profiles to the origin $`x=0`$ for the contracting case.

## 5 Summary

We have investigated the formation of P-Cygni type Ly$`\alpha `$ by a supershell surrounding a star-forming region and outflowing in a bulk manner with expansion velocities of $`V_{\mathrm{exp}}100300\mathrm{km}\mathrm{s}^1`$. We have modified the Monte Carlo code developed in the previous paper (Ahn et al. 2001), and studied in detail the Ly$`\alpha `$ line transfer in the expanding supershell consists of neutral hydrogen gas. The Ly$`\alpha `$ line-centre optical depth of the supershell is very large and of order $`\tau _010^510^7`$. Its turbulence motion is assumed to be about $`40\mathrm{km}\mathrm{s}^1`$ in accordance with observations.

We found that the back-scatterings play a crucial role in the formation of asymmetric Ly$`\alpha `$ profiles. Most of Ly$`\alpha `$ photons escape the system efficiently through back-scatterings. Whenever it is scattered back by the expanding supershell, photon’s frequency is redshifted successively. This process promotes the redshift of Ly$`\alpha `$ photons. It is also naturally expected that the collapsing supershell promotes blueshift of Ly$`\alpha `$ photons, and the emergent profile for that case is simply mirror-symmetric to that of the expanding case. The emission peak formed by the $`N`$-th back-scattering with the supershell has a frequency of $`(2N3)x_{\mathrm{exp}}`$, where $`N>1`$. The number of peaks and their relative flux ratios in emergent profiles are determined by the interplay between the line-centre optical depths and the expansion velocities of supershells.

Recently optical nebular lines of about a dozen of Lyman break galaxies have been obtained by using near-IR spectrographs and large telescopes (Pettini et al. 2001; Teplitz et al. 2000a; Teplitz et al. 2000b). The Ly$`\alpha `$ spectra of a large number of primeval galaxies can be found in the literature (Steidel et al. 1996; Steidel et al. 1999; Lowenthal et al. 1997; Frye et al. 2002), which are reviewed in Stern & Spinrad (1999). Their UV spectra and Ly$`\alpha `$ emission lines show similar characteristics to those of nearby starburst galaxies. Their P-Cygni type Ly$`\alpha `$ emission lines are always accompanied by nebular lines blueshifted with respect to Ly$`\alpha `$. We have showed that these asymmetric Ly$`\alpha `$ can be formed by back-scatterings of Ly$`\alpha `$ photons by an expanding supershell.

The quantitative measurement of star-formation rates of starburst galaxies is one of the most crucial issues in modern cosmography. It is suggested that the expansion of surrounding media can help Ly$`\alpha `$ photons escape the star-forming regions more freely. In light of our works, it is very interesting to see that nearly all the observed P-Cygni type Ly$`\alpha `$ emissions from starburst galaxies do not show any higher order peaks which are formed by back-scatterings. This can certainly be attributed to dust extinction. However, it is very subtle problem if we want to quantify the escape fraction of Ly$`\alpha `$ photons for the case of the dusty and expanding supershell. Hence, we believe that very careful consideration of both abundance and spatial distribution of dust is needed to establish a realistic model describing the star-forming activity in starburst galaxies whether it is nearby or primeval. We are now undertaking this topic, and its results will be reported in a forthcoming paper.

## Acknowledgments

This work was done as a part of doctoral dissertation of SHA at the School of Earth and Environmental Sciences of Seoul National University financially supported by Brain Korea 21 of the Korean Ministry of Education. Also this work is completed with financial support of Korea Institute for Advanced Study, which is gratefully acknowledged. HML was supported by KOSEF Grant No. R01-1998-00023.


text:
# CHAOTIC DYNAMICS OF FOREST FIRES

## 1 Introduction

The problem of forest fires has got a wide audience  for its percolation origin  and rich connections to basic problems in computational physics, e.g. to the self-organized criticality. On the other hand, the problem is not entirely virtual, although experiments are not eagerly accepted by the society. The data can be found in annual statistics of forest rangers and insurance companies. Both institutions are obviously interested in predicting of forest fires, at least in the statistical sense. However, there are some indications that such predictions are not possible in long time scale — the problem is chaotic .

The present work is a continuation of Ref. . In that paper, a map $`p_np_{n+1}`$ which describes the forest density $`p`$ in subsequent time steps $`n`$ was proposed. The construction was based on the idea of reducing a probabilistic cellular automaton to the map, which is equivalent to the mean field approximation of the automaton rules . The map was designed not so as to reproduce the dynamics of the fire front, but rather to simulate the statistics of the forest density in subsequent years $`n`$. The forest was represented by a square lattice $`L\times L`$, occupied by trees with probability $`p`$. The above map was equivalent to a superposition of two functions, where each function transforms the forest concentration $`p`$. The growing process was represented by a parabolic function, $`p_{n+1}=p_n+rp_n(1p_n)`$, where $`r`$ is a parameter. Next, a season of fires was modeled by a reduction of $`p_{n+1}`$ by what was burned. This amount was demonstrated to be equal to the mean squared size of clusters of trees. By definition, a tree belonged to a cluster if it was a nearest neighbor of another tree which belongs to it. The neighborhood was of von Neumann kind.

The prescription of calculating the burned area contained the probability that a cluster is ignited. This probability can be defined “per area” or “per cluster”. We have obtained two kinds of statistics, which can be interpreted as follows: In the former case, a lightning strikes at a cell, occupied by a tree or not. The probability that the ignited tree belongs to a cluster is proportional to the size of this cluster, i.e. the number of trees in it. The burned area is also proportional to this size. Then, the totally burned area is the mean squared size of the cluster, as stated above. The weight of the cluster of size $`s_i`$ is $`w_i^{}=s_i/L^2`$. In the latter case, the fire is set intentionally, what means that a tree is ignited with probability one. The rest of argumentation is the same as in the former case. The appropriate weight in a such situation is

$$w_i^{\prime \prime }=s_i/\underset{i}{}s_i.$$
(1)
The former case ($`w^{}`$) was discussed thoroughly in Ref. . In the present work we focus on the latter case, where the weights are given by Eq. (1).

The paper is organized as follows: The next section is devoted to the calculation of the map $`p_np_{n+1}`$, which is obtained by an application of the HKA . In Section 3 we demonstrate the results on the Lyapunov exponent and the symbolic dynamics, which allows to obtain the cycle of the length three in the time evolution of the trees density $`p`$. There, we also refer to the experimental data. Our conclusions are given in the last section.

## 2 The model

The map $`p_np_{n+1}`$ is a superposition of two functions. The first one, $`pp+rp(1p)`$, describes the season of growing. As a result, trees appear in empty cells with the probability $`rp(1p)`$, where $`r`$ is a parameter from the range $`[0,1]`$. Afterwards, the burned area of the forest is calculated with the weights $`w_i^{\prime \prime }`$, i.e.

$$A(p)=\frac{1}{L^2}\underset{i}{}w_i^{\prime \prime }s_i$$
(2)
and it is subtracted from trees concentration after growth season, i.e. $`ppA(p)`$.

The cluster size distribution is obtained by means of HKA. With HKA  we are able to label all occupied sites on a $`L\times L`$ large square lattice in such a way that the sites with the same label belong to the same cluster and different labels are associated to different clusters. As a result, we obtain the cluster size distribution for $`1001`$ values of $`p`$, $`p_i=10^3i`$, $`i=0,1,2,\mathrm{},10^3`$. The simulation is carried out on a square lattice $`10^3\times 10^3`$ large, with a fraction of $`p+rp(1p)`$ occupied sites for different values of the initial trees concentration $`0p1`$. The map is got for 1001 equidistant values of $`p`$ as an average over hundred trial lattice configurations for each value of $`p_i`$. For other values of $`p`$ we use straight lines joining neighboring points $`p_i`$ and $`p_{i+1}`$. The obtained map is shown in Fig. 1 for three different values of the parameter $`r`$.

The method above suffers from the finite size of the lattice and the finite grid. In particular, we observe some shift of the maximum of the map for $`r=0`$ with respect to the percolation threshold for the square lattice. Instead to occur at $`p_c=0.59273`$ , we find the maximum near $`p=0.583`$. The rounding effect of the finite size reduces the map to the well-known unimodal map, where the period-doubling bifurcation must occur, at least at some range of the parameter $`r`$. That is why we extrapolate the shape of the curve to a tentative limit of infinite size (the so-called thermodynamic limit). The result for $`r=0`$ is

$$p_{n+1}=\{\begin{array}{cc}p_n\hfill & p_np_c,\hfill \\ p_c\mathrm{exp}(12.735\sqrt{p_np_c})\hfill & p_n>p_c.\hfill \end{array}$$
(3)
## 3 Results

The symbolic dynamics is known to be the most rigorous method to investigate chaos in one-dimensional maps . The basic tool is to assign symbolic words to the periodic cycles. Example giving, the cycle of length two is represented by CR, of length three — by CRL, etc. The appearance of the cycle three in the time evolution of a given variable is known to be a proof of the presence of cycles of all possible lengths. In our case the cycles of the variable $`p`$ are not superstable, because the curve given by Eq. (3) is not differentiable at $`p=p_c`$. In fact, the cycles are even not stable: we do not observe any windows of stability in the bifurcation diagram. Still, it is not difficult to find a cycle three by the method of “shooting”. Fixing the value of $`p_1`$, we change the parameter $`r`$ until $`p_4=p_1`$. An example is the cycle $`p=0.340000,0.419593,0.505974`$ for $`r=0.354694`$.

In Fig. 2 we show the values of the Lyapunov exponent as dependent on the parameter $`r`$. These results are obtained for the extrapolated map given by the Eq. (3) and by the method described in Ref.  (Fig. 2). To obtain the Lyapunov exponent we have investigated the subsequent differences between two initially nearby trajectories. The slope of the time dependence of the logarithm of this difference (see Fig. 3) gives the exponent $`\lambda `$. The results for the finite lattice are qualitatively the same (Fig. 2).

In Fig. 4 we show two plots: calculated time dependence of $`A(p)`$ for the lattice of $`50\times 50`$ trees (Fig. 4) and the statistical data on the forest fires in Yukon  (Fig. 4). We have also checked that the phase portraits $`A(p_{n+1})`$ vs. $`A(p_n)`$ are more or less the same as the appropriate phase portraits for the experimental data  (Fig. 5). As we can see, small fires occur more often than large ones. Both kinds of the phase portraits show some anticorrelations between the forest fires in subsequent years: large burned area in a given year is likely to be followed by a small area in the subsequent year. This remains true also for the data on the forest fires in other provinces of Canada in years 1970–2000.

The bifurcation diagram for the map obtained by HKA (shown in Fig. 1) presents itself as a homogeneous spot (see Fig. 6). The same is true for the bifurcation diagram for the analytic map defined by Eq. (3). The histogram of the forest fires for the map obtained by HKA is presented in Fig. 7, together with the histogram for the analytic map. The difference is that in the latter case, small fires do not occur. Still however, the events without fires are the most frequent in all the cases. Both histograms show large amount of events without fires, and the exponential shape in the range of intermediate and large fires. In this region $`N(A)`$ increases with $`A`$, on the contrary to the predictions of the models which display the self-organized criticality . This can be interpreted as a consequence of our assumption, that large clusters are ignited more likely. The obtained curve is shifted to right when the parameter $`r`$ increases, as it is shown in Fig. 7 for the analytic map.

## 4 Conclusions

The results presented allow to state that it is not possible to predict the statistics of the forest fires in long time scale. This conclusion is drawn from our much simplified model of the fires, where the forest is reduced to a lattice of cells, and a tree ignition — to the nearest-neighbor relation with a tree which already burns. Obviously, almost all technical details of a real fire  are omitted in this approach. Still we believe that our conclusion is justified. The argument is as follows: suppose that a much simplified model predicts chaos. Suppose that the model has been improved by introducing several complex details. Is it possible that the new version gives a regular and predictable behavior? The answer is: no. The chaotic character could be eliminated by a noise, but the predictability cannot be improved .

## Acknowledgements

The authors thank to Professor Dietrich Stauffer for valuable comments. The simulations were carried out in ACK-CYFRONET-AGH. The time on SGI 2800 machine is financed by the Polish Committee for Scientific Research (KBN) with grant No. KBN/SGI2800/022/2002.