text:
# Modelling of Intermediate Age Stellar Populations: III- Effects of Dust-Shells around AGB Stars

## 1 Introduction

Stars in the AGB phase are a dominant source of the near-infrared light of intermediate-age (0.1 Gyr$``$age$`\mathrm{\hspace{0.17em}2}`$Gyr) stellar populations. As an example, AGB stars (defined as stars with M$`{}_{bol}{}^{}3.6`$) contribute up to 60% of K-band light at 1 Gyr (Ferraro et al. 1995, Mouhcine & Lançon 2002a).

AGB stars develop at the end of their life a strong mass loss. The optical light of AGB stars with the highest mass loss rate is almost entirely absorbed by their dusty circumstellar envelopes and re-emitted at longer wavelengths. These obscured objects become very bright infrared sources with particularly red near-infrared colours. Omont et al. (1999) have shown that at least 25% of the RGB stars in the galactic bulge are losing mass, based on their very red ISO colours.

Another characteristic of AGB stars is the formation of carbon stars (i.e., AGB stars with C/O$`>\mathrm{\hspace{0.17em}1}`$). It has been well known from counts in the Magellanic Cloud clusters (e.g. Aaronson & Mould 1980, Frogel et al. 1990, Rebeirot et al. 1993) that carbon stars represent a significant fraction of the AGB population of intermediate age clusters, at least at sub-solar metallicity. Hence, they are responsible of a significant fraction of the integrated light of intermediate age stellar populations (Costa & Frogel 1996). In addition, carbon stars tend to be the intrinsically brightest of the AGB stars present, as they usually correspond to later evolutionary stages than their M-type counterparts. It is well established observationally, at least for simple stellar populations, that carbon stars are redder than M stars (Persson et al. 1983).

Newly large and homogeneous statistical sample of photometric data sets in the near-infrared (2MASS, Skrutskie 1998, and DENIS, Cioni et al 1997) and in the mid-infrared (ISOGAL, Omont et al 1999), and in the optical (Zaritsky et al. 1997) give us a unique opportunity to study the late-type stellar content of our galaxy and the Magellanic Clouds, and to draw conclusions about their formation history.

The observed structure of the Hertzprung-Russell (HR) diagram of resolved galaxies is a convolution of two functions. One function describes the rate of the formation of the stars which are present today in a galaxy. The second function describes the evolution of a star in the observed HR diagram (i.e., lifetime of different stellar phases, the location of those phases in the HR diagram). At identical star formation histories, the apparent morphology is controlled mostly by the second function. Accurate stellar evolutionary models that cover the necessary range of ages, metallicities and evolutionary phases are needed to retrieve information about the star formation history of a galaxy. Unfortunately, the actual state of our understanding of late-type stellar evolution is still far from being reliable enough to perform such a goal accurately. The work presented in this series of papers aims at improving our knowledge of the effect of late type stars on stellar population properties.

Effects of the dusty envelopes surrounding AGB stars, in addition to those related to the formation of carbon stars, on (i) the morphology of the AGB sequence in the observational HR diagram of resolved stellar populations and (ii) integrated properties of unresolved populations have been neglected for a long time. Mouhcine & Lançon (2002a) have presented a new spectral library of single stellar populations where the formation of carbon stars was taken into account. They showed that the presence of carbon stars in the stellar populations lead to redder integrated near-infrared colours. Bressan et al (1998) have presented single burst population models taking into account the effect of the formation of dusty envelopes associated with AGB mass loss. They have used a radiative transfer model to correct the stellar spectra predicted by standard evolutionary tracks, assuming that all circumstellar envelopes have the same chemistry (i.e., silicate or carbonaceous grains). Hence when they consider carbon rich dust, the whole dust-enshrouded stellar population have the same properties, even low core mass stars that never form carbon stars. In addition, they have not taken into account the spectral dichotomy between oxygen rich and carbon stars, considering giant star spectra to be representative of the spectral energy distribution of AGB stars.

Based on the models of Mouhcine & Lançon (2002a), we constructed a set of theoretical models which account simultaneously for the formation of carbon stars and the dust shells around thermally pulsating-AGB (TP-AGB hereafter) stars, and obtained a new set of isochrones suited for the analysis of the near-infrared data. Our main goal, as noted above, is to study the evolution of the AGB sequence morphology and integrated near-infrared properties.

In Sect. 2 we present the grid of stellar evolutionary tracks which we use in order to calculate isochrones in the theoretical plane (i.e. luminosity vs. effective temperature). We recall briefly the modelling of the TP-AGB evolutionary phase, which is included by means of a synthetic model including all relevant physical processes that control the evolutionary properties that are related to their effects on integrated properties.

We describe the dusty envelope model used to calculate the effects of circumstellar shells around late-type stars on the energy distribution, and calculate the isochrones in the observational plane (i.e., magnitude vs. colour).

In Sect. 3 we study the effects of various TP-AGB subtypes on the evolution of near-infrared properties of both resolved and unresolved stellar populations. We discuss how they depend on age and metallicity. Finally, in Sect. 4 our conclusions are drawn.

## 2 Population synthesis modelling

### 2.1 AGB stellar evolution

Synthetic evolution models represent an attractive tool to investigate the evolutionary and chemically aspects of AGB stars. The evolution of TP-AGB properties are followed by means of such models. We present briefly the stellar evolution models that we used to follow the formation and the evolution of TP-AGB stars. The reader is referred to Mouhcine & Lançon (2002a) for more details.

The model includes different physical processes affecting the evolution of AGB stars along the TP-AGB phase and playing a dominant role in determining the lifetime, the extent of nuclear processing and the chemical abundances, like the envelope burning, the third dredge-up, the high rate mass loss at the tip of the TP-AGB phase. We start from the last model calculated along the early-AGB evolutionary tracks. The total mass, core mass, effective temperature, bolometric luminosity of the models, and carbon to oxygen ratio in the stellar envelope are let to evolve according to semi-analytical prescriptions.

Synthetic modelling of the TP-AGB phase requires the knowledge of: (i) the critical core mass M$`{}_{}{}^{min}{}_{c}{}^{}`$ above which dredge-up is triggered, (ii) the dredge-up efficiency (defined as $`\lambda =\mathrm{\Delta }M_{dredge}/\mathrm{\Delta }M_c`$, where $`\mathrm{\Delta }M_{dredge}`$ and $`\mathrm{\Delta }M_c`$ are, respectively, the amount of mass dredged-up to the envelope after a pulse and the increase in core mass during the time between two successive TPs), (iii) the critical envelope mass above which the envelope burning is effective, and (iv) the element abundances in the inter-shell material.

The basic physical prescriptions which the models stand on are the following:

1. the core mass/luminosity relation;
2. the core mass/interpulse-period relation;
3. the core mass evolutionary rate;
4. an algorithm to calculate the effective temperature;
5. CNO burning possibly occurring in the convective envelope base;
6. prescription to calculate the mass loss rate by stellar wind as function of the stellar parameters;
7. the third dredge-up and its efficiency;
8. the composition of the third dredge-up material.

The prescriptions used account for different physical processes affecting the core mass/luminosity and the core mass/interpulse-period relations, and do not rely on the assumption that the classical linear relations are valid. Marigo et al. (1999) have shown that the actual state of semi-analytical modelling account for both breaks of the core mass-luminosity relation, for low-mass stars when diffusive overshooting is introduced (Herwig et al. 1998), and for massive stars due to the envelope burning (Blöcker 1995). The adopted relations are rather good representations of the outputs of the detailed numerical models (Wagenhuber & Groenewegen 1998).

The theoretical values of the dredge-up parameters are uncertain. Standard numerical models predict low dredge-up efficiencies (Boothroyd & Sackmann 1988), while models assuming diffusive overshooting predict highly efficient dredge-up events (Herwig et al. 1998). Actually both parameters $`\lambda `$ and M$`{}_{}{}^{min}{}_{c}{}^{}`$ affect the lifetime of carbon-rich phase and the age interval for which carbon stars are present in a stellar population. Both parameters are also expected to depend on properties such as the stellar metallicity or its instantaneous mass and structure, but the available information on these relations is neither complete nor easy to extrapolate (e.g. Herwig 2000, Marigo et al. 1999). Most of stellar evolution calculations fail to predict carbon stars at luminosities as low as observed, and are only able to obtain dredge-up for relatively high initial stellar masses (Forestini & Charbonnel 1997). The models presented here assume constant $`\lambda `$ and $`M_c^{min}`$. With this assumption, reproducing the luminosity function of carbon stars in the LMC requires intermediate dredge-up efficiencies and leads us to adopt $`\lambda =0.75`$ and $`M_c=0.58`$ (Groenewegen & de Jong 1993, Marigo et al. 1996).
The physical conditions needed to trigger the envelope burning, and its efficiency, are still uncertain and matter of debate. We assume that the envelope burning is effective only if the envelope mass exceeds a critical value (see Marigo, 1998, for alternative approach). To estimate this critical value, we use the analytical relation presented by Wagenhuber & Groenewegen (1998).
The chemical composition of the dredged-up material is taken from Boothroyd & Sackmann 1988.

The final calculations are performed using a mixing length parameter $`\alpha =\mathrm{\hspace{0.17em}2}`$, Blöcker’s (1995) mass loss prescription with a mass loss efficiency of $`\eta =\mathrm{\hspace{0.17em}0.1}`$ (Groenewegen & de Jong 1994). The end of TP-AGB phase is determined by the ejection of the stellar envelope under the influence of the stellar wind. The models provide the lifetime of the whole TP-AGB phase and its subtypes phases, in addition to the effective temperature and luminosity, as functions of time and initial mass. All this information is used to determine the spectral type of each star in the HR diagram. A crucial key that determines the nature of intermediate age populations, is the sensitivity of the TP-AGB star properties to the initial metal content. This will affect in particular the formation of carbon stars. We would like to note that the models discussed here have shown their ability to reproduce relevant constraints relative to intermediate age stellar populations as observed in the solar neighbourhood, in Magellanic Cloud globular clusters, or in Local Group galaxies (Mouhcine & Lançon 2002 a,b). Fig. 1 shows the predicted AGB sequence in the HR diagram for the stellar populations dominated by AGB stars at solar metallicity.

Let us discuss first what these models tell us about the effects of carbon stars and dust-enshrouded stars on the properties of single stellar populations. To account for both AGB sub-populations we need to know, for each stellar initial mass, the ages marking the transition from (i) oxygen-rich chemistry to carbon-rich chemistry, (ii) optically bright source to optically obscured source, and (iii) TP-AGB phase to Post-AGB phase. To estimate the age at which an AGB stars becomes obscured, we assume that the stars are not longer observable when the optical depth at $`1\mu `$m ($`\tau _{1\mu m}\dot{M}\mathrm{\Psi }`$, where $`\tau _\lambda `$ is the optical depth, $`\dot{M}`$ the mass-loss rate, and $`\mathrm{\Psi }`$ the dust-to-gas ratio) is higher than a critical value ($`\tau _{1\mu m}\tau _{crit}`$). See Sec.2.2 for more informations about the envelope models used to estimate the optical depth.

Fig. 2 shows the evolution of the fraction of the whole TP-AGB lifetime stars spend as obscured objects, as a function of their initial mass for Z=0.02 and Z=0.008 metallicities, for two different values of the critical optical depth at $`1\mu `$m ($`\tau _{1\mu m}=1and3`$). This choice is justified by the fact that the ingredients (grain properties, mass-loss rate determination, condensation temperature, … etc) used to construct the envelope models are known within a factor of few (3-5). The plot shows clearly that massive stars (i.e., M$`{}_{init}{}^{}\mathrm{\hspace{0.17em}4}3.5`$M) that suffer from envelope burning, spend a large fraction of their TP-AGB lifetime as obscured objects, due to their high mass-loss rate (see below). On the other hand, the plot shows that this fraction decreases rapidly for stars having initial masses below this limit. This plot shows again the effect of the envelope burning on intermediate age stellar population properties, and emphasizes the importance of its inclusion in stellar population models.

Figure 3 shows the evolution of the bolometric magnitude as a function of the initial mass for the three characteristic stages of the AGB phase. The lower and the upper solid lines represent the magnitudes when a star enter and quit the TP-AGB phase respectively, the dashed line represents the transition magnitude of a TP-AGB star from oxygen-rich to carbon-rich chemistry, while the dashed-dotted line represents the transition magnitude of a stars from being optically bright to being dust-enshrouded. The bolometric luminosity distribution as a function of mass shows that when both oxygen-rich and carbon stars coexist in a single age stellar population, carbon stars are the brightest, as they are more evolved objects than oxygen-rich stars. However, for single stellar populations with a turn-off mass of M$`{}_{TO}{}^{}\mathrm{\hspace{0.17em}4}M_{}`$ (i.e., stellar population age about $`\mathrm{\hspace{0.17em}0.2}`$Gyr), for which the competition between the envelope burning and the third dredge-up starts to vanish and the latter process becomes more effective than in more massive stars, oxygen-rich stars could be brighter. Note that at that age, most of the TP-AGB stars are losing mass at a relatively large rate.

Note that massive/luminous carbon stars (i.e., M$`{}_{init}{}^{}\mathrm{\hspace{0.17em}4}3.5`$M) spend their entire lifetime as optically invisible object. The TP-AGB evolution of those stars is controlled mainly by the interplay between the envelope burning, the mass-loss, and the third derdge-up. The overluminosity produced by the envelope burning, trigger the super-wind ($`\dot{M}\mathrm{\hspace{0.17em}10}^610^4`$) very early, and the hot temperature at the base of the envelope prevent the formation of carbon stars. This scheme is valid as long as the envelope mass is high enough to maintain a hot temperature at the base of the envelope. When the mass-loss has dramatically reduce the mass of the envelope, lowering the temperature at the base of the envelope, the envelope burning stops and the star may become carbon-rich, but being already embedded in dense circumstellar shell.

The bolometric luminosity distribution of dust-enshrouded M-type stars is expected to peak at brighter luminosities than the one of dust-enshrouded carbon stars, which, in turn, is expected to peak at higher luminosities than the luminosity distributions of both optically bright carbon stars and M-type stars respectively. Optically bright carbon stars are expected to populate the AGB mostly between M$`{}_{bol}{}^{}4.3`$ and $`5.5`$, with a peak around M$`{}_{bol}{}^{}4.8`$, whilst dust-enshrouded carbon stars may be as bright as M$`{}_{bol}{}^{}6.5`$. Nevertheless, independently of the initial metallicity dust-enshrouded carbon stars as faint as M$`{}_{bol}{}^{}4.5`$ may be found (with progenitor masses about M$`{}_{init}{}^{}\mathrm{\hspace{0.17em}1.5}M_{}`$). Note that the peak of the bolometric luminosity distribution of optically bright carbon stars is not expected to vary with the metallicity, because the preferred mass range for carbon star formation is relatively insensitive to metallicity (for more details, see Mouhcine & Lançon 2002b). The luminosity distribution of dust-enshrouded M-type stars is expected to have a peak at M$`{}_{bol}{}^{}6.5`$ to $`6.8`$. These stars are predicted to live $`\mathrm{\hspace{0.17em}4}9\times \mathrm{\hspace{0.17em}10}^4`$ yr in agreement with van Loon et al. (1999).

Note that Fig. 3 is constructed neglecting the luminosity modulation due to the recurrence of third dredge-up events. This point will be discussed in more detail in Sec. 3.1.

Olivier et al. (2001) have used a sample of oxygen rich and carbon rich dust-enshrouded AGB stars in the solar neighbourhood to estimate empirically the duration of the high mass-loss rate AGB phase. They show that the stars in their sample have a main-sequence progenitor with $`1`$M$`{}_{}{}^{}`$M$`\mathrm{\hspace{0.17em}2}`$M (with an average progenitor mass of $`\mathrm{\hspace{0.17em}1.3}`$ M) have an average dust-enshrouded AGB phase of ($`3.7\pm \mathrm{\hspace{0.17em}1.9}`$)$`\times \mathrm{\hspace{0.17em}10}^4`$ years. For the same parameters (i.e. M$`{}_{init}{}^{}\mathrm{\hspace{0.17em}1.3}`$M and Z=Z), our models predict the high mass-loss AGB phase have a duration of $`5.5\times \mathrm{\hspace{0.17em}10}^4`$ years, in agreement with the lifetimes derived by Olivier et al (2001) (see Habing 1996 for a discussion of other determination of the duration of the mass-loss phase).

### 2.2 Dust envelope models

Extensive work was done to study the properties of circumstellar dust shells around late-type stars taking into account radiative transfer and chemistry (Hummer & Rybicki 1971, Bedijn et al. 1978, Rowan-Robinson & Harris 1983, Martin & Rogers 1984, Justtanont & Tielens 1992). However, as our principal goal is to study the dust-enshrouded stellar populations as a whole, we need a dust envelope model, which is consistent with the observational constraints, and makes us able to figure out (i) the extinction that affects the stellar spectrum due to the dust shell, (ii) how much light is emitted in the near-infrared as function of the fundamental stellar parameters, and (iii) accounts for the average properties of the dust-enshrouded stellar populations.

To construct dusty envelope models around late type AGB stars, the following assumptions are made: (i) spherical symmetry, and (ii) constant properties of the dust grains at different ages and metallicities. We stress that our purpose is to have a good description of the average properties of grains in circumstellar envelopes. For a given mass loss rate and from the equation of continuity, the dust density outside the condensation radius, where the dust forms, scales as $`r^2`$:

$$\rho _d(r)=\frac{\dot{M_d}}{4\pi r^2v_d}$$
(1)
where $`\dot{M_d}`$ is the dust mass-loss rate, and vd is the dust outflow velocity. The grain emission being negligible in the near-infrared, the wavelength range of interest in this paper, the spectral energy distribution of mass-losing star is primarily determined by the dust optical depth. Note that for some extreme stars, grain emission may contributes to a relatively large fraction, 20%-50%, of K-band light. However these objects have very short lifetime and consequently have small effects on stellar population properties: a very small number of such stars may present in a stellar population, and contribute to a small fraction of integrated light budget.
The dust optical depth is defined as:

$`\tau _\lambda `$ $`=`$ $`{\displaystyle _{R_c}^{R_{out}}}\pi a^2Q_{ext}(\lambda )n_d(r)𝑑r`$ (2)
$`{\displaystyle \frac{3}{16\pi }}{\displaystyle \frac{\dot{M}\mathrm{\Psi }}{R_cv_d}}{\displaystyle \frac{Q_{ext}(\lambda ,a)/a}{\rho _{gr}}}`$
where Q$`{}_{ext}{}^{}(\lambda `$)=Q$`{}_{abs}{}^{}(\lambda )`$+Q$`{}_{sca}{}^{}(\lambda )`$ is the dust extinction, absorption and scattering coefficients respectively, nd is the dust grain number density, $`a`$ is the grain size, $`\dot{M}`$ is the total mass-loss rate, $`\mathrm{\Psi }`$ the dust-to gas ratio, $`\rho _{gr}`$ is the grain density, Rout is the outer radius of the circumstellar shell, and Rc the dust condensation radius. We will assume that the dust velocity equals the terminal velocity of the gas. Groenewegen (1993) has shown that this approximation is correct to within 10%. To derive Eq. (2) we have neglected 1/Rout with respect to 1/Rc.

For the near–infrared region we have assumed that the absorption and scattering efficiencies of the grains follow a power-law, Q($`\lambda `$) = Q($`1\mu `$m)$`(\lambda /\mu `$m$`)^\beta `$ (Jura 1983, Bedijn 1987). The grain size of both carbon rich dust and oxygen-rich dust was taken to be a = 0.1 $`\mu `$m, and the grain density was taken to be $`\rho _{gr}=\mathrm{\hspace{0.17em}2.5}`$gr cm-3 for silicates and $`\rho _{gr}=\mathrm{\hspace{0.17em}2.26}`$ gr cm-3 for amorphous carbon grains. The absorption efficiencies at the reference wavelength of 1 $`\mu `$m, $`Q(1\mu `$m), were taken from Jones & Merrill (1976) for silicates and from Rowan-Robinson (1986) for amorphous carbon grains. The value of $`\beta `$ which fits the energy distribution of carbon stars in the region of our interest is $`\beta =\mathrm{\hspace{0.17em}1.0}`$ (e.g. Rouleau & Martin 1991, see also Rengarajan et al. 1985, Le Bertre et al. 1995, Suh 2000). The same value was taken for silicate grains (Schutte & Tielens 1989). We consider silicate grains when dealing with oxygen rich circumstellar envelopes, and amorphous carbon grains when dealing with carbon rich circumstellar envelopes. We assumed that the properties of all the grains in the dusty envelope around a TP-AGB star change instantaneously when the star becomes carbon rich. The condensation radius was taken to vary like:

$$\frac{R_c}{R_s}\frac{1}{2}\left[\frac{T_{eff}}{T_c}\right]^{(4+\beta )/2}$$
(3)
where Tc is the condensation temperature at the inner radius and Rs is the stellar radius (i.e., $`R_sL^{1/2}T_{eff}^2`$). Groenewegen & de Jong (1994) have shown, using detailed radiative transfer calculations, that the evolution of the condensation radius as function of the stellar parameters, at least for carbon stars, is well represented by Eq. (3). The condensation temperature of both carbon-rich and oxygen-rich dust is taken to be $`\mathrm{\hspace{0.17em}1000}`$K (Le Bertre 1986, Le Sidaner & Le Bertre 1996).

To achieve our envelope models, we need to evaluate the dust-to-gas ratio $`\mathrm{\Psi }`$ and the expansion velocity $`v_{exp}`$ for stars evolving along the TP-AGB phase. To evaluate $`v_{exp}`$, the period-mass-radius relation of fundamental mode pulsators from Wood (1990) is used in combination with the expansion velocity-period relation from Vassiliadis & Wood (1993). Regarding $`\mathrm{\Psi }`$ scaling arguments show that $`\mathrm{\Psi }v_{exp}^2L^{1/2}`$. Using detailed modelling Habing et al. (1994) have shown that the exact dependence of $`\mathrm{\Psi }`$ on luminosity and vexp deviates slightly from the scaling formula. We will use the formula of Habing et al. Note that we do not introduce any additional explicit dependence on metallicity in our semi-analytical envelope model.

## 3 Effects of obscured AGB stars.

In the following subsections, we will discuss how dust-enshrouded AGB stars affect near-infrared properties of both resolved galaxies and unresolved stellar populations. For this purpose, we coupled the dusty envelope model with synthetic AGB models. Regarding resolved populations, we study the sensitivity of the AGB sequence morphology in the near-infrared HR diagram to the dominant AGB subtype. For unresolved stellar populations we will examine how optical/near-infrared broad-band colours are sensitive to the presence of obscured AGB stars.

Transformations from the theoretical to the observational plane were performed by making use of the stellar spectral library of Lançon & Mouhcine (2002). The latter library is based on stars on the solar neighbourhood. To account for a dependence on the metal content even for M giants, we adopted the same library at different metallicities but assigned the spectral class, identified by the (I-K) colour, adopting the (I-K)-Teff relation of Bessell et al. (1989) which depends on the metallicity. (R-H) was taken to be an effective temperature indicator of carbon stars (Loidl et al. 2001). The sensitivity of carbon star spectra to metallicity is thought to be small (Gautschy 2001). Hence carbon star spectra of solar neighbourhood stars were used at all metallicities.

To account for the effect of the dust along the TP-AGB, we combine the AGB star synthetic evolution models and the semi-analytical envelope models presented in Sec. 2.2. The optical depth is a combination of two effects. The first one (i.e., Q($`\lambda `$)/a $`\rho _{gr}`$) describes the circumstellar shell chemistry and the second one describes the effect of the evolutionary status of a TP-AGB star on its circumstellar shell structure. At fixed age, we use stellar evolutionary parameters of the TP-AGB stellar population that populate the isochrone (i.e., Teff, L, Minit, Z/Z) to derive the optical depth for each star, and hence to compute the near-infrared spectral energy distribution of mass-losing stars (Volk & Kwok 1988, Bressan et al 1998).

### 3.1 Resolved stellar populations: Isochrones in the Infrared

First, we examine the effect of carbon stars on the morphology of the upper cool/luminous part of the colour-magnitude diagram without considering obscuration. Fig 4 shows selected isochrones with Z=0.008 and Z=0.02 in the MK vs. (J-K) diagram.

Continuous lines show isochrones where oxygen-rich and carbon stars are accounted for. Dashed lines show isochrones for which only oxygen-rich stars are considered. The plot illustrates nicely the effect of carbon stars on the morphology of the upper colour-magnitude diagram of intermediate-age stellar populations. When carbon stars are neglected, the models predict that AGB star populations should appear as a single and continuous ”plume” in the MK vs. (J-K) diagram, with (J-K)$`\mathrm{\hspace{0.17em}0.8}1.5`$. However, when carbon star formation is followed, a carbon star branch is expected. This branch will appear in the colour-magnitude diagram as a second plume where (J-K) may be as red as 2 or redder for the same range of K-band magnitudes than the models that neglect the formation of carbon stars. The models predict that the first plume is occupied by oxygen-rich stars. Isochrones younger than $`\mathrm{\hspace{0.17em}0.2}0.3`$Gyr are not affected by the inclusion of carbon stars in the modelling for the simple reason that carbon stars are not present yet in the stellar population. The models predict that the carbon stars branch in the colour-magnitude diagram is more pronounced at Z=0.008 than at Z=0.02, because of the higher efficiency of carbon star formation in metal-poor systems.

Fig. 5 shows selected isochrones with Z=0.008 and Z=0.02 in the MK vs. (J-K) diagram. Now the formation of both carbon stars and circumstellar shells around the TP-AGB stars are accounted for to construct isochrones. A new morphological feature appears in the colour-magnitude diagram as an additional plume. The plot shows the location of three different populations in the MK vs. (J-K) diagram. The first stellar sub-population have $`0.6(JK)\mathrm{\hspace{0.17em}1.4}`$. This region of the HR diagram is occupied by Early-AGB stars and oxygen-rich TP-AGB stars. The magnitude extension of this sub-population is quite large because it is populated by stars with different progenitor initial mass (from 5-6 M to $`\mathrm{\hspace{0.17em}1}`$M), with quite different magnitude at the tip of the AGB. The second stellar sub-population have $`1.4(JK)\mathrm{\hspace{0.17em}2}`$, and it contains optically bright carbon-rich TP-AGB stars. The luminosity extent of this population is significantly narrower than that of the first population. The extent of the morphological feature related to this population is less pronounced at solar metallicity than at Z=0.008 metallicity, because of the larger efficiency of carbon star formation at low metallicity. The third stellar sub-population has extremely red colours (i.e., J-K$`\mathrm{\hspace{0.17em}2}`$), and decreasing K-band luminosity as the colour becomes redder. This region contains extremely red dust-enshrouded TP-AGB stars, both oxygen-rich and carbon-rich. Their large (J-K) and decreasing K-band luminosity are a consequence of large extinction due to dusty circumstellar envelopes.

The predicted morphology of this region of the HR diagram is in good agreement with the HR diagram morphology of the LMC stellar populations as observed by 2MASS (Nikolaev & Weinberg 2000).

To check further the behavior of our models in the near infrared region of the spectrum, we have compared our isochrones with samples of optically bright and dust-enshrouded stars. Fig. 6 (right panel) shows the selected isochrones at solar metallicity in the two colour diagram (J-H) vs. (H-K). Also shown is the location of M giants (area delineated by dashed line) and Miras (area delineated by continuous line) as observed in the Southern Polar Cap (Whitelock et al. 1994, 1995), in addition to carbon (filled circles) and oxygen-rich (open circles) dust-enshrouded AGB stars in the solar neighbourhood (Olivier et al. 2001). The left panel shows the same isochrones at Z=0.008 metallicity, compared to a sample of dust-enshrouded AGB stars in the LMC (Zijlstra et al. 1996), and the location of LMC/2MASS sources with (J-K)$`\mathrm{\hspace{0.17em}2}`$ (Nikolaev & Weinberg 2000).

Both the bulk of the M giants at (J-H)$`\mathrm{\hspace{0.17em}1}`$ and (H-K)$`\mathrm{\hspace{0.17em}0.3}`$ and of Miras samples are well reproduced by the models. As the star reaches the TP-AGB phase, the effect of the mass-loss increases and the circumstellar shells modify the energy distributions. Isochrones with different ages and metallicities reach different values of extinction. In the framework of our semi-analytical envelope model, the optical depth is an increasing function of mass and metallicity. The larger the metallicity and the lower the isochrone age, the larger the optical depth, and the further the isochrone extends in the colour-colour diagram. TP-AGB stars losing mass at high rates populate an extended red arm in the (J-H) vs. (H-K) two colour diagram. Fig. 6 shows that this red arm splits into two arms occupied by carbon stars and oxygen rich stars respectively. Dust-enshrouded carbon stars exhibit different colours from dust-enshrouded M-stars, indeed they have a redder (H-K) colour than M-type stars at fixed (J-H). This predicted behavior is in agreement with the observational constraints in the solar neighbourhood (Olivier et al. 2001). This behavior stresses again the importance of the inclusion of both AGB subtypes simultaneously to model intermediate-age stellar populations.

We have to mention that a fraction of carbon star subpopulation may occupy bluer and fainter regions in the HR diagram, than what the bulk of the subpopulation does. These relatively blue carbon stars may be low-mass carbon stars evolving through the faint/hot luminosity-dip phase of the quiescent hydrogen burning phase that follows immediately a thermal pulse. It is known that this luminosity-dip phase affects strongly the faint tail of carbon stars luminosity function (Marigo et al. 1996, see also below). If the star formation in a galaxy proceeds continuously, one may expect that a fraction of carbon stars should mix with oxygen-rich stars at bluer colours; this was shown by Kontizas et al. (2001) in their survey of carbon stars in the LMC where they found that a significant number of spectroscopically selected carbon stars are bluer and fainter than ”carbon star’s region” defined photometrically by Nikolaev & Weinberg (2000). The implication of the presence and the properties of this subpopulations on the the link between carbon star statistics and the star formation history, and hence the usage of carbon stars as quantitative star formation tracer are left to a forthcoming paper, where we present comparisons of our models to data (Mouhcine et al. 2002, in preparation).

The results and discussions presented until now were drawn out using smooth tracks; i.e. tracks where the stellar luminosity modulation by thermal pulses was averaged out and the resulting brightening rate of carbon stars is nearly constant. In fact TP-AGB stars undergo a long-lived underluminous stage, soon after the occurrence of a thermal pulse, with deeper depth and longer duration for low-mass stars (M$`{}_{init}{}^{}\mathrm{\hspace{0.17em}3}`$M). This effect plays a crucial role in determining the extension of the low-luminosity tail of the luminosity function of carbon stars (Boothroyd & Sackmann 1988, Groenewegen & de Jong 1993).

Fig. 7 shows bolometric and K-band luminosity functions of carbon stars that exist in simple stellar populations of selected ages, where roughly all carbon stars have the same initial mass of the progenitor noted in each panel. From this figure, it turns out that the faint tail of the bolometric luminosity function extends $`\mathrm{\hspace{0.17em}0.8}1`$ mag fainter than what is expected from the smoothed tracks. Stars with initial masses M$`{}_{init}{}^{}\mathrm{\hspace{0.17em}2.5}`$M have luminosity functions that peak almost at the same magnitude (see also Marigo et al. 1999). Note that few stars brighter than the tip of the AGB magnitude exist in the stellar population. Those stars evolve on the rapid luminosity peaks related to the thermal pulses. Additional large-amplitude pulsations of the AGB stars could affect the magnitude extent of different AGB sub-populations. Bolometric amplitudes of galactic carbon stars are $`\mathrm{\hspace{0.17em}0.6}`$ mag (Le Bertre 1992) and may be similar to those of galactic OH/IR stars (Le Bertre 1993).

The right panel of Fig. 7 shows comparisons between K-band luminosity functions of carbon stars of populations where all carbon stars are considered to be optically bright (continuous line), and stellar populations where both optically-bright and dust-enshrouded carbon stars coexist (dashed-dotted line). Two features arise. (i) When both sub-types of carbon stars are considered, the peak of the luminosity functions shifts to fainter magnitudes in comparison to what is expected when the formation of shells around carbon stars is ignored, and (ii) a faint and extended luminosity function tail (M$`{}_{K}{}^{}5`$) is predicted. The stars that populate this faint wing are mainly stars close to the end of their AGB phase having high mass-loss rate and extreme red colours (J-K$`\mathrm{\hspace{0.17em}3}`$), and most of their light is emitted at wavelengths longer than the K-band. As a consequence of the formation of this faint tail, the main body of the K-band luminosity function becomes narrower than what is predicted when the formation of dust-enshrouded carbon stars is neglected.

### 3.2 Unresolved stellar populations

In this section we will examine the effect of the dusty envelopes around TP-AGB stars on the integrated near-infrared properties of intermediate-age stellar populations. We will restrict these investigations to the case of instantaneous burst models (occurring at t=0).

Once the optical depth of a TP-AGB star is calculated, its ”reddened” spectrum is added to the stellar population spectrum weighted by the initial mass function (hereafter IMF). In the rest of the paper we adopt a Salpeter IMF between the lower and upper cut-off masses, M$`{}_{low}{}^{}=\mathrm{\hspace{0.17em}0.1}M_{}`$ and M$`{}_{up}{}^{}=\mathrm{\hspace{0.17em}120}M_{}`$ respectively ($`\varphi (m)m^\alpha ,\alpha =\mathrm{\hspace{0.17em}2.35}`$). Fig. 8 shows the temporal evolution of selected optical/near-IR colours when the formation of dust enshrouded AGB stars is taken into account at solar and Z=0.008 metallicities. For comparison, the evolution of optical/near-IR colours neglecting the effects of circumstellar shells is also plotted. The most striking feature is that models that account for the formation of dust-shells around the TP-AGB stars have bluer colours than models that do not account for them (see also Bressan et al. 1998). This behavior is due to the large obscuration affecting the brightest TP-AGB stars leading to a decrease of their near-IR flux. Keep in mind that the optical light is dominated by other stellar populations than AGB stars. Using these models for age-dating of stellar populations will lead to assigning older ages than what is predicted by models assuming pure photospheric emission. As the age of stellar population proceeds, the departure from the predictions of pure photospheric emission models decreases translating the fact that red giant branch stars become the dominant source of the near-IR spectral energy.

Metallicity effects are depicted in Fig. 9 where we plot the temporal evolution of differences between the synthetic colours of models that differentiate optically-bright and dust-enshrouded stars and models that do not for Z=0.02 and Z=0.008. For young stellar populations (i.e. age$`\mathrm{\hspace{0.17em}0.25}`$Gyr), when carbon stars do not yet exist, higher metallicity models are more affected by the inclusion of dust shells around TP-AGB stars because of their low effective temperatures which lead to higher mass loss rates and hence higher optical depth. The bump-like shape that appears around $`\mathrm{\hspace{0.17em}0.25}`$Gyr and extends until $`\mathrm{\hspace{0.17em}1}`$Gyr is related to the presence of carbon stars in the stellar populations. During this age range the effect of dust shells is larger at low metallicity even though metal-poor stars have lower mass loss rates than metal-rich ones. This is explained by the fact that at lower metallicity the fraction of carbon stars is higher, and that the extinction coefficients of carbonaceous grains found around carbon stars are larger than those of silicate grains surrounding oxygen-rich stars. The latter dominate the light of the stellar populations at the same age interval at high metallicity. In addition, AGB stars become brighter as the stellar metallicity decreases, increasing the optical depth. Hence carbon stars are affected more severely than oxygen-rich stars having the same location location in the HR diagram leading to a bigger departure from the prediction of pure photospheric emission models. For older stellar populations, as the contribution of AGB stars goes down, no envelope-related metallicity effect is seen.

## 4 Conclusions

In this paper we have modelled the effects of both carbon stars and of the presence of a steady-state outflow of matter surrounding TP-AGB stars on the optical and near-IR emission of intermediate-age stellar populations. To this aim we coupled stellar evolution models and dusty shell models, in order to derive the evolution of the extinction due to circumstellar shells along the TP-AGB phase. Effects on the luminous and cool region of the HR diagram of resolved stellar populations and on integrated near-infrared properties of unresolved intermediate-age stellar populations were investigated.

The effects of both AGB sub-populations are predicted to be strong on the distribution of AGB stars in the colour-magnitude diagram. Indeed, the AGB sequence is expected to split into three morphological features occupied respectively by oxygen-rich Early-AGB/TP-AGB stars, carbon stars, and dust-enshrouded stars with little overlap, in good agreement with the feature-rich 2MASS LMC colour-magnitude diagram.

As a test of our models, we compared our isochrones with a sample of galactic M stars, Mira stars, and extreme dust-enshrouded objects in the Galaxy and in the LMC. In the near infrared colours the models reproduce the location of both optically bright and dust-enshrouded stars in a (J-H) vs. (H-K) diagram. The models predict that the dust-enshrouded carbon star population populates a sequence with redder (H-K) than dust-enshrouded M-star population at fixed (J-H).

We show that the K-band luminosity function of carbon stars is sensitive to the inclusion of the effect of dust-shells around TP-AGB stars. The peak of the luminosity function is expected to shift to fainter magnitudes and faint tails populated by cool, luminous but obscured stars close to the end of their AGB phase are predicted.

Regarding integrated properties of stellar populations, the inclusion of dust shells around AGB stars reduces the contribution of these stars to near-IR light because of the large obscuration affecting them. This will lead to bluer optical/near-IR colours in comparison to what is predicted by models where the formation of circumstellar shells is neglected. This effect decreases as the stellar population grows older.

###### Acknowledgements.

I would like to thank A. Lançon, and C. Loup for their insightful remarks on earlier version of this paper, and the anonymous referee for her/his constructive report.


format: text

text:
# Mean effects of turbulence on elliptic instability in fluids

## Abstract

Elliptic instability in fluids is discussed in the context of the Lagrangian-averaged Navier-Stokes-alpha (LANS$`\alpha `$) turbulence model. This model preserves the Craik-Criminale (CC) family of solutions consisting of a columnar eddy and a Kelvin wave. The LANS$`\alpha `$ model is shown to preserve the elliptic instability for the inviscid case. However, the model shifts the critical stability angle. This shift increases (resp. decreases) the maximum growth rate for long (resp. short) waves. It also introduces a band of stable CC solutions for short waves.

Elliptic, or tilting, instability is a fundamental phenomenon in fluids that results from parametric resonance. This is the mechanism by which vortex stretching creates three dimensional instabilities in swirling two dimensional flows. Specifically, the energy in an elliptical columnar eddy may be transferred to a propagating Kelvin wave kelvin  by this mechanism. A breakthrough in the study of this problem occurred when Craik & Criminale craik:crim:86  showed that superposing a columnar eddy and a Kelvin wave yields an exact solution of the Euler equations. Thus, the elliptical instability can be treated as a modulation of the Craik-Crimale (CC) family of solutions by using Floquet analysis as was first done by Bayly bayly:86 . The history of discovery and rediscovery of elliptic instability in laminar fluids is reviewed by Kerswell ker:02 . Here we address the mean effects of turbulence on the elliptic instability. The CC family of superposed solutions is preserved by the closure model we shall consider. For this model, turbulence is shown to enhance the growth rates of elliptic instability for Kelvin wavelengths that are longer than the turbulence correlation length. Conversely, turbulence is found to suppress elliptic instability at the shorter wavelengths and to create a band of stable CC flows with nonzero eccentricities.

The turbulence model we shall consider is the Lagrangian-averaged Navier-Stokes-alpha (LANS$`\alpha `$) model chen:foias:holm:olson:titi:wynne , whose equations are:

$$\begin{array}{c}\frac{𝐯}{t}+(𝐮)𝐯+(u)^T𝐯\hfill \\ \hfill +\left(p\frac{1}{2}|𝐮|^2\frac{1}{2}\alpha ^2|𝐮|^2\right)=\nu \mathrm{\Delta }𝐯\end{array}$$
(1)
together with $`𝐮=0`$. In this model, the mean fluid velocity $`𝐮`$ is related to the mean momentum $`𝐯`$ via the Helmholtz operator $`(1\alpha ^2\mathrm{\Delta })`$ as $`𝐯=𝐮\alpha ^2\mathrm{\Delta }𝐮`$. This Helmholtz filtering of the fluid velocity introduces the length scale $`\alpha `$ as a parameter in the model. The LANS$`\alpha `$ model preserves the fundamental transport theorems for circulation and vorticity dynamics of the NS equations. Direct numerical simulations of the LANS$`\alpha `$ model for forced homogeneous turbulence show it to be considerably less computationally intensive than the exact NS equations while preserving essentially the same behavior as NS at length scales larger than alpha chen:holm:mar:zha:99 . The unforced, inviscid Euler$`\alpha `$ form of these equations first appeared in the context of averaged fluid models holm:mar:rat . The basic properties of the LANS$`\alpha `$ model, its comparison with experiment, and its early development are reviewed in Ref. foias:holm:titi:01, . See also Refs. newalpha,  for additional results for this model. As discussed in Ref. holm:mar:rat  the LANS$`\alpha `$ turbulence equations formally coincide in the inviscid limit with a classic rheological model known as the 2nd-grade fluid grade2 . Thus, the present results for the inviscid elliptic instability apply to both the LANS$`\alpha `$ turbulence model and to rheology of 2nd-grade fluids.

We construct an exact solution to Eq. (1) with zero divergence of the form $`𝐮_0(𝐱,t)=𝒮(t)𝐱+𝐔(t)`$, where $`𝒮𝐱`$ is the action of the matrix $`𝒮`$ on the vector $`𝐱=[x,y,z]^T`$ from the left. The matrix $`𝒮`$ is a time dependent matrix with zero trace such that $`d𝒮/dt+𝒮^2=(t)`$, where $`(t)`$ is a symmetric matrix which contains the contributions of $`𝐔(t)`$. The corresponding pressure $`p_0`$ is obtained from $`(t)`$; see Ref. craik:crim:86,  for details. We nondimensionalize the system using the variables $`𝐱^{}=𝐱/l`$, $`t^{}=\omega t`$, $`𝐮^{}=𝐮/\omega l`$, $`𝐯^{}=𝐯/\omega l`$, $`\alpha ^{}=\alpha /\sqrt{l}`$, where $`l`$ is a typical length scale and $`\omega =\frac{1}{2}|\times 𝐮_0|`$. The resulting equation with the prime notation suppressed is Eq. (1) with $`\nu `$ replaced by $`\nu /\omega `$. We construct a second solution to Eq. (1) of the form $`𝐮_0+𝐮_1`$ with corresponding pressure $`p_0+p_1`$, where

$`𝐮_1`$ $`=\mu 𝐚(t)\mathrm{sin}(\beta \psi (𝐱,t)),`$ (2)
$`\begin{array}{cc}\hfill p_1& =\mu \widehat{p}_{11}(t)\mathrm{cos}(\beta \psi (𝐱,t))\hfill \\ & +\mu ^2\widehat{p}_{12}(t)\mathrm{cos}^2(\beta \psi (𝐱,t)),\hfill \end{array}`$ (3)
$`\psi (𝐱,t)=𝐤(t)𝐱+g(t)`$, and $`\mu `$ and $`\beta `$ are scaling factors so that we can choose the initial conditions $`|𝐚(0)|=1`$ and $`|𝐤(0)|=1`$. The unknown phase $`\psi (𝐱,t)`$ and the amplitudes $`𝐚(t)`$, $`\widehat{p}_{11}(t)`$, and $`\widehat{p}_{12}(t)`$ are to be determined. Such flows which are the sum of a ‘base flow’ $`𝐮_0`$ and a ‘disturbance’ $`𝐮_1`$ are called $`\alpha `$-CC flows. The incompressibility condition gives

$`𝐚𝐤=0.`$ (4)
The evolution equations for the amplitudes and phase are

$`(_t+𝒮𝐱)\psi +𝐔𝐤`$ $`=0,`$ (5)
$`\begin{array}{cc}\hfill (_t+𝒮𝐱)((1+\mathrm{\Gamma })𝐚)+\mathrm{\Gamma }& 𝒮^T𝐚+𝒮𝐚\hfill \\ \hfill (\beta \widehat{p}_{11}\beta ^2\alpha ^2𝐚𝒮𝐤)𝐤& ={\displaystyle \frac{(1+\mathrm{\Gamma })\nu \beta ^2}{\omega }}|𝐤|^2𝐚,\hfill \end{array}`$ (6)
$`\widehat{p}_{12}\mathrm{\Gamma }|𝐚|^2`$ $`=0.`$ (7)
Here $`\mathrm{\Gamma }=\alpha ^2\beta ^2|𝐤|^2`$. Note that the amplitude scaling $`\mu `$ is immaterial. The parameter $`\alpha `$ couples various terms throughout the system. Moreover, this coupling in $`\alpha `$ appears only in the combination $`\mathrm{\Gamma }`$, which is proportional to wavenumber-squared. Consequently, this coupling affects the high wavenumber behavior of the solution for $`\alpha 0`$. Equation (5) states that the phase is advected with the base flow. Only two free parameters remain ($`\mathrm{\Gamma }`$ and $`E_\omega `$) upon introducing the vorticity based Ekman number $`E_\omega =\nu \beta ^2/\omega `$. Without loss of generality, we set $`g/t+𝐤𝐔=0`$. Denoting the material derivative as $`\mathrm{d}_t=(_t+𝒮𝐱)`$ and taking the gradient of Eq. (5) reduces Eqs. (5)-(6) to a system of ordinary differential equations:

$`{\displaystyle \frac{\mathrm{d}𝐤}{\mathrm{d}t}}+𝒮^T𝐤`$ $`=0,`$ (8)
$`\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}((1+\mathrm{\Gamma })𝐚)}{\mathrm{d}t}}+(1+\mathrm{\Gamma })𝒮^T𝐚& \\ \hfill +2\varpi \times 𝐚\stackrel{~}{P}𝐤& =(1+\mathrm{\Gamma })E_\omega |𝐤|^2𝐚,\hfill \end{array}`$ (9)
where $`\stackrel{~}{P}`$ is the coefficient of $`𝐤`$ in Eq. (6), $`\varpi =\frac{1}{2}\times 𝐮_0`$ is the (normalized) vorticity of the base flow and $`(𝒮𝒮^T)𝐚=2\varpi \times 𝐚`$. We eliminate the pressure term by taking the dot product of Eq. (9) with $`𝐤`$ and by using $`\mathrm{d}𝐚/\mathrm{d}t𝐤=𝐚\mathrm{d}𝐤/\mathrm{d}t=𝒮𝐚𝐤`$, the first of which follows from Eq. (4) and the second from Eq. (8):

$`\stackrel{~}{P}={\displaystyle \frac{1}{|𝐤|^2}}\left\{(1+\mathrm{\Gamma })(𝒮+𝒮^T)𝐚𝐤+2\varpi \times 𝐚𝐤\right\}.`$ (10)
In summary, we have obtained a new exact incompressible solution to Eq. (1). The variables are amplitude $`𝐚(t)`$ and wave vector $`𝐤(t)`$. Once these are determined, the pressure terms follow from Eqs. (7) and (10). Note that $`𝐮_0`$ and $`𝐮_0+𝐮_1`$ are exact solutions to the nonlinear equations, but $`𝐮_1`$ by itself is only a solution to Eq. (1) linearized about $`𝐮_0`$. The construction described above also can be applied to Eq. (1) expressed in a rotating coordinate system in which an $`\alpha `$-CC flow can still be found. The effects of rotation will be discussed elsewhere. Finally, we emphasize that the operator $`\mathrm{d}/\mathrm{d}t+S^T`$ acting on a vector represents the complete time derivative of that quantity in a Galilean frame moving with $`𝐮_0`$.

Insight into the dynamics of the problem can be gained by examining Eq. (9) in the asymptotic regimes $`\mathrm{\Gamma }1`$ and $`\mathrm{\Gamma }1`$, where $`\mathrm{\Gamma }=\alpha ^2\beta ^2|𝐤|^2`$. (We assume that $`|𝐤(t)|`$ remains bounded and never vanishes.) For $`\mathrm{\Gamma }1`$, Eq. (9) becomes

$$\begin{array}{c}\frac{\mathrm{d}𝐚}{\mathrm{d}t}+𝒮^T𝐚=E_\omega |𝐤|^2𝐚2\varpi \times 𝐚+\frac{2𝒮𝐚𝐤}{|𝐤|^2}𝐤\hfill \\ \hfill +\mathrm{\Gamma }\left(2\varpi \times 𝐚+\frac{2}{|𝐤|^2}\left\{(𝒮𝐤𝐤)𝐚(\varpi \times 𝐚𝐤)𝐤\right\}\right)\\ \hfill +O\left(\mathrm{\Gamma }^2\right).\end{array}$$
(11)
Combined with Eq. (8), this equation preserves $`𝐚𝐤=0`$ at each order. The $`O(1)`$ term in the above equation is exactly the expression for the amplitude of the modulated traveling wave in the CC flow for the classical Euler and NS equations. This, of course, is expected since Eq. (1) reduces to the NS equations for $`\alpha =0`$. Thus, to leading order, the amplitude decays with viscosity, stretches with the base shear and rigidly rotates with the vorticity of the base flow. For $`\mathrm{\Gamma }1`$, Eq. (9) becomes

$$\begin{array}{c}\frac{\mathrm{d}𝐚}{\mathrm{d}t}+𝒮^T𝐚=E_\omega |𝐤|^2𝐚+\frac{2𝒮𝐚𝐤}{|𝐤|^2}𝐤\hfill \\ \hfill +\frac{2}{|𝐤|^2}\left\{(𝒮𝐤𝐤)𝐚(\varpi \times 𝐚𝐤)𝐤\right\}+O\left(\frac{1}{\mathrm{\Gamma }}\right).\end{array}$$
(12)
Again, this equation preserves $`𝐚𝐤=0`$ at each order Thus, as $`\mathrm{\Gamma }\mathrm{}`$ (corresponding to either $`\alpha \mathrm{}`$ or $`\beta \mathrm{}`$), the amplitude no longer rigidly rotates with the vorticity of the base flow.

As an example, we examine the stability of a rotating column of fluid with elliptic streamlines whose foci lie on the $`y`$-axis and vorticity $`\varpi =\omega \widehat{𝐞}_z`$:

$`𝐮_0=\omega L𝐱,L=\left(\begin{array}{ccc}0& 1+\gamma & 0\\ 1+\gamma & 0& 0\\ 0& 0& 0\end{array}\right).`$ (13)
Here, $`0\gamma <1`$ is the eccentricity of the ellipses, and the pressure is $`p_0=\frac{1}{2}\omega ^2(1\gamma ^2)(x^2+y^2)`$. Equation (8) with $`𝒮=L`$ is analytically solvable:

$$\begin{array}{c}𝐤=[\mathrm{sin}\theta \mathrm{cos}(t\sqrt{1\gamma ^2}),\hfill \\ \hfill \kappa \mathrm{sin}\theta \mathrm{sin}(t\sqrt{1\gamma ^2}),\mathrm{cos}\theta ]^T\end{array}$$
(14)
where $`\kappa ^2=(1\gamma )/(1+\gamma )`$ and $`\theta `$ is the polar angle that $`𝐤`$ makes with the axis of rotation. In summary, we have a four parameter problem in $`\mathrm{\Gamma }`$, $`E_\omega `$, $`\gamma `$, and $`\theta `$. Eq. (9) has the form $`\mathrm{d}𝐚/\mathrm{d}t=𝒩(t)𝐚`$, where the elements of the matrix $`𝒩(t)`$ are periodic with period $`\tau =2\pi /\sqrt{1\gamma ^2}`$. Therefore, the system can be analyzed numerically using Floquet theory yaku:star:76 . We compute the monodromy matrix $`𝒫`$, that is, the fundamental solution matrix with identity initial condition evaluated at $`t=\tau `$. Equation (9) will have exponentially growing solutions if $`\mathrm{max}_i|\mathrm{}(\rho _i)|>1`$, where $`\rho _i,i=1,2,3`$ are the eigenvalues of $`𝒫`$, with corresponding Lyapunov-like growth rates given by

$`\sigma ={\displaystyle \frac{1}{\tau }}\mathrm{ln}\{\underset{i}{\mathrm{max}}|\mathrm{}(\rho _i)|\}.`$ (15)
Thus, we can simulate numerically the solution to Eq. (9) over one period and indisputably determine the exponential growth rates. We can be certain that at least one of the eigenvalues will always be unity because of the incompressibility condition (4) and that the remaining two eigenvalues appear as complex conjugates on the unit circle or as real valued reciprocals of each other.

The present investigation considers the case of inviscid flow, i.e. $`E_\omega =0`$. Viscosity, which only slightly modifies the inviscid results, will be discussed elsewhere. For flows with circular streamlines ($`\gamma =0`$), the monodromy matrix can be analytically computed. It follows from Eq. (14) that $`|𝐤(t)|=1`$. Then, $`\mathrm{\Gamma }`$ is constant in time (denoted by $`\mathrm{\Gamma }_0=\alpha ^2\beta ^2`$) and Eq. (9) has three linearly independent solutions:

$`𝐚_1(t)`$ $`=\mathrm{cos}(\xi (t)+\varphi )𝐤_1+\mathrm{sin}(\xi (t)+\varphi )𝐤_2`$ (16)
$`𝐚_2(t)`$ $`=\mathrm{sin}(\xi (t)+\varphi )𝐤_1\mathrm{cos}(\xi (t)+\varphi )𝐤_2`$ (17)
$`𝐚_3(t)`$ $`=\widehat{𝐞}_z,`$ (18)
where $`\xi (t)=2t\mathrm{cos}\theta /(1+\mathrm{\Gamma }_0)`$, $`𝐤_2=[\mathrm{sin}t,\mathrm{cos}t,0]^T`$ and $`𝐤_1=[\mathrm{cos}\theta \mathrm{cos}t,\mathrm{cos}\theta \mathrm{sin}t,\mathrm{sin}\theta ]^T`$ are vectors orthogonal to $`𝐤`$, and $`\varphi `$ is an arbitrary phase. Clearly the first two solutions $`𝐚_1`$ and $`𝐚_2`$ satisfy Eq. (4). The monodromy matrix can be constructed from these three solutions:

$`𝒫=\left(\begin{array}{ccc}\mathrm{cos}(\xi (2\pi ))& \mathrm{cos}\theta \mathrm{sin}(\xi (2\pi ))& 0\\ \mathrm{sin}(\xi (2\pi ))/\mathrm{cos}\theta & \mathrm{cos}(\xi (2\pi ))& 0\\ \mathrm{tan}\theta (1\mathrm{cos}(\xi (2\pi )))& \mathrm{sin}\theta \mathrm{sin}(\xi (2\pi ))& 1\end{array}\right).`$
The three eigenvalues are $`\rho _{1,2}=\mathrm{exp}(\pm i\xi (2\pi )),\rho _3=1`$. All of the eigenvalues lie on the unit circle, from which it follows that all solutions in the inviscid case for $`\gamma =0`$ are stable. The values of $`\mathrm{cos}\theta `$ for which $`|\rho _i|=1,i=1,2,3`$ are called ‘critically stable’ and are given by $`\xi (2\pi )=n\pi `$, $`n=0,\pm 1,\pm 2,\mathrm{}`$, corresponding to $`\mathrm{cos}\theta =n(1+\mathrm{\Gamma }_0)/4`$. At these parameter values an exponentially growing solution can appear (together with an exponentially decaying one) as $`\gamma `$ increases from zero. Since $`\mathrm{\Gamma }_00`$, the only values of interest are $`n=0,\pm 1,\pm 2,\pm 3`$, and, for the case $`\alpha =0`$, $`n=\pm 4`$. Bayly bayly:86  argues that the evenness of $`\stackrel{~}{P}𝐤`$ as a function of $`𝐤`$ implies that the eigenvalues, if real and unequal, must be positive. This dismisses the $`n=\pm 1`$ and $`n=\pm 3`$ choices. The cases $`n=0`$ and $`n=\pm 4`$ preserve the two-dimensional structure of the base flow and thus should be stable under small perturbations in the eccentricity. The remaining value, $`\mathrm{cos}\theta =\frac{1}{2}(1+\mathrm{\Gamma }_0)`$ is the critical parameter value at which $`𝐚(t)`$ suffers exponential growth as $`\gamma `$ increases from zero. We conclude that introducing $`\alpha `$ preserves the existence of elliptic instability, though it shifts the angles at which elliptic instability arises to $`\mathrm{cos}\theta =(1+\mathrm{\Gamma }_0)/2`$. In addition, for $`\mathrm{\Gamma }_0>1`$, the LANS$`\alpha `$ model stabilizes Bayly’s elliptic instability.

Additional understanding of this result emerges by following the analysis of Waleffe wale:90  and Kerswell ker:02 . By taking the dot product of Eq. (9) with $`𝐚`$, we obtain (for all $`\gamma `$)

$`{\displaystyle \frac{\mathrm{d}\left(\frac{1}{2}|𝐚|^2\right)}{\mathrm{d}t}}=2\gamma a_1a_2+{\displaystyle \frac{4\gamma \mathrm{\Gamma }}{1+\mathrm{\Gamma }}}{\displaystyle \frac{k_1k_2}{|𝐤|^2}}|𝐚|^2.`$ (19)
One can determine an exponential growth rate to leading order in $`\gamma `$ by inserting the zeroth order solutions for $`𝐤`$ and $`𝐚_1`$ into the right hand side of this equation:

$$\begin{array}{c}\sigma \frac{1}{|𝐚|^2}\frac{\mathrm{d}\left(\frac{1}{2}|𝐚|^2\right)}{\mathrm{d}t}=\frac{\gamma }{4}[(1\mathrm{cos}\theta )^2\mathrm{sin}(2(\xi _++\varphi ))\hfill \\ \hfill (1+\mathrm{cos}\theta )^2\mathrm{sin}(2(\xi _{}+\varphi ))\\ \hfill 2(1\mathrm{cos}^2\theta )\mathrm{sin}(2t)]+\frac{2\gamma \mathrm{\Gamma }_0}{1+\mathrm{\Gamma }_0}\mathrm{sin}^2\theta \mathrm{sin}(2t),\end{array}$$
(20)
where $`\xi _\pm =\xi (t)\pm t`$. Upon averaging over a period of $`𝐚_1`$, this quantity will vanish except when $`\xi _\pm =0`$, corresponding to $`\mathrm{cos}\theta =\mathrm{\Gamma }_0/2`$. The maximum values for $`\sigma `$ will occur at $`\varphi =\pi /4`$ for $`\xi _\pm =0`$, respectively, with growth rate

$`\sigma ={\displaystyle \frac{(3+\mathrm{\Gamma }_0)^2}{16}}\gamma +O\left(\gamma ^2\right),`$ (21)
valid for $`\mathrm{\Gamma }_01`$. Thus, we see that the angle of critical stability is again $`\mathrm{cos}\theta =(1+\mathrm{\Gamma }_0)/2`$. Furthermore, we see that the maximum growth rate increases as a function of $`\mathrm{\Gamma }_0`$ due to the $`\mathrm{\Gamma }_0`$ dependence of the critical stability point up to a maximum of $`\sigma =\gamma `$ at $`\mathrm{\Gamma }_0=1`$, after which a set of stable solutions emerges in a band of nonzero eccentricities. See Fig. 1.

For nonzero values of $`\gamma `$, we must investigate the system numerically. Figure 2 shows the evolution of the critical instability surface as a function of $`\alpha ^2\beta ^2`$. For $`\alpha ^2\beta ^2<<1`$, there is little change in the critical instability surface as predicted by Eq. (11). For $`\alpha ^2\beta ^2>0`$, all angles of incidence for the traveling wave are unstable in a neighborhood of $`\gamma =1`$. The maximum growth rate in the $`\gamma `$-$`\mathrm{cos}\theta `$ plane increases as a function of $`\alpha ^2\beta ^2`$ and shifts to the corner $`\gamma =1`$, $`\mathrm{cos}\theta =1`$ by $`\alpha ^2\beta ^2=0.1`$. When $`\alpha ^2\beta ^2`$ exceeds unity, the flow stabilizes. For a given set of parameters $`(\gamma ,\mathrm{cos}\theta )`$, one of following three situations will occur as shown in Fig. 2: the flow is stable for all $`\alpha ^2\beta ^2`$; the flow is unstable for $`0\alpha ^2\beta ^2<\alpha _1^2\beta _1^2`$ and stable for $`\alpha ^2\beta ^2\alpha _1^2\beta _1^2`$; or the flow is stable for $`0\alpha ^2\beta ^2\alpha _2^2\beta _1^2`$, unstable for $`\alpha _2^2\beta _2^2<\alpha ^2\beta ^2<\alpha _1^2\beta _1^2`$, and stable again for $`\alpha ^2\beta ^2\alpha _1^2\beta _1^2`$.

Thus, the LANS$`\alpha `$ turbulence model enhances the growth rates of the elliptic instability for long waves with $`\alpha ^2\beta ^2<1`$ while it shifts the angle of critical stability along the cusp rising diagonally in Fig. 2. It also stabilizes the elliptic instability for short waves with $`\alpha ^2\beta ^2>1`$ as seen in Fig. 2b. Finally, for any $`\alpha ^2\beta ^20`$, this turbulence model modifies the region in $`(\gamma ,\mathrm{cos}\theta )`$ parameter space where the elliptic instability occurs, as also shown in Fig. 2.

In principle, one can now examine the stability properties of this new family of exact $`\alpha `$CC solutions of the LANS$`\alpha `$ model. This would be a secondary stability analysis of the rotating base flow. Work of this type was carried out by Lifschitz and collaborators lif:fab  for the classical CC flows under high-frequency, short wavelength perturbations. A similar perturbation analysis for the $`\alpha `$-CC flow will be carried out elsewhere.

The authors are indebted to A. Lifschitz-Lipton for stimulating our original interest in CC solutions and to the Center for Scientific Computing at Southern Methodist University for use of their computing resources. Furthermore, BF thanks the Theoretical Division at the Los Alamos National Laboratory for their hospitality. The solutions to $`𝐚(t)`$ were simulated using the variable coefficient ODE solver DVODE dvode .