text:
# The dynamical evolution of the circumstellar gas around low-and intermediate-mass stars I: the AGB

## 1 Introduction

After the completion of hydrogen and helium core burning, low- and-intermediate mass stars ascend the Asymptotic Giant Branch (AGB) in the HR diagram. At the end of the AGB phase, during the so-called thermal pulsing AGB stage (TP-AGB), these stars experience high mass-loss rates modulated by the occurrence of thermal pulses. A circumstellar envelope (CSE) of gas and dust is formed which later on will be ionized by the hot stellar remnant, forming a planetary nebula (PN).

The mass-loss during the AGB plays a decisive role in the subsequent evolution of the star. The star detaches from the AGB when the stellar envelope mass falls below a critical value. Moreover, it is during the AGB phase when heavy elements which have been processed in the stellar interiors are returned to the interstellar medium (ISM). The mass-loss itself also determines the properties of the AGB population, the distribution of white dwarf masses, the maximum white dwarf progenitor mass (hence the minimum supernova type II progenitor mass), and the PN structure.

Mass-loss is a crucial process in the formation of PNe and in the evolution of stars, however, it cannot be calculated from first principles. The origin of the high mass-loss rates (the so-called superwind (SW); Renzini 1981) needed to remove the stellar envelope at the tip of the AGB is still a matter of controversy (Mowlavi 1998). The scenario generally accepted nowadays involves two processes which work together: shock waves caused by the stellar pulsation and the acceleration of dust by radiation pressure. The stellar pulsation in AGB stars (Mira-type stars) creates shock waves which propagate through the stellar atmosphere. The dissipation of the mechanical energy associated with these shocks leads to the levitation of the upper layers of the atmosphere, where the gas becomes sufficiently cool (by expansion and by dilution of the stellar radiation field) and dense to allow heavy elements to condense into grains. As grains nucleate and grow they experience the force exerted by the stellar radiation pressure and thus are accelerated. The momentum coupling between gas and dust drives the outflow. Pioneering work on dynamical calculations to drive the stellar wind were first made by Wood (1979) and Bowen (1988). However, it is possible that other mechanisms such as mass-loss induced by rotation (Dorfi & Höfner 1996), acoustic waves in the stellar atmospheres (Pijpers & Habing 1989), radiation pressure in molecules (Maciel 1976), turbulent pressure (Jiang & Huang 1997), Alfven waves (Hartmann & MacGregor 1980) and magneto-hydrodynamic processes (Pascoli 1994), may also play an important role in the mass-loss process. In addition, large-scale convection and thermal expansion should also need to be taken into account in the mass-loss mechanism.

Dynamical model atmospheres are needed for the calculation of mass-loss. However, problems arise when time-dependent dynamics (shock waves and winds), radiation transfer (strong variable stellar radiation field) and dust and molecular formation processes need to be considered together. Computational limitations and our current knowledge of the fundamental physical data make even the state-of-the-art models far from ideal in several ways. For one, they are restricted to grey radiative transfer. Winters et al. (1994) computed models that give mass-loss rates as a function of the fundamental stellar parameters for stationary atmospheres while dynamical model atmospheres have been computed by Bowen (1988), Fleischer, Gauger & Sedlmayr (1992), Arndt, Fleischer & Sedlmayr (1997) and Höfner et al. (1998). Höfner et al. (1998) demonstrated how changes in the micro-physics, in particular, the use of more realistic opacities in the atmosphere, may result in considerable lower mass-loss rates. The problem is that obtaining reliable quantitative theoretical mass-loss rates from these models is still not possible (Olofsson 1999).

Due to the inherent difficulties in the calculations of mass-loss, and in recovering the history of mass-loss from observational studies, its exact evolution still remains unknown. It is this basically unknown mass-loss history experienced by the star during the AGB which determines the circumstellar structure at the end of the TP-AGB and, moreover, provides an important clue for understanding PNe formation. Stellar evolutionary models are able to provide the temporal behavior of the mass-loss during the AGB and are available in the literature (Vassiliadis & Wood 1993; Blöcker 1995; Schröder et al. 1999). Although the mass-loss rates are not derived from first principles in these models, and most of them rely either on the dynamical model atmosphere calculations of Bowen (1988) and Arndt et al. (1997), or on the semi-empirical mass-loss rates formula derivations of Wood (1990), they do provide a unique opportunity to study the extensive history of mass-loss on the AGB and beyond.

Hydrodynamical models of the evolution of a dusty circumstellar shell in the final stages of the AGB evolution have been performed by Steffen, Szcerba & Schönberner (1998) and Steffen & Schönberner (2000). They used the long term variations of the mass-loss for a 3 M$`_{}`$star during the AGB predicted by Blöcker (1995) and a two component radiation hydrodynamical code to drive the wind. Due to the inherent difficulties of this kind of approach they did not solve the energy equation for the gas component. The grid sizes they used do not allow to study the full mass-loss process experienced by the star, they were limited to $`1\times 10^{18}`$ pc (Steffen et al. 1998) and $`5\times 10^{17}`$ pc (Steffen & Schönberner 2000). The influence of the ISM could not be studied in these models because the ISM is pushed out of the computational domain during the evolution .

Multiple dynamical interactions are expected through the wind which is strongly modulated by the thermal pulses. Our aim is to study the structures generated by the full process of mass-loss during the AGB for different progenitors. For this we have ensured that our grids are large enough to allow us to fully follow the structures generated by the stellar ejecta. Thus the influence of the ISM on the shell formation process can be addressed.

We present a grid of models describing the dynamical gas evolution during the AGB, encompassing the range of initial masses of PN progenitors. We use the wind temporal behavior predicted by Vassiliadis (1992) for stars with solar metallicity. In $`\mathrm{\S }2`$ we describe our numerical method, the initial and boundary conditions adopted for the simulations, and the assumptions on which the models are based. In $`\mathrm{\S }3`$ we describe separately the results obtained for each stellar mass. The discussion of the evolution during the AGB with mass-loss is presented in $`\mathrm{\S }4`$, together with a comparative analysis for the different progenitors, the ISM influence and a limited comparison of our results with the observations.

## 2 The numerical method and computational details

We have performed hydrodynamical simulations of the wind evolution starting at distances from the star where the wind has reached its terminal velocity. For the numerical simulations we have used the multi-purpose fluid solver ZEUS-3D (version 3.4), developed by M.L. Norman and the Laboratory for computational Astrophysics (LCA). This is a finite-difference, Eulerian, fully explicit code which works efficiently in one or two dimensions by choosing the appropriate symmetry axis (for further details about the numerical methods used in the code see Stone & Norman 1992a; Stone & Norman 1992b; Stone, Mihalas & Norman 1992).

We have performed the simulations in 1D using spherical coordinates. The radial sizes of the grids as well as the cell sizes for each model are specified in Table 1. Column (1) gives the (initial) mass of the models. The TP-AGB evolutionary time and the time until the onset of the superwind are given in columns (2) and (3). We have adopted a resolution of 1000 zones in the radial coordinate of the grid for all the models, therefore, the cell size changes depending of the radial extension. We are interested in probing the effects of the mass-loss rate adopted by the stellar models on the circumstellar structure, thus we should follow the whole stellar ejecta. To this end, the radial extensions have been chosen to be large enough to avoid the mass-loss flow out of the grid. Since the wind has reached its terminal velocity at the innermost radial zone in our simulation which lies at 8$`\times 10^{16}`$ $`\mathrm{cm}`$ from the central star, we do not need a two component (dust and gas) hydrodynamical code to drive the wind. For the cooling curves given by Dalgarno & McCray (1972) and MacDonald & Bailey (1981) a gas with solar composition has been used.

### 2.1 The boundary conditions: the stellar evolutionary models

The inner boundary condition for our simulations is determined by the way the stellar wind flows from the star. The gas density at the inner boundary is given by

$$\rho (t)=\frac{\dot{M}(t)}{4\pi r^2v_{\mathrm{}}(t)}$$
(1)
where the temporal dependence of the mass-loss rate (Ṁ(t)) and terminal wind velocity ($`v_{\mathrm{}}`$(t)) is given by the stellar evolutionary models of Vassiliadis (1992) for stars with solar composition. We have chosen these models because they provide with the temporal evolution of Ṁ and $`v_{\mathrm{}}`$ during the AGB which is needed to study the long-term dynamical evolution of the gas. The parameterization made for the mass-loss and the wind terminal velocity is based on semi-empirical relations.

Vassiliadis & Wood (1993) (hereafter VW93) computed the evolution of stars from the Main Sequence through the Red Giant phase, to the end of the AGB. In their formulation, the mass-loss during the AGB is derived from an empirical formula computed from observations of Mira variables and pulsating OH/IR stars in the Galaxy and the Large Magellanic Cloud. Based on the empirical relation between the stellar pulsation and the period given by Wood (1990), they established that the mass-loss rate is a function of the pulsational period. They then derive Ṁ from fundamental physical parameters such as the mass and the radius of the star by assuming that Miras, and the AGB stars surrounded by dust are pulsating in the fundamental mode. They consider two distinct phases for mass-loss: for stellar masses below 2.5 M, and for periods less than $``$500 days the mass-loss increases exponentially with the period (Mira phase), and for stars with mass greater than 2.5 M$`_{}`$they delay the onset of the superwind to periods of $``$ 750 days. For periods beyond $``$ 500 days, the mass-loss rate is essentially constant, and lies within a factor of $``$ 2 of the value given by the radiation pressure limit. The terminal wind velocity ($`v_{\mathrm{}}`$) is computed assuming the empirical relation for the pulsation period derived by Wood (1990), and it is restricted to lie in the range 3.0-15.0 km s$`^1`$.

We take the temporal evolution of the wind parameters for stellar masses 1, 1.5, 2, 2.5, 3.5, and 5 M$`_{}`$with solar metallicity (Z = 0.016) from the stellar evolutionary models computed by Vassiliadis (1992), and set them in the five innermost zones of the grid for the AGB evolution of the star. In Fig. 1 we show the evolution of the mass-loss and terminal wind velocity used as the input for our models. We have labeled the main wind modulations (e.g I, II) in order to allow a better description through the text of the features they generate in the gas. The wind temperature is assumed to be the effective temperature of the star. A free-streaming boundary condition was set at the outer radial direction.

Fig. 1 shows that both Ṁ and $`v_{\mathrm{}}`$ suffer from strong modulations. The wind density is a bimodal function of these parameters, and therefore whenever one changes so does the density. The density will be softened when Ṁ and $`v_{\mathrm{}}`$ change in such a way that they neutralize each others effect.

### 2.2 The initial conditions

The mass-loss experienced by low-and intermediate-mass stars during the main sequence and red giant phases ($`10^8`$ M yr-1) is negligible compare to the mass-loss during the AGB (up to $`10^4`$ M yr-1). The local ISM is not significantly modified by the wind of previous evolutionary stages, due to the low mass-loss rates and velocities involved. Therefore, in our study, we have considered an homogeneous ISM as the initial condition.

The main component of the ISM is hydrogen in both molecular and atomic form, and dust, although the amount of mass residing in the form of dust is much smaller than that residing in gas (estimates of the dust-to gas mass ratio are $`3\times 10^3`$; Burton, Elmegreen & Genzel 1992). Although molecular hydrogen dominates over all gaseous material in the temperature range from a few to about 30 K, it is confined to clouds (small, massive and compact regions) which have an annular distribution in the inner Galaxy (from 3 to 8 kpc). The diffuse ISM is thought to consist of four major components: the cold neutral medium (CNM), the warm neutral medium (WNM), the warm ionized medium (WIM), and the hot ionized medium (HIM). Neutral atomic hydrogen (HI) is the main observed constituent of the ISM; estimates of its filling factor ranges from 20% to 90% (Burton 1988). The WNM is the ISM component which occupies most of the space (Kulkarni & Heiles 1988) in an “intercloud medium” or surrounding the cold clouds. Since information about the WNM mainly comes from 21 cm line emission data, direct measurements of temperatures (estimated to be in the range 5000 K to 8000 K) or densities (estimated to be between 0.1 and 1 $`cm^3`$) are difficult to obtain. Therefore, the WNM is the least well understood component of the ISM. A typical density of 0.4 $`cm^3`$ characterizes much of the ISM (Burton 1988).

We assume that the ISM has the characteristics of the WNM because it is the main observed constituent of the ISM. We have considered that the ISM pressure is simply the standard gas kinetic pressure, $`P=nKT`$ (where $`n`$ is the number of particles per unit volume, $`K`$ is the Boltzman’s constant and $`T`$ is the gas temperature). Spitzer (1978) demonstrated that the contribution of cosmic rays, and magnetic pressures can be comparable to the gas kinetic pressure. As we have not included magnetic fields or turbulent motions in our simulations, we therefore, use large densities (1 $`cm^3`$ instead of the typical 0.4 $`cm^3`$) to compensate for the cosmic rays, magnetic ($`B^2/8\pi `$) and turbulent pressure components (assuming all of them contribute to an isotropic pressure).

In our models, at time zero the grid is filled homogeneously with neutral atomic hydrogen with a density of 1 $`\mathrm{cm}^3`$, a temperature of 6000 K, and zero macroscopic velocity. Hence, we have to consider wind propagation in an homogeneous (cloudless), stationary ($`v=0`$) and warm ($`T=6000\mathrm{K}`$) ISM (the isothermal sound speed is $`c_\mathrm{s}=7`$ km s$`^1`$). The interaction of the stellar wind with the ISM influences the subsequent evolution of the flow. The adopted ISM pressure of 6000 $`\mathrm{K}\mathrm{cm}^3`$ lies in the range of the ISM pressure which allows pressure equilibrium between the different ISM components. In particular, thermal equilibrium between the WNM and the CNM is only possible over a range of pressure between $``$2500 and $``$6500 (in units of the the Boltzman’s constant). The gas evolution considering an ISM with pressure ten times lower is also presented in $`\mathrm{\S }4`$.

In order to keep the ISM at its original temperature, no source of gas heating is included. The temperature of the ISM decreases due to radiative cooling during the evolution a negligible amount ($``$ 500 K).

## 3 Results

The beginning of the evolution is the same for all the stellar models. A wind with a constant mass-loss rate ($`10^8`$ M yr-1) and velocity ($``$3 km s$`^1`$) starts to interact with the ISM. The duration of this stage of constant mass-loss rate at constant velocity depends on the stellar mass (see third column in Table 1)111We considered that the so-called superwind phase starts when the mass-loss increases to values greater than $`2\times 10^7`$ M yr-1.. The result is a compression wave propagating through the ISM at the speed of sound (7 km s$`^1`$for $`T_{ISM}=6000\mathrm{K}`$). It should be noted that since the wind temperature is assumed to be the $`T_{eff}`$ of the star ($``$ 2500 K) a supersonic interaction is expect to occur for lower velocity values. Therefore, despite the low values reached by the wind parameters at the wind modulations labeled as I in Fig. 1 for the 1 M$`_{}`$model, and at I and II in the 1.5, 2 and 2.5 M$`_{}`$models, all of them are able to develop shock waves initially. The subsequent evolution of the flow is related to the wind history of each stellar model and are described below.

The minimum and maximum values reached by the wind temperature during the evolution are 2630 K and 3160 K respectively. The relative importance of the thermal pressure of the wind in the shell morphology, although small at the end of the evolution still affects the final size of the shells.

In the following, we present the results of the hydrodynamical simulations of the evolution of the stellar wind ejected during the TP-AGB for each stellar mass considered. We have selected the most representative outputs in order to describe the evolution, avoiding the transient features appearing in the flow that rather complicates the description.

### 3.1 The 1 M$`_{}`$star case

The dynamical evolution of the gas is described by relating it with the evolution of the central star wind. As explained in $`\mathrm{\S }2`$, the time-dependent values of the mass-loss rate, velocity and wind temperature are used as the inner boundary condition (see top of Fig. 1 for the 1 M$`_{}`$values). The mass-loss reaches the radiation pressure limit (i.e the maximum mass-loss value) only during the last four thermal pulses. Between the pulses the mass-loss rate and velocity decreases, due to the extended luminosity dips which characterize the thermal pulses for low-mass stars.

In Fig. 2, we have summarized the AGB evolution of the flow by plotting five radial density and velocity profiles at different stages during the evolution. From top to bottom, the plots have been selected at 78, 84, 90, 97, and 99 per cent of the total AGB time, which in Fig. 1 corresponds to the gap between pulse II and III, to the beginning of III, to 20,000 yr before the end of III, to the gap between III and IV, and to the middle of IV respectively. We have also marked each density peak and its corresponding velocity.

The main major changes in the gas structure appear as a consequence of the wind modulations associated with pulse II in Fig. 1. The increase in mass-loss rate and velocity during this pulse causes the pile up of material identified as feature a222Hereafter, shell features are singly called with the letters they are referred to in the figures., (i.e. a shock structure). At the end of II, the wind velocity drops to nearly half of its value during the pulse, while the mass-loss rate remains constant. A drop in velocity at constant mass-loss rate gives rise to an increase in density in the inner parts of the grid. Denser material from the inner regions moves outwards at lower velocities than the material ejected previously, hence forming a gap in the density structure. The shell identified as b is been formed by this process. When the high mass-loss rate during pulse II ceases, c is formed. In the top panel of Fig. 2 we show the density structure $``$ 7,400 yr after the formation of c. Both b and c are accelerated by pressure gradients, whilst a is decelerated by the ISM.

About 32,000 yr later (second panel from top) only one density peak, ab, formed when b reaches a and they merge, is visible. Due to velocity gradients c has broadened and has disappeared. At the time shown, pulse III has already started, and therefore the density and velocity at the inner parts of the grid have increased (see also the third plot from the top 32,000 yr later). Due to pulse III, denser material moving at a higher velocity reaches the previously ejected material. A shock develops in the region of interaction between the mass lost during pulse II and that occurring during III. This shock region will give rise to the formation of the shell identified as x (see fourth plot from top). Shells d and e are formed by the same process which formed shells b and c, i.e. a gradual velocity drop at the end of III (forming d) and the subsequent drop in mass-loss rate which forms e.

The gas structure 7,600 yr later is shown at the bottom panels of Fig. 2. In this case, it is d, formed by the velocity drop at the end of the pulse, which broadens and disappears, whereas e last for a longer time because it is not captured by denser material leaving the star. The gas density reaches higher values during pulse IV due to the small increase in the wind velocity.

In Fig. 3 we show the temporal behavior in terms of radius, velocity and density of the shells formed in the CSE during the TP-AGB evolution for a star with 1 M$`_{}`$. Note that only shell x, formed as a consequence of a shock, and shell a, due to the capture of b, increase in density. Note that all the shells are decelerated during the evolution.

### 3.2 The 1.5 M$`_{}`$star case

The temporal behavior of the wind parameters used as the inner boundary condition for the study of the 1.5M$`_{}`$model are shown in the second panel of Fig. 1. The circumstellar gas evolution during the TP-AGB is plotted in Fig. 4. The five outputs shown in Fig. 4 have been selected at times which corresponds in Fig. 1 to 12,300 yr after the end of III, 95,600 yr since the begining of IV, 15,000 yr before and after the end of IV and 43,000 yr since the beginning of V. The percentage of the total AGB time-life at which the density and velocity profiles are shown is specified in the upper right corner of the density plots.

Despite the different wind history experienced by stars with 1 and 1.5 $`M_{}`$, the general behavior of the circumstellar gas is very similar. The density peaks a, and b seen in the top panel of Fig. 4 are formed by the mass accumulated during III, and by the reduction of mass-loss and velocity in the aftermath of this pulse respectively. Pressure gradients are responsible for the acceleration of b, while a is decelerated by the pressure exerted by the ISM gas which lies ahead. Later on ($``$30,000 yr), these two shells merge and ab is then formed.

The density and velocity profiles in the third plot of Fig. 4 have been disrupted by the higher mass-loss and velocity wind arising from the star, i.e. pulse IV has begun. As in the 1 M$`_{}`$ case, a shock develops at the position where the mass lost by IV encounters the previously ejected material which has been slowed down during the evolution. Later on, x will appear at this position. Shell c is formed when IV ends $`30,000`$ yr later (fourth plot from top). This shell is smoothed rapidly due the large velocity gradients. The transient feature labeled d is formed by the mass-loss rate and velocity increases at V. We can also see how, as a consequence of the shock, matter in x is compressed.

The shell x1, visible in the bottom plot of Fig. 4, is formed in the interaction region between the material which has been ejected during pulse V and the one ejected during pulse IV.

### 3.3 The 2 M$`_{}`$star case

The evolution of the mass lost by a star with a main sequence mass of 2 M$`_{}`$is shown at ten different times in Figs. 5 and 6. The gas density and velocity radial profiles have been selected at times which correspond in Fig. 1 to 23,000 yr after the beginning of II, 4,000 yr after the end of II, 16,000 yr after the end of II, 3,000 yr after the end of III and 6,000 yr after the beginning of IV.

Both a and b in Fig. 5 are caused by II, the former by the pile up of material during this pulse, and the latter when the pulse has finished. They merge $``$ 18,000 yr later, forming ab. Shell ab propagates outwards, and finally disappears only slightly modifying the surrounding gas. Shells c and d have the same formation process as shells a and b mentioned previously, but in this case their origin is associated with III. Again, because c is being decelerated, d catches it up in $`3\times 10^4`$ yr and hence cd appears.

In Fig. 6 we have shown the gas density and velocity evolution at times which correspond in Fig. 1 to 4,500 and 3,200 yr before and after the end of IV respectively and 4,500, 16,500 and 71,300 yr after the beginning of V. At the top of Fig. 6 we have marked the density peak e which is formed by IV. A strong shock develops in the density valley formed between cd and e, leading to the formation of x. When pulse IV finishes, shell f appears. With increasing time, all the density peaks (cd, e, f) disappear because they feed x. Only x remains as an identifiable shell until the end of the AGB. Pulse V is responsible for the increase labeled as g which completely disappears 4,400 yr after when the density falls as a power law of the type $`r^2`$ law. This is what is expected for a constant mass-loss rate at a steady wind velocity. The disruption of the density structure seen in the bottom plot of Fig. 6 is caused by the velocity decrease produced during the middle of pulse V.

### 3.4 The 2.5 M$`_{}`$star case

Of all the cases considered this model has the longest evolution during the AGB. The wind parameters undergo three main modulations in mass-loss rate and velocity, which occur during the last 300,000 yr (see Fig. 1). The evolution of the gas is displayed at five snapshots during the TP-AGB in Fig. 7. The times selected have a relation to the stellar evolution as follows: the top plot corresponds to 6,500 yr after the end of II, the second and third plots have been selected at 39,000 and 8,000 yr after the begining and end of III respectively, and the fourth and fifth plots at 45,500 and 12,300 yr after the begining of IV.

The small modulations in velocity experienced by the wind before pulse II do not have any effect on the gas structure. The first significant effect arises as a consequence of II when a high density shell is formed. Due to its high thermal pressure this shell expands splitting the density and velocity into two peaks. This peculiar two peak structure lasts for around 22,000 yr, disappearing as the gas propagates outwards, resulting in the formation of a (see top plot in Fig. 7). Shells b and b1 are formed by the gradual end of II: first, both the wind velocity and the mass-loss rate are reduced by factors of 3 and 5 respectively (b is formed), and second, the mass-loss rate drops again a factor of 1.5, meanwhile the wind velocity grows in III and b1 is formed. Shell b1 survives for only a short time (8,000 yr), it is accelerated by the incoming wind and quickly merges with b.

A density peak around 0.1 pc appears in the second plot of Fig. 7 due to the fast increase in mass-loss and velocity experienced by the wind towards the end of III, c is formed when III ends. A shock where the gas is compressed (x in the fourth plot) is produced in the region of interaction between structures ab and c.

All the modulations produced during IV give rise to the formation of transient density peaks, such as d, which are smoothed out during the evolution of the flow. Only x remains at the end of this stage (see bottom plot of Fig. 7).

### 3.5 The 3.5 M$`_{}`$star case

Both the wind velocity and mass-loss rate for the 3.5 M$`_{}`$model grow continuously during the TP-AGB, suffering only from small sudden drops in intensity (see Fig. 1).

The gas evolution is shown in Fig. 8, where the plots have been selected in order to show the formation of structures that appear after the three major changes in the wind parameters. Shell a is formed at 220,000 yr when the mass-loss rate reaches a value of $`5\times 10^7`$ M yr-1and the wind velocity is 8 km s$`^1`$. This shell is fed continuously by the stellar wind. There are three important drops in the wind velocity and mass-loss rate that give rise to the formation of the transient features labeled b, c and d respectively. Only one shell (abcd) is visible in the bottom plot of Fig. 8 which is the result of the accumulation of matter during the AGB evolution of the star.

### 3.6 The 5 M$`_{}`$star case

The relatively simple behavior of the input wind parameters for the star with 5 M$`_{}`$(without abrupt changes in mass-loss or velocity, see Fig. 1) is translated directly into a simple evolution of the circumstellar gas. As expected, a $`\rho r^2`$ law characterizes the radial density profile of the gas up to the position of an external shell. This cold external shell is formed and fed with the high mass-loss rates from the star and propagates outwards with little disruption. In Fig. 9 the evolution of this shell since its formation is shown at five different times. The top plot in Fig. 9 corresponds to 80,000 yr after the beginning of the AGB evolution and the following plots have been chosen at time intervals of $`1\times 10^4`$, $`5.4\times 10^4`$, $`5.5\times 10^4`$, and $`5.2\times 10^4`$ yr respectively.

## 4 Discussion

In the previous section, we summarized the circumstellar gas evolution during the TP-AGB phase and related it with the stellar wind behavior. The dynamical evolution of the circumstellar gas is highly non-linear, and although the stellar wind on the TP-AGB experiences huge modulations, the circumstellar gas structure does not retain this information. The subsequently formed shells are either accelerated during the evolution, or decelerated by the ISM, most of them disappearing when they merge with other shells or when they reach pressure equilibrium with the surroundings. The shells formed by a shock in the interaction region between two subsequents enhancements of mass-loss and velocity are the most stable structures, and remain visible at the end of the AGB phase.

In table 2 we have summarized, for each of the shells marked in Figs. 2, 4, 5, 6, 7, 8 and 9 the radius (determined at the position of the density peak), density and velocity for the times shown in the figures. The formation process of the shell and its life-time ($`L_t`$) are also given, where “$`\mathrm{}`$” means that the shell remains until the end of the TP-AGB evolution. The largest radius corresponds to the shells with the longest lifetime. It should be noted that although shells ab for both the 1 and 1.5 M$`_{}`$models have the largest radii, they are not the main observable feature. Shell x is the main observable feature due to its highest density. The radii show a range between $``$0.07 pc, for shells that have just been formed, up to $``$1.8 pc for evolved shells at the end of the AGB. The shells velocities range between 2 km s$`^1`$up to 19 km s$`^1`$. Although the maximum wind velocity was 15 km s$`^1`$we can see how many of these shells have been accelerated by thermal pressure.

In Fig. 10 we show the circumstellar gas density and velocity structure at the end of the AGB for all the initial masses considered. For the 1 M$`_{}`$model, the bottom panel of Fig. 2 was selected very close to the end of the AGB stage, and therefore the gas radial profiles are very similar to the ones shown in Fig. 10. In the case of the 1.5 M$`_{}`$ model only one external thick shell remains visible which is formed when shell x merges with shell ab. A double-peaked structure can be seen in the outermost shell of the 2 M$`_{}`$model where the outermost peak correspond to shell x (see Fig. 6) and the innermost peak is produced by a shock region which develops when the matter ejected later on from the star reaches shell x. Therefore this double-peaked structure is not caused by an expansion due to a pressure difference. The density structure of the rest of the models has not changed significantly since the times shown in Figs. 7, 8, and 9.

It has been widely assumed in the literature that a law of the type $`\rho r^2`$ describes the gas structure resulted from the AGB evolution (e.g. Okorokov et al. 1985; Schmidt-Voigt & Köppen 1987a, 1987b; Marten & Schönberner 1991; Mellema 1994; Frank & Mellema 1994). In order to compare our density structures with a constant mass-loss rate at constant velocity density law we have overplotted in Fig. 10 a density law of the form $`\mathrm{log}\rho =2\mathrm{log}(r)+C`$ where $`C=\mathrm{log}(\dot{M}/4\pi v_{\mathrm{}})`$. To compute $`C`$ we have used the values of Ṁ and $`v_{\mathrm{}}`$ reached by the wind at the end of the AGB. The 1 M$`_{}`$model is different to the others in the sense that at the end of its evolution, when pulse IV occurs, the velocity remains very low (in contrast to what happens during previous pulses). This results in a huge increase in density in the inner parts of the grid (see Fig. 10). Two fits have been made for this model; the dotted line represents the fit for the real values of Ṁ and $`v_{\mathrm{}}`$ reached by the wind at the end of AGB phase, the dashed-dotted line is for a wind with a terminal velocity of 15 km s$`^1`$and the Ṁ experienced by the star. The values for the constant $`C`$ adopted to fit the constant mass-loss rate, constant velocity law are plotted in Fig. 11. The relationship of $`C`$ with stellar mass can be described by a linear fit of the form $`C=23.79+0.16M(M_{})`$ (for $`M>1`$ M$`_{}`$), which is shown in Fig. 11. The same values of $`C`$ which fit the density profiles are obtained when other ISM densities are used in the simulations.

The common factor in all these models, with the exception of the 1 M, is that the complicated wind history experienced during the AGB results in the formation of simple structures. The density distribution can be described by a $`r^2`$ density law except for the outermost shell. The gas structure produced by the evolution does not retain the information about the wind history experienced by the stars. The low-mass stars undergo four or five significant modulations in the wind parameters during the AGB. However, except for the 1 M$`_{}`$star model, only one shell is visible in the density structure. Therefore, according to our models, is not possible to recover the history of mass-loss experienced by the star from the density structure at the end of the AGB. Moreover, the prominent peaks are not signatures of mass-loss rate enhancements but are the consequence of the development of shocks in the interaction region between two important mass-loss events. Another important fact is that the shells disperse on timescales of the order of 20,000 yr (much shorter than the characteristic AGB evolutionary time), independent of the formation process. The only exceptions are the shells formed by shocks x which remain as stable structures in the density profiles at the end of the AGB.

According to the formation process of the outer shell we can classify the models into two mass groups. Firstly, the group formed by the stellar masses 1, 1.5, 2 and 2.5 M$`_{}`$for which the outer shell is formed by a shock in the interaction region between two subsequents mass-loss rate events, and secondly, the group formed by stellar masses 3.5 and 5 M, for which the external shell is formed by the accumulation of mass during the course of the evolution. Thus, low- and intermediate-mass stars have different outer shell formation processes.

The stellar models with masses 2, 2.5, 3.5 and 5 M$`_{}`$share similar radial velocities profiles. All of them are characterized by a continuous increase (with a small slope) up to the position of the external shell, where the velocity reaches its highest value of $`20`$ km s$`^1`$. Moving outwards, the sharp fall off in the velocity is interrupted by the presence of a small plateau at the position of the outermost density peak (i.e. external shell).

### 4.1 The role of the ISM

A novel aspect that makes our study different from previous ones (e.g. Steffen et al. 1998, Steffen & Schönberner 2000) is the use of grids large enough to prevent the flowing out of stellar ejecta. Thus, we can quantify the effect of the ISM on the gas structure, since we have ensured that the ISM is not completely pushed out of the computational domain. When the wind flows from the star it encounters a non-negligible pressure provided by the ISM. The ISM is pushed and compressed while the expansion of the stellar wind is slowed down. However, we should keep in mind that the pressure exerted by the ISM can change radically depending on the position of the star within the Galaxy. In particular, the ISM density fall off exponentially in the direction perpendicular to the Galactic plane from a value of $`2\mathrm{cm}^3`$ (Spitzer 1978) and a scale height of 100 pc (Mihalas & Binney 1981).

In order to consider the evolution of stars at higher Galactic latitude, we have also performed the simulations using an ISM ten times less dense than the one used in the models described in the previous sections. This gives us an external pressure ten times smaller, because only the thermal pressure component is considered in this study. The gas structure at the end of the AGB as a result of the evolution through this low density ISM is shown in Fig. 12 for all the initial stellar masses considered. The dotted line represents a fit using a $`r^2`$ density law. A first glance at Fig. 12 reveals that the external shell has lost density contrast compared to the models evolving though a medium with higher pressure.

#### 4.1.1 The mass of the circumstellar envelope

The dynamical interaction of the stellar wind with its environment is important to understand the mixing of the ejecta with the ISM (fundamental for chemical abundance studies). As the stellar ejecta evolves it sweeps up ISM material that has not been enriched in the stellar interiors. In order to address the relative importance of this mixing effect we compute the fraction of the mass contained in the circumstellar shell that belongs to the ISM. This can be done in our case since we know the amount of mass returned to the ISM by the star. We integrate the total amount of mass in the grid which is above the ISM density at the end of the AGB for each stellar mass. This is the mass we might expect to observe above the background level. This is done for the two ISM densities considered and is shown as a function of the stellar mass in the left panel of Fig. 13. The long-dashed and dotted lines represent the integrated mass for ISM densities of 1 and 0.1 $`cm^3`$ respectively. The solid line represents the total mass lost by the star during the course of the AGB evolution. The integrated mass in Fig. 13 approaches the total mass lost by the star when the density of the ISM decreases.

We compute the fraction of the mass residing in the CSE that belongs to the ISM by subtracting from the integrated mass above the ISM level the mass lost by the star, and then dividing by the initial stellar mass. This is shown in the right panel of Fig. 13 for the two ISM density values considered. As expected, the percentage of ISM material in the CSE decreases with the ISM density. However, even for a low density medium this fraction is not a negligible amount. For a star with 1 M$`_{}`$evolving through a medium with density of 1 $`\mathrm{cm}^3`$, approximately 70 % of the matter in the CSE belongs to the ISM. This amount decreases to $``$ 20 % when the medium is 10 times less dense. The swept up ISM matter also decreases with increasing main sequence stellar mass. If a low mass progenitor evolves in a high density medium, up to $`70\%`$ of the measured mass is ISM matter which has not been processed in the stellar interior. This can have important consequences for abundance analysis, since the observed gas is mixed with non chemically enriched ISM matter.

The presence of this non negligible amount of matter from the ISM in the CSE leads to another important problem. The mass of the progenitor derived from observations of the CSE is based on the assumption that all the observed matter has been lost by the star. We have seen that the mass in the CSE increases due to the swept up ISM material, and this effect should be taken into account in order to not overestimate the mass of the progenitor from observations of CSE.

#### 4.1.2 The ISM pressure

The size of the CSE reached at the end of the TP-AGB depends on the external pressure, the momentum of the wind, and on the time they are allowed to evolve. In order to quantify the relative importance of all these quantities, we have plotted in the left panel of Fig. 14 the radius of the CSE at the end of the TP-AGB as a function of the initial stellar mass for the two ISM densities considered ($`1\mathrm{cm}^3`$ long-dashed line and $`0.1\mathrm{cm}^3`$ dotted line). The radius has been determined as the largest distance from the star where the ISM has been modified, we can see that the radius increases by a factor of $`1.5`$ when the ISM density decreases by a factor of ten. In the right panel of Fig. 14, we show the TP-AGB evolutionary time versus the stellar mass (solid line), and the kinematical ages of the envelope for the two cases considered of ISM densities, $`1\mathrm{cm}^3`$ and $`0.1\mathrm{cm}^3`$ (long-dashed line and dotted line respectively). The kinematic age of the envelope has been computed by using the radii plotted in the left panel of Fig. 14 and assuming a constant outflow velocity of 15 km s$`^1`$(a common value usually reported in the literature; i.e. Young, Phillips & Knapp 1993a).

We can see from the right panel of Fig. 14 that the kinematic ages computed are almost the same for all the stellar masses, even when the differences in the evolutionary times are as large as a factor of ten. A possible consequence of this results is that the CSE kinematical ages are not ideal to constrain the evolutionary age of the star, thus the progenitor mass. The largest radius is reached by the star with initial mass 2.5 M, i.e. the star that has the longest evolutionary time. The secondary maximum in radius (for the 1.5 M$`_{}`$model) is related with the total time the wind is supersonic. Pressure equilibrium with the ISM is not reached in any of the two cases of ISM pressures considered.

In Table 3 we have listed the sizes, densities and velocities of the circumstellar shells at the end of the AGB, where cases A and B represent the two ISM densities considered, 1 and 0.1 $`\mathrm{cm}^3`$ respectively. These values have been determined at the position of the outermost density maximum, and are aimed to illustrate the effect of the ISM pressure on the CSE properties. Both values considered for the thermal pressures of the ISM, 6,000 and 600 $`\mathrm{cm}^3`$ K (where K is the Boltzman’s constant) are able to confine the shells (note the density enhancement at the border of the density profiles). An ISM with a thermal pressure ten times smaller (case B in Table 3) results in external shells with average radii 1.1 times larger, shell densities 0.2 times smaller and expansion velocities which are 1.5 times faster. Moreover, the kinematical ages derived for the ISM with less pressure are on average 0.8 times smaller.

The net effect the ISM has on the formation of the CSE is to decrease its size, stall its expansion, and increase its mass. Moreover, and even more important, a non-negligible amount of ISM material is mixed in the CSE with the chemical enriched material ejected from the star. Thus, the effect of the ISM is very important and should not be neglected in further studies.

### 4.2 Comparison with the observations

The CSE can be observed at infrared wavelengths, where the emission comes mainly from the heated dust, and in the radio where it is produced mainly by molecules such us CO and OH. Because the molecular line emission arises from photodissociation, its extent is limited to $`10^{16}`$ cm for thin shells, although it can be larger by up to few times $`10^{17}`$ cm for thicker envelopes (Mamon et al. 1988; Letzelter et al. 1987). Therefore, using observations of the molecular emission from CSE, one can only probe the recent mass-loss history experienced by the star (Olofsson et al. 1993). From an observational point of view, tracing the mass-loss history as far back as $`10^5`$ yr needs to be done at infrared wavelengths, but the dust has low emissivity at the temperatures and distances involved.

According to our models, large CSE are formed around AGB stars, with shells sizes ranging from 1.7 to 2.6 pc (see Table 3) depending on the progenitor mass and on the ISM pressure.

A large CSE surrounding the oxygen-rich ‘red giant’ W Hydrae was discovered with IRAS observations by Hawkins (1990). Emission was detected out to 40  (up to 35 at the 6 sigma level) which gives a diameter larger than 1.6 pc, assuming a distance of 135 pc. W Hydrae is one of the nearest red giant star and it is found at high galactic latitude in an area of low infrared cirrus emission that might otherwise confuse its surrounding extended emission. Gillett et al. (1986) reported the observation of another large CSE around R Coronae borealis. It has an extended shell with dimensions $`r_{\mathrm{min}}=0.65`$ pc and $`r_{\mathrm{max}}=4.3`$ pc (assuming a distance of 1.6 kpc). The galactic latitude of R CrB is 51 , which places it at about 1.24 kpc above the Galactic plane.

Of the 76 CSE resolved in the 60 $`\mu m`$ band in the IRAS survey of Young, Phillips & Knapp 1993b), 26 have outer radii larger than 0.7 pc, and 12 of these have radius between 1.2 and 2.5 pc, in excellent agreement with our predictions. The shell average radius for the whole Young et al. (1993a) sample is 0.74 pc. Three carbon stars W Pic, RY Dra and R CrB were found to have shells with radii of 3 pc or more.

Recently, Speck, Meixner & Knapp (2000) have found highly extended dust shells around two well known proto-PN, AFGL 2688 (the Egg nebula), and AFGL 618. In both cases the dust emission extends out to a radius of 300-400 ′′, whilst the previously reported radius of the optical reflection nebula is at least 10 times smaller. The size of the molecular emission region is of the same order of magnitude as that observed in the optical. The circumstellar shells emit in the far-infrared because the emission is mainly due to cool dust ($``$ 30-50 K). In the case of the egg nebula, the far-infrared ISOPHOT emission is matched by a uniform extended source of dust with size $``$ 2 pc and two denser dust shells at radii $``$ 0.87 and $``$ 1.7 pc (with the adopted distance of 1.2 Kpc). The emission of AFGL 618 is matched by a weak extended source with radius $``$ 3.3 pc and two enhancements in emission at radii $``$ 1.3 and $``$ 2.3 pc (with the adopted distance of 1.7 Kpc). The sizes of the observed shells are in agreement with the size predictions of our models (see Table 3).

All of these objects have spherical symmetric shells which implies that no source of asymmetry is needed in the mass-loss rate process from the AGB, and confirms our predictions that large CSE should be present around AGB stars. Probably, the existence of these huge dust shells is related to the mass-loss suffered by the stars during the TP-AGB. The formation of structures with sizes as large as $``$ 3.3 pc, like those found around AFGL 618, W Pic, RY Dra or R CrB, is not predicted by any of our models, however the distance errors and a lower density for the ISM (those are objects at high Galactic latitude), should be taken into account.

With regards to the emission enhancements in the PPNe AFGL 618 and AFGL 2688, Speck et al. (2000) related their existence to the occurrence of thermal pulses. By comparing the dynamical timescales of the shells with the thermal pulses timescale they derived the mass of the star. This kind of analysis is widely done in the literature. However, the non-linear evolution of the shells generated by the thermal pulses, and the deceleration introduced by the ISM cause that the models with different evolutionary AGB give rise to the formation of structures with similar sizes. According to our models, does not seem possible to relate the dynamical timescales derived from CSE with the theoretical timescales associated with the thermal pulses.

The only known CSE surrounding an AGB star which shows multiple shell structures is IRC+10216. IRC+10216, is the CSE of the long period Mira-type variable carbon star CW Leo has been widely observed at infrared and millimeter wavelengths. Its distance is estimated to be in the range 120 (Loup et al. 1993) to 170 pc (Winters et al. 1994) and is considered to be at the end of its AGB evolution with a mass-loss rate of $`1.5\times 10^5`$ M yr-1(Huggins 1995) and a terminal wind velocity of $``$ 14-15 km s$`^1`$ (Olofsson et al. 1993). It is the nearest carbon star that has a high mass-loss rate and consequently is one of the best observed AGB stars. Recently, deep B and V-band and archival HST WFPC2 images revealed the existence of an extended circular halo with multiple shell structures in the extended CSE IRC+10216 (Mauron & Huggins 2000). The envelope is detected out to 200 ′′. The brightness of the shells is explained in terms of the illumination of dust-scattered ambient Galactic light. The observed intensity is proportional to the column density along any line of sight for an optically thin limit, when the source of illumination is the interstellar radiation field.

We have computed the optical intensity from our models assuming that it is proportional to the column density along the line of sight. We show two examples that could match the observed structure found in IRC+10216 (see Fig. 7 from Mauron & Huggins 2000). The left panel of Fig. 15 corresponds to the 1 M$`_{}`$stellar models at 97% of its AGB evolution and the right panel to the 1.5 M$`_{}`$stellar model at its 93% of its AGB evolution. These are the only two models that at the end of their AGB evolution have enough density peaks for account to the ones observed in IRC+10216 (for a direct comparison see Fig. 15 and Fig. 7 from Mauron & Huggins 2000). The central star in both models suffers from mass-loss rates consistent with the rates observed in IRC+10216 ($`0.35\times 10^5`$ and $`0.8\times 10^5`$ M yr-1respectively). For the 1.5 M$`_{}`$model the terminal wind velocity is exactly what is observed, but for the 1 M$`_{}`$star model it is a factor of two lower. By confronting the predictions of nucleosynthesis models with the measurement of isotopic ratios in the CSE of CW Leo, Kahane et al. (2000) infered a low progenitor mass for the star, M $``$ 2 M. In the last few years models have appeared in the literature to explain the origin of multiple shell events in CSE or post-AGB stars. Some of them invoke the presence of a binary companion with an eccentric orbit (Harpaz, Rappaport & Soker 1997) or a detached binary companion (Mastrodemos & Morris 1999). All the models involving binary companions predicts strict regularity in the arc spacing between the shells which are not observed. The presence of a binary interaction together with a dynamo mechanism in the AGB star, such as that experienced by the Sun has been proposed by Soker (2000).

The shells observed in IRC+10216 shows trends such as increases in shell thickness with increasing radius, and a wide range of shell spacing, which cannot be accounted by the current models proposed for the formation of discrete shells in CSE envelopes. Although, we do not reproduce either the number of shells or the global size of the CSE envelope, the properties of the shells, such as the thickness increase with increasing radius, and the wide range of shell spacing are well characterized by our models. Evidence for shell interactions have been also pointed out by Mauron & Huggins (2000).

The shells found around IRC+10216 do not show many of the characteristics found in the shells around PPN and PNe. Moreover, by the end of the AGB evolution our models predict that the shells disappear. Therefore, if our proposed mechanism for IRC+10216 is correct, it cannot be the same one that forms the shells observed around PPN. These arcs found around proto-PN must have a more recent origin (see García-Segura et al. 2001).

### 4.3 The imprints of stellar evolution in the shells

Since we know the evolution of the wind from the star, our aim now is to investigate what are the characteristics that remain recorded in the final observable structure.

As an example we have made a simple test. We have chosen the 1 M$`_{}`$model at 3.8$`\times 10^5`$ yr since the begining of the evolution. The density structure consists of three well defined shells. Let’s now examine the processes by which these three shells are formed. The outer shell is formed at 1.8$`\times 10^5`$ yr and is continuously fed by the mass-loss from the star during II. The intermediate shell is caused by the decrease in velocity experienced at the end of II, while the inner shell is produced later by the decrease in the mass-loss associated with the end of the same pulse (see Fig. 1). The difference in time between the formation of the inner and the intermediate shell is 4,000 yr while between the intermediate and the outer it is 190,600 yr.

Let’s now assume that we are able to observe this density structure and try to infer the mass-loss modulation which gives rise to its formation. In principle we should measure shell separations of 0.147 and 0.37 pc respectively. If we assume a typical outflow speed of about 15 km s$`^1`$then this gives us timescales of 9,590 yr and 18,100 yr for the formation of the inner to intermediate shell and the intermediate to outer shell respectively. Therefore, we could conclude that the shells have been formed by mass-loss modulations separated by these timescales. We therefore derive timescales which are more than 2 times longer for the time interval between the inner and intermediate shells and more than 10 times shorter for the time interval between the intermediate and outer shells. If we include the gas velocity information the situation becomes even worse; we derive timescales of 65 and 14,500 yr respectively.

We find that velocity gradients play a key role in the evolution of the shells and that we cannot infer the mass-loss modulations that give rise to their formation from observations. Observations do not provide much information about the mass-loss history.

### 4.4 Implications for the stellar evolution models

Our simulations provide an external shell slightly larger than what is observed on average. According to our results the size of the external shell is mainly determined by the AGB evolutionary time. If this is the case, shorter AGB evolutionary times than those obtained by VW93 could connect our models to the existing data. It appears that the observed AGB evolutionary time (obtained from HR diagrams) is shorter than that computed by stellar evolution models. The AGB evolutionary time could be reduced by increasing the mass-loss rate during the thermal pulses. Moreover, the maximum mass-loss rate in VW93 models was limited to the radiation pressure value which corresponds to the maximum momentum transfer due to simple ‘scattering’ from the electrons to the gas. It is known that multiple scattering can significantly increase this amount (by factors from 5 to 10) and it could be that the maximum mass-loss rates are underestimated. If the increase in mass-loss rate considering multiple scattering does not shorten the AGB evolution sufficiently, other mechanisms together with radiation pressure on dust grains should be invoked in the mass-loss rate processes (i.e magnetic fields as proposed by Pascoli 1994). It has also been argued that in order to account for the formation of molecular species such as titanium carbides observed in CSE around AGB stars (Molster 2000) higher mass-loss rate than those predicted by current models should be invoked. The problem is that observations of the CSE cannot constrain the numerical simulations. The difficulties in detecting these huge circumstellar shells at infrared wavelengths bias the observations and hence the constrains we can place on the stellar evolution calculations using our models.

## 5 Conclusions

We have studied the response of the circumstellar gas to the large temporal variations of the wind during the TP-AGB given by VW93 taking into account the influence of the external ISM. We shown how the mass-loss and velocity variations associated with the thermal pulses lead to the formation of transient features in the density structure and that these features do not survive at the end of the TP-AGB evolution. The evolution of the circumstellar gas during the AGB phase can be said to be highly non-linear due to the hydrodynamical processes which the gas experiences.

The formation of stable shell structures occurs via two main processes: shocks, between the shells formed by two consecutive enhancements of the mass-loss, for stellar masses 1, 1.5, 2 and 2.5 M$`_{}`$, and for stellar masses 3.5 and 5 M$`_{}`$by continuous accumulation of the material ejected by the star in the interaction region with the swept up ISM gas.

According to our models a large fraction of the mass of the external shells (up to $``$ 70 % for the 1 M$`_{}`$stellar model) belongs to the ISM matter that had been swept up by the stellar wind. This amount decreases according to the assumed density of the ISM, and according to the progenitor mass; the higher the progenitor mass the lower the fraction of the ISM material residing in the shells.

We quantified the effect of the ISM pressure on the circumstellar gas structure, and found that the presence of the ISM provides a non-negligible source of pressure that decelerates and compresses the external shells. The velocity of the external shell increases by a factor of $`1.5`$ and the density decreases by a factor of $`0.2`$ when the interstellar medium provides ten times less thermal pressure.

From our simulations we find that the history of mass-loss history experienced by the star is not recorded in the density or velocity structure of the gas. Therefore, does not seem possible to use observations to recover information on the later stages of the stellar evolution, such as the number of thermal pulses experienced by the star or the time between consecutive mass-loss events.

According to our simulations, large regions of neutral gas surrounding the molecular envelopes of AGB stars should be expected. These large regions of gas are formed from the mass-loss experienced by the star during the AGB evolution. These large shells should be detectable at infrared wavelengths.

We thank M. L. Norman and the Laboratory for Computational Astrophysics for the use of ZEUS-3D. We also want to thank Emanuel Vassiliadis for his fruitful comments at the begining of this work and for providing us with some of his models. EV is grateful to Tariq Shahbaz and Letizia Stanghellini for their careful reading of the manuscript and their valuable comments. The work of EV and AM is supported by the Spanish DGES grant PB97-1435-C02-01. GGS is partially supported by grants from DGAPA-UNAM (IN130698, IN117799 & IN114199) and CONACyT (32214-E).


format: text

text:
# Direct Measurement of the Photon Statistics of a Triggered Single Photon Source

## Abstract

We studied intensity fluctuations of a single photon source relying on the pulsed excitation of the fluorescence of a single molecule at room temperature. We directly measured the Mandel parameter $`Q(T)`$ over 4 orders of magnitude of observation timescale $`T`$, by recording every photocount. On timescale of a few excitation periods, subpoissonian statistics is clearly observed and the probablility of two-photons events is 10 times smaller than Poissonian pulses. On longer times, blinking in the fluorescence, due to the molecular triplet state, produces an excess of noise.

Over the past few years, there has been a growing interest for generating a regular stream of single photons on demand. This was mainly motivated by applications in the field of quantum cryptography Gisin et al. (2001). An ideal single photon source (SPS) should produce light pulses containing exactly one photon per pulse, triggered with a repetition period $`\tau _{\mathrm{rep}}`$, and delivered at the place of interest with 100% efficiency. For any given measurement time $`T`$, this source would emit exactly $`N=T/\tau _{\mathrm{rep}}`$ photons, so that the standard deviation $`\mathrm{\Delta }N\sqrt{N^2_TN_T^2}=0`$ ($`_T`$ has to be understood as a mean value over a set of measurements lasting $`T`$). Such a source would then be virtually free of intensity fluctuations, therefore corresponding to perfect intensity squeezing Reynaud (1990).

A first category of SPSs already realized consists of sources operating at cryogenic temperature. They rely on optically Kim et al. (1999); Michler et al. (2000); Santori et al. (2001); Moreau et al. (2001) or electrically Yuan et al. (2002) pumped semiconductor nanostructures or on the fluorescence of a two level system coherently prepared in its excited state Brunel et al. (1999). A one–atom micromaser has also been used to prepare arbitrary photon number states on demand Brattke et al. (2001). However the collection efficiency of photons is barely higher than a few $`10^3`$ in these experiments. Due to this very strong attenuation, the intensity statistics are very close to a Poisson law at the place where the stream of photons is available. Another route is to realize SPSs at room temperature. In this case higher collection efficiency (around 5%) is achieved. The existing room–temperature SPSs rely on the pulse saturated emission of a single 4–levels emitter Martini et al. (1996); Brouri et al. (2000).

When the pulse duration $`\tau _\mathrm{p}`$ is much shorter than the dipole radiative lifetime $`\tau _{\mathrm{rad}}`$, such a single emitter can only emit one photon per pulse. This temporal control of the dipole excitation allows therefore to easily produce individual photons on demand Lounis and Moerner (2000); Beveratos et al. (2002). However, in previous SPS realizations, little attention has been paid to analyse their intensity fluctuations. To address this problem we realized a room–temperature SPS relying on the pulsed saturation of a single molecule embedded in a thin polymer film Treussart et al. (2001).

The samples are made of cyanine dye DiIC18(3) molecules dispersed at a concentration of about one molecule per 10 $`\mu m^2`$ into a 30 nm thick PMMA layer, spincoated over a microscope coverplate. The fluorescence from the single molecule is excited and collected by the standard technique of scanning confocal optical microscopy Nie and Zare (1997). The molecules are non–resonantly excited at 532 nm, with femtosecond pulses generated by a Ti:Sapphire laser and frequency doubled by single pass propagation into a LiIO3 crystal. The repetition rate, initially at 82 MHz, is divided by a pulsepicker. The energy per pulse $`E_\mathrm{p}`$ is adjustable by an electro–optic modulator. The pulse duration is about $`\tau _\mathrm{p}100`$ fs. The excitation light is reflected by the dichroic mirror of an inverted microscope, and then focused by an oil–immersion objective ($`\times 60`$, NA=1.4), to form a spot of $`1\mu `$m2 surface area. The fluorescence light from the sample is collected by the same objective and then focused into a 30 $`\mu `$m diameter pinhole. After recollimation, a holographic notch filter removes the residual pump light. A standard Hanbury Brown and Twiss (HBT) setup is then used to split the beam and detect single photons on two identical avalanche photodiodes. Glass filters are placed onto each arm to suppress parasitic crosstalk Kurtsiefer et al. (2001) between the two photodiodes.

In order to rapidly identify single molecule emission, we first measure the intensity autocorrelation function $`I(t)I(t+\tau )`$ of the fluorescence light by the standard Start–Stop technique with a time–to–amplitude converter Brunel et al. (1999). When a single emitter is addressed, there is virtually no event registered at $`\tau =0`$, since a single photon cannot be simultaneously detected on both sides of a beamsplitter Grangier et al. (1986). The histogram shows a peak pattern at the pulse repetition period $`\tau _{\mathrm{rep}}`$. As explained in Ref.Brunel et al. (1999), the peaks’ areas allow one to infer the probabilities $`P_\mathrm{S}(n)`$ for the source (S), of giving $`n=0,1,2`$ photocounts per excitation pulse, where 2 photons counts are due to deviation from the ideal SPS. Nevertheless, this technique can hardly be used to extract the intensity fluctuations on timescale longer than a single pulse. We have therefore chosen to record each photodetection event with a two–channel Time Interval Analyser computer board (GuideTech, Model GT653). Since each detection channel has a deadtime of 250 ns the excitation repetition rate was chosen to be 2 MHz.

In a typical experiment, we first raster scan the sample at low excitation energy per pulse (0.5 pJ). When a single molecule is located, we apply the excitation energy ramp $`E_\mathrm{p}(t)`$ shown on inset of Fig.1, and simulteaneously record the fluorescence counts on a 1 ms integration time. Fig.1 displays the fluorescence counting rate $`R`$ vs. $`E_\mathrm{p}`$. The large intensity fluctuations are due to triplet state excursion of the molecule (see Fig.2). If this state is not taken into account, the molecular energy levels can be modelised by a 2-level system assuming a very fast non-radiative relaxation between the two higher and the two lower energy states. The excited state population $`\sigma `$ at the time $`\tau _\mathrm{p}`$ after the pulse arrival is then

$$\sigma =\frac{E_\mathrm{p}/E_{\mathrm{sat}}}{(1+E_\mathrm{p}/E_{\mathrm{sat}})}\left(1e^{\frac{\tau _\mathrm{p}}{\tau _{\mathrm{rad}}}\left(1+\frac{E_\mathrm{p}}{E_{\mathrm{sat}}}\right)}\right),$$
(1)
where $`\tau _{\mathrm{rad}}2.8`$ns for the cyanine considered. The data $`R(E_\mathrm{p})`$ are fitted by the function $`R=R_0\times \sigma `$ in a two steps procedure. After a first fit of the raw data, all the points below this fit, which are attributed to triplet state excursion, are removed. The fit of the remaining set of data yields $`R_0=160\times 10^3`$ counts/s and $`E_{\mathrm{sat}}=5.6\times 10^5`$ pJ.

In order to optimize the number of emitted photons and avoid rapid photobleaching, we then set $`E_\mathrm{p}^{\mathrm{max}}`$ to 5.6 pJ. Such a value would correspond to $`\sigma =97\%`$, for the molecule studied in Fig.1. During the constant maximum pumping energy period of the excitation ramp, $`10^4`$ detection events are typically recorded before photobleaching. Thanks to the high stability of the period of the pulsed laser, this set of times can be synchronized on an excitation timebase. We then build the table of the number of photocounts $`n_i=0,1,2`$ for each excitation pulse $`i`$. Photons which are delayed by more than $`10\times \tau _{\mathrm{rad}}`$ are considered to come from the dark counts of the two photodiodes, and are therefore rejected.

The data considered hereafter corresponds to a molecular source (S) which survived during 319769 periods (about 160 ms) yielding 14928 recorded photons including 14896 single photon events, 16 two-photons events. We deduced $`P_\mathrm{S}(1)=0.0466`$ and $`P_\mathrm{S}(2)=5.0\times 10^5`$ and a mean number of detected photon per pulse $`\overline{n}_\mathrm{S}=0.0467`$ (see Table 1). The real source is considered as the superposition of an attenuated ideal SPS with an overall quantum efficiency $`\eta `$, and a coherent source simulating the background, which adds a mean number of detected photon per pulse $`\gamma `$. From the measured values of $`P_\mathrm{S}(1)`$ and $`P_\mathrm{S}(2)`$, we infer $`\eta 0.0445`$ and $`\gamma 2.2\times 10^3`$. This leads to a signal-to-background ratio of about 20.

We also compared experimentally our SPS to a reference source (R) made of attenuated pump laser pulses, with approximately the same mean number of detected photons per pulse. Quantitative tests of this reference source and of the detection setup are however necessary. A particular care has to be taken to the bias of photocount statistics due to the detection deadtimes on both channels. For a coherent state of light (C) containing $`\alpha `$ photon per excitation pulse, one can readily calculate the counting probability distribution and show that $`P_\mathrm{C}(0)=e^\alpha ,P_\mathrm{C}(1)=2e^{\alpha /2}(1e^{\alpha /2})`$, $`P_\mathrm{C}(2)=(1e^{\alpha /2})^2`$ and $`\overline{n}_\mathrm{C}=2(1e^{\alpha /2})`$, where $`\overline{n}_\mathrm{C}`$ is the mean number of detected photons per pulse. For the reference source (R), we measured $`\overline{n}_\mathrm{R}=0.0462`$, $`P_\mathrm{R}(1)=0.0452`$, $`P_\mathrm{R}(2)=50\times 10^5`$, whereas one predicts, for $`\overline{n}_\mathrm{C}=\overline{n}_\mathrm{R}`$, $`P_\mathrm{C}(1)=0.0451`$ and $`P_\mathrm{C}(2)=53\times 10^5`$. The measured values are in good agreement with the predictions, which proves that the faint Ti:Sapphire pulses make a good calibration source for Poisson statistics. We then infer the ratio $`P_\mathrm{S}(2)/P_\mathrm{R}(2)=0.10`$, which tells that the number of two photons pulses in our SPS, is 10 times smaller than in the reference poissonian source (R).

In a first attempt to estimate the fluorescence intensity fluctuations per pulse, we considered samples of the data made of $`W`$ excitation cycles. We introduced a normalized variance $`V_\mathrm{W}`$ defined, on the sample, by $`V_W(\mathrm{\Delta }n)^2_\mathrm{w}/n_\mathrm{w}`$, with $`(\mathrm{\Delta }n)^2_\mathrm{w}_{i=1}^W(n_in_\mathrm{w})^2/W`$, where $`n_i`$ is the number of detected photons for the pulse $`i`$ and $`n_\mathrm{w}`$ is the mean number of detected photon per pulse in the sample considered. In the very few samples for which $`n_\mathrm{w}=0`$, $`V_W`$ is not defined and is set to 1. For a Poisson distribution of photocounts $`V_W=1`$, whereas $`V_W<1`$ for subpoissonian distribution.

In order to follow the time evolution of intensity fluctuations, we then extract from the set of $`\{n_i\}`$ all the successive samples of photocounts measurements of size $`W=100`$, separated by a single pulse period. Fig.2(a) displays the mean number of detected photons per pulse vs. time. We clearly see random intermittency in the fluorescence of the molecule, due to the presence of a dark triplet state $`T_1`$ in the molecular energy levels diagram (see inset of Fig.2(a)). At each excitation cycle the molecule has a small probability to jump into this non–fluorescent state, where it stays for a time much longer than the repetition period. Fig.2(b) shows, in parallel, the timetrace $`V_W(t)`$ of the normalized variance. During an emission period, $`V_W`$ stays below 1, and the statistics of the number of the detected photons per pulse is subpoissonian. On the other hand, when the molecule stops to emit, the background light yields $`V_W1`$. If we now consider the whole set of data, our measurements yields a single value for the variance $`V=0.9555`$. In this intensity fluctuation analysis at the level of a single pulse, this value of $`V`$ is also directly related to the Mandel paramater Loudon (2000) $`Q(\mathrm{\Delta }n)^2/n1`$ by $`Q=V1=0.0445`$. Let us point out that due to the photodetection deadtime, the triggered reference source (R) also yields a subpoissonian counting statistics. More precisely, for the coherent source (C) giving about the same mean number $`\overline{n}_\mathrm{C}`$ of photons per pulse than our SPS, one predicts a value $`Q=\overline{n}_\mathrm{C}/2=0.0231`$. This is confirmed by our measurements on the reference source (see Table 1). Nevertheless, the fluctuations of the number of detected photons per pulse coming out of our SPS show a clear departure from the reference coherent source. Albeit still limited by the quantum efficiency $`\eta `$, this direct measurement of $`|Q|`$ is larger by more than one order of magnitude compared to previous experiments Short and Mandel (1983); Diedrich and Walther (1987). For such a solid state SPS like ours, any improvement achieved in the light collection efficiency would therefore yield higher values of this subpoissonian character. We indeed observed, in preliminary experiments, an increase of the collection efficiency by placing the molecule at a controlled distance of a metallic mirror.

However, the leak in the dark triplet state induces correlations between consecutive pulses. The measurement of the variance $`V`$ of the detected photon number per pulse is therefore insufficient to characterize the noise properties of our SPS. Whereas such a characterization is usually infered from the record of noise power spectra, our photocount measurements are performed in the time domain. We therefore introduce, as a new variable, the number $`N(T)=_{i=1}^kn_i`$ of the detected photons during any period of observation $`T=k\tau _{\mathrm{rep}}`$. The analysis of the fluctuations of the variable $`n_i`$ can be generalized to the variable $`N(T)`$, by using the time dependent Mandel parameter Mandel (1979) $`Q(T)(\mathrm{\Delta }N)^2_T/N_T1`$. We can also define a Mandel parameter $`Q_\mathrm{s}(T)`$ for the number of photons *emitted* by the source in the same period of time $`T`$. In the case of an ideal SPS, we have $`Q=\eta \times Q_\mathrm{s}`$ Abate et al. (1976). For such a source, $`Q_\mathrm{s}=1`$, and therefore $`Q(T)=\eta `$, for any value of $`T`$.

Fig.3 shows that we did observe subpoissonian intensity fluctuations on timescales from $`T=1\times \tau _{\mathrm{rep}}`$ to $`T=8\times \tau _{\mathrm{rep}}`$, with the minimum value $`Q(\tau _{\mathrm{rep}})=0.0445`$ achieved on a single pulse timescale, as explained above. When we consider the number of detected photons on timescales larger than $`10^5`$ s, the intensity fluctuations exhibit a superpoissonnian behaviour ($`Q(T)>0`$) as shown on inset of Fig.3. This is a direct consequence of the bunching due to the triplet state Bernard et al. (1993). We developped a simple model to account for the intermittency of the SPS emission. In this model, the molecule is either available for fluorescence and is said to be in a ON state, or it is in its triplet OFF state and does not fluoresce. Let us note $`p`$, the probability per unit of time to make a ON $``$ OFF transition, and $`q=1/\tau _\mathrm{T}`$ the one to make the reverse OFF $``$ ON transition, where $`\tau _\mathrm{T}`$ is the lifetime of the triplet state. Note that $`p\tau _{\mathrm{rep}}=𝒫_{\mathrm{ISC}}`$ is the intersystem crossing probability per excitation pulse. From measured values at the single molecule level with DiIC18(3) cyanine dye Veerman et al. (1999), $`p\tau _{\mathrm{rep}}10^41`$ and $`q\tau _{\mathrm{rep}}2.5\times 10^31`$. In this limiting case, the Mandel parameter of the source is

$$Q_\mathrm{s}(k\tau _{\mathrm{rep}})=\frac{2\times 𝒫_{\mathrm{ISC}}}{\beta ^2}\left\{1\frac{1}{k\beta }\left[1(1\beta )^k\right]\right\}1$$
(2)
where $`\beta (p+q)\tau _{\mathrm{rep}}=𝒫_{\mathrm{ISC}}+\tau _{\mathrm{rep}}/\tau _\mathrm{T}`$. The Mandel parameter of the detected photon counts is then $`Q(T)=\eta \times Q_\mathrm{s}(T)`$. As shown on Fig.3, our data for $`Q(T)`$ are well fitted by eq.(2) over more than 4 orders of magnitude, with $`\eta =0.0445`$ (measured) and the free parameters $`p`$ and $`q`$. The fit yields $`p\tau _{\mathrm{rep}}2\times 10^4`$ and $`\tau _\mathrm{T}250\mu `$s, in good agreement with Ref.Veerman et al. (1999).

As a conclusion, the record of every photocount time allows one to make a direct time domain fluctuation analysis, as presented in this Letter. This technique can be straightforwardly applied to other SPS, and is also suited to investigate photochemical propreties at the single molecule level.

###### Acknowledgements.

We are grateful to Carl Grossman and Philippe Grangier for help and stimulating discussions. This work is supported by France Telecom R&D and ACI “jeunes chercheurs” (Ministère de La Recherche et de l’Enseignement Supérieur).


text:
# 1 Centroid-shift plot from 𝛾–𝛾–t coincidences in 98Sr (filled circles). The scattering of the points gives an indication of the systematical errors. The prompt curve (dashed line) is obtained by an interpolation using transitions in the A=97 and 98 chains (crosses).

A tentative 4- isomeric state in 98Sr
B. Pfeiffer1, G. Lhersonneau1,a, H. Gabelmann1,b, K.–L. Kratz1, and the ISOLDE-Collaboration2

1 Institut für Kernchemie, Universität Mainz

2 CERN, Geneva, Switzerland

In the frame of our systematic investigations of nuclear structure properties of neutron-rich A$``$100 nuclei at the isotope separator ISOLDE at CERN, an isomeric level at 1619 keV with the surprisingly long $`t_{1/2}`$ = 85 ns was observed in 100Sr . In analogy to a $`K^\pi `$=4 rotational band observed in the isotone 102Zr , the half-life could be interpreted as a $`\mathrm{\Delta }`$K = 4 hindered decay from the two quasi-neutron configuration ($`\nu `$3/2$``$$`\nu `$5/2)4-. These configurations are supported by quantum Monte-Carlo deformed shell calculations, further indicating a strong weakening of the pairing residual interaction .

This initiated a reexamination of $`\gamma `$-$`\gamma `$-t coincidence data on the $`\beta `$-decay of 98Rb collected at the former OSTIS separator at the ILL (Grenoble) and during the experiment on 100Rb decay. As a first result, the halflife of the 0$`{}_{2}{}^{}{}_{}{}^{+}`$ level of 98Sr at 215 kev could be remeasured with higher precision . In addition, an evaluation of the centroid shifts of the intense transitions revealed the existence of lifetimes related to the transitions at 140, 145 and 1693 keV, respectively (see Fig. 1 of Ref. ). The delay of $`t_{1/2}`$ = 7 ns for the 140-1693 keV time distribution was attributed to the 1838 keV level (see Fig. 2 of Ref. ). However, $`I^\pi `$=3+ was assigned to this level after a new analysis of the data, leading to a revision of the earlier level scheme of 98Sr . In accordance with quantum Monte-Carlo pairing calculations, this isomeric state could now be explained as the deformed two quasi-neutron configuration ($`\nu `$9/2$``$$`\nu `$3/2)3+ .

In order to search for an analogue to the 4- state at 1619 keV in 100Sr and 1821 keV in 102Zr, respectively, centroid shifts of higher energy transitions are displayed in Fig. 1. This figure is an extension to Fig. 3 of Ref.  and detailed information can be found therein. Apart from the 1693 kev line, only two transitions show a detectable delay: 2498 keV ($`t_{1/2}`$ = 3$`\pm `$1 ns) and 2526 keV ($`t_{1/2}`$ = 5$`\pm `$1 ns), respectively. Whereas the line at 2498 keV was placed in the decay scheme of 98Sr, the line at 2526 keV could (up till now) not be attributed to any of the nuclei in the A=97 and 98 decay chains following the $`\beta `$\- and $`\beta `$-delayed neutron decay of 98Rb. K-hindered decays should exist in the strongly deformed A=98 Sr and Y isotopes, but feeding from $`\beta `$-decay of the 0+ ground-state of Sr ought not feed high-spin states in Y. A placement of the 2526 keV line as a ground-state transition in 98Sr, possibly fed by a hinderd low-energy transition, seems a plausible explanation.

The 2498 keV line is in coincidence with the 129 and 289 keV transitions in 98Sr, exclusively deexciting the 2932 keV level to the 4+ member of the ground-state band. This single decay branch and the hindrance of the 2498 keV transition suggest an analogy with the 1619 keV level in 100Sr . The higher excitation energy is in the range expected for such a level due to the lowering of the Fermi level for neutrons with respect to 100Sr. Consequently, we propose the ($`\nu `$3/2$``$$`\nu `$5/2)4- configuration for the 2932 keV level.

Further studies with new technologies are needed to firmly (dis)prove these assumptions. New measurements ought to be performed using the higher beam intensities and better detection systems now available at CERN-ISOLDE. With the installation of big Ge-arrays, the studies of prompt deexcitation of fission products have gained new interest , offering a complementary approach to decay studies. In the particular case discussed here, they could be used to search for the band expected to be built on the isomer.

Present address:
a INFN Laboratori Nazionali di Legnaro, Legnaro, Italy
b KSM-Analytik, Mainz