text:
# The red tail of carbon stars in the LMC: Models meet 2MASS and DENIS observations

## 1 Introduction

It has been already known for a long time that carbon (C-type) stars occupy a particular region of near-infrared colour-colour diagrams (e.g. Richer et al. 1979; Cohen et al. 1981; Frogel et al. 1990; Hughes & Wood 1990). Specifically, C-type stars exhibit systematically redder colours than M-type stars, and particular colour-colour diagrams (e.g. $`JH`$ vs. $`HK`$) display such a sharp dichotomy between the two spectral classes that they are often used diagnostic tools to identify carbon star candidates from observed samples.

Recent wide-area near-infrared surveys like the Two Micron All Sky Survey (2MASS; Skrutskie et al. 1997) and the Deep Near-Infrared Southern Sky Survey (DENIS; Epchstein et al. 1997) have confirmed these findings. By providing photometric data for the complete AGB population of the Magellanic Clouds, they reveal the presence of a red plume of C-stars in a striking way (Cioni et al. 1999; Nikolaev & Weinberg 2000). This feature can be appreciated in the $`K_\mathrm{s}`$ vs. $`(JK_\mathrm{s})`$ diagrams of Fig. 1: they show a marked almost-vertical sequence of red giants at $`(JK_\mathrm{s})1.1`$, from which an inclined branch departs at $`K_\mathrm{s}11`$ towards redder colours, extending up to $`(JK_\mathrm{s})2.0`$. This inclined branch is what we refer to as “the red tail”. Nikolaev & Weinberg (2000) and Cioni et al. (2001) demonstrate that the red tail is populated by non-obscured C stars, whereas the almost-vertical feature at $`(JK_\mathrm{s})1.1`$ is composed basically of O-rich giants (of spectral types from late-K to M). Moreover, the luminosities of red tail stars indicate that they belong to the thermally pulsing asymptotic giant branch (TP-AGB) phase, and hence they are expected to be the result of carbon surface enrichment by the recurrent third dredge-up process.

As for the modelling of this feature, the present-day situation is as follows. Most of the stellar isochrones for old and intermediate ages available in the literature do not include the TP-AGB phase, as they usually extend up to tip of the red giant branch (RGB) or the Early-AGB (e.g. Lejeune & Schaerer 2001; Domínguez et al. 1999; Weiss & Salaris 1999; Bergbusch & VandenBerg 2001 and references therein). Two main reasons for this discrepancy are that i) fully modelling the AGB phase is indeed complex and time-consuming, and ii) most stellar models fail to explain the formation of carbon stars except for those at higher luminosities.

To our knowledge the only available isochrone sets that include TP-AGB stars are the Padova ones (Bertelli et a. 1994; Girardi et al. 2000; Salasnich et al. 2001; Marigo & Girardi 2001). In these models the entire TP-AGB evolution is followed – with the aid of a synthetic approach to various degrees of detail – up to ejection of the stellar envelope. However, all of them fail to predict the red tail of C stars, and in fact the reddest TP-AGB tracks hardly reach $`(JK)`$ colours as large as 1.3. This is illustrated in Fig. 2. If in the older isochrone sets (Bertelli et al. 1994; Girardi et al. 2000) the basic problem was the lack of the C-star phase in the stellar models, the same does not affect the most recent models (Marigo & Girardi 2001), that account for the formation of C-stars via the third dredge-up.

What is then the problem with the colours of theoretical C-star models? Why do they not reach $`(JK)`$ as red as $`2`$ ? Since in Marigo & Girardi (2001) isochrones the third dredge-up process is calibrated so that C-stars appear at the right luminosities (see Marigo et al. 1999 for details), the problem might be related either to (i) an improper transformation from $`T_{\mathrm{eff}}`$ to the $`(JK)`$ colour, which introduces errors of several tenths of a magnitude for C-type stars, and/or (ii) an improper modelling of the C-star radii and effective temperatures, with typical errors of some hundreds of degrees Kelvin.

In case of alternative (i), a simple revision of the $`T_{\mathrm{eff}}`$–colour relations, via e.g. the use of suitable relations for M- and C-type stars, would solve the problem. In the case of option (ii), the implications would be much deeper, since the $`T_{\mathrm{eff}}`$ values of Marigo et al. (1999) models are typical of most TP-AGB models in the literature.

Exploring both possibilities is the main scope of this paper, which is organised as follows. In Sect. 2 we start by recalling the main prescriptions in our synthetic TP-AGB models, that also account for the transition from O-rich (M-type) to C-rich (C-type) stars (Sect. 2.1). A brief summary of the present status of other AGB evolution models available in the literature is given in Sect. 2.2, putting particular emphasis on the role of molecular opacities.

Section 3 describes our simulations of the 2MASS data in the $`K`$ vs. $`(JK)`$ colour magnitude diagram (CMD). If several 2MASS features are reproduced in detail, the red tail of C-stars is missed by our initial models. Section 4 investigates the possible reasons for it. Finally, our main conclusions are summarised in Sect. 5.

## 2 Models for C-stars

### 2.1 Our prescriptions: old and new

In this study the formation and evolution of C-stars is described with the aid the synthetic TP-AGB model developed by Marigo et al. (1996, 1998, 1999), and Marigo (2002) to whom the reader is referred for all the details. For the sake of clarity let us briefly outline the main structure of the model and its basic prescriptions.

* The TP-AGB evolution – of a star of given initial mass and metallicity – is followed from the first thermal pulse up to the complete ejection of the envelope by stellar winds. This is performed by suitably coupling (i) analytic prescriptions (e.g. the core mass-luminosity relation, the core mass-interpulse relation, the growth rate of the core mass during the quiescent stages), and (ii) a complete static envelope model to be integrated from the atmosphere down to the core.
* The initial conditions at the first thermal pulse are derived from sets of stellar tracks for low- and intermediate-mass stars by Girardi et al. (2000), which cover the previous evolutionary phases, i.e. from the zero-age main sequence up to the beginning of the TP-AGB phase. In this way consistency and homogeneity are guaranteed in the construction of the corresponding stellar isochrones (see Sect. 3.2).
* Third dredge-up episodes possibly take place at thermal pulses and eventually lead to the transition from M- to C-class if the surface C/O ratio increases above 1. The third dredge-up is described as a function of two parameters, namely: the efficiency $`\lambda `$, and the minimum temperature $`T_\mathrm{b}^{\mathrm{dred}}`$ at the base of the convective envelope for dredge-up to take place. The previous calibration of these quantities (Marigo et al. 1999) on the observed luminosity functions of C-stars has yielded ($`\lambda =0.50,\mathrm{log}T_\mathrm{b}^{\mathrm{dred}}=6.4`$) for the LMC and ($`\lambda =0.65,\mathrm{log}T_\mathrm{b}^{\mathrm{dred}}=6.4`$) for the SMC.
* Hot-bottom burning occurring in the most massive AGB stars ($`M4.5M_{}`$), is followed in detail with the aid of complete envelope integrations (Marigo et al. 1998, Marigo 1998). The possible break-down of the core mass-luminosity relation and the CNO-cycle nucleosynthesis (possibly preventing the formation of luminous C-stars) are both taken into account.
* Mass loss is included according to Vassiliadis & Wood’s (1993) formalism, based on observations of pulsating AGB stars. In our calculations we assume that pulsation occurs with periods ($`P`$) corresponding to either the fundamental mode ($`P_0`$), or the first overtone mode ($`P_1`$).
* Of great importance for the analysis developed in this work are the opacity prescriptions to be used in the envelope model. At high temperatures ($`T10000`$ K) we adopt Iglesias & Rogers (1993) for O- and C-rich mixtures. At low temperatures ($`T<\mathrm{10\hspace{0.17em}000}`$ K) we use either the solar-scaled opacity tables by Alexander & Ferguson (1994; hereafter also $`\kappa _{\mathrm{fix}}`$), or the routine developed by Marigo (2002; hereafter also $`\kappa _{\mathrm{var}}`$) that predicts molecular opacities for any chemical composition of the gas. It should be remarked that before Marigo (2002), all TP-AGB tracks by Marigo et al. were calculated with $`\kappa _{\mathrm{fix}}`$.

In summary, for the purpose of our present work we will make use of two groups of synthetic TP-AGB calculations, namely:

* Old TP-AGB models with $`\kappa _{\mathrm{fix}}`$ (Marigo et al. 1999);
* New TP-AGB models with $`\kappa _{\mathrm{var}}`$ (Marigo 2002).

For the remaining prescriptions, we refer to the aforementioned model outline unless otherwise specified.

### 2.2 The formation of C-stars in other models and the role of molecular opacities

Before proceeding with our analysis it is worth recalling briefly the current status of AGB models available in the literature and the importance of a consistent choice in the use of molecular opacities.

First of all, the lack of efficient dredge-up in low-mass AGB models (say with $`M2M_{}`$) has been a difficulty shared by all full stellar evolution codes for a long time (the so-called “carbon star mystery”, as referred to by Iben (1981)). More recently this difficulty seems to be overcome by a few groups (e.g. Herwig et al. 1997; Straniero et al. 1997; Karakas et al. 2002), due to different numerical treatments in describing the convective boundaries.

However, even if some progress is achieved by full stellar evolution models in relation to the formation of low-mass C-stars, the distance between theory and observations is still large. Currently there is no complete set of evolutionary tracks – based on full stellar models – extending up to the end of the AGB phase for a sufficient coverage of stellar masses and metallicities. This prevents a systematic test of stellar models by comparison with observations, e.g. the luminosity functions of M- and C-stars, and simulations of CMDs including the most evolved AGB stars. The main reason is the heavy computational effort required to fully model the AGB phase. On the other hand, this limit can be totally removed if one allows for a certain loss of detail and opts for the more agile synthetic approach (see e.g. Groenewegen & de Jong 1993, Marigo et al. 1996, 1998 for recent works).

Another heavy inadequacy of present stellar models (both full and synthetic) that has been recently highlighted by Marigo (2002) is the use of fixed solar-scaled molecular opacities in AGB models experiencing the third dredge-up.

In fact, this prescription (previously defined as $`\kappa _{\mathrm{fix}}`$ case) is still the standard choice common to most full and synthetic evolutionary calculations111Actually, non-solar scaled molecular opacities were already taken into account in a simple way by Vassiliadis & Wood (1993), but no significant difference emerged in comparison to models with solar-scaled opacities, as the authors did not practically find carbon dredge-up in their AGB calculations. of the AGB phase (e.g. Mouhcine & Lançon 2002; Chieffi et al. 2001; Herwig 2000; Ventura et al. 1999; Wagenhuber & Groenewegen 1998; Forestini & Charbonnel 1997; Straniero et al. 1997).

Marigo (2002) has proved, instead, that properly coupling the molecular opacities to the actual surface chemical abundances ($`\kappa _{\mathrm{var}}`$) brings along such important consequences that the standard evolutionary scenario for the AGB evolution may be significantly affected. The reader is referred to that work for a detailed analysis.

One major effect is that the surface enrichment in carbon that leads to C/O$`>1`$ – i.e. the formation of carbon stars – causes a significant increase in opacities, which in turn is responsible for a sudden and marked cooling of the stellar tracks in the H-R diagram. The sizeable excursion towards lower $`T_{\mathrm{eff}}`$, displayed by the $`\kappa _{\mathrm{var}}`$ C-star tracks, sets the first significant difference compared to $`\kappa _{\mathrm{fix}}`$ C-star models that, instead, do not show any change in the slope of their evolutionary tracks while climbing the AGB at increasing luminosities.

A closer comparison between models with $`\kappa _{\mathrm{fix}}`$ and $`\kappa _{\mathrm{var}}`$, and the impact of the new opacities prescription will be discussed in Sect. 4.2.

## 3 The $`K`$ vs. $`(JK)`$ diagram towards the LMC

In this section, we aim to model the red tail of carbon stars in the LMC by using the TP-AGB models at our disposal. Our study is limited to a particular CMD – the $`K`$ vs. $`(JK)`$ one – that reveals the red tail in all its prominence (see Fig. 1). Needless to say, the same feature is displayed in a variety of different CMDs involving near-infrared passbands, as can be seen e.g. in van der Marel & Cioni (2001, their figure 2).

### 3.1 The selected data

In the following, we base our comparisons on 2MASS data towards the LMC, that present an extremely good photometric quality. More specifically, we use the point-source catalog of the 2MASS Second Incremental Data Release (see Cutri et al. 2002), which includes $`10\sigma `$ detections as faint as $`16.3`$, $`15.3`$, $`14.7`$ mag in $`J`$, $`H`$ and $`K_\mathrm{s}`$, respectively. As it can be seen in Fig. 3, these data wholly cover the LMC except for two limited strips of right ascension. In order to reduce the fore-/back-ground of Galactic stars and distant galaxies, we select only the stars located closer than 4 degrees (projected) to the LMC bar optical center (at $`\alpha =5^\mathrm{h}24^\mathrm{m}`$, $`\delta =69^{}44\mathrm{}`$, cf. Nikolaev & Weinberg 2000). This region, also considering the empty strip, translates into to a total sky area of about 48.2 squared degrees. The corresponding CMD is shown in Fig. 4.

### 3.2 Simulations

The $`K`$ vs. $`(JK)`$ diagram has been simulated by using a population synthesis code (Girardi et al., in preparation). In short, the code randomly generates stars following a given star formation rate (SFR), age-metallicity relation (AMR) and initial mass function (IMF). The stellar intrinsic properties (luminosity $`L`$, effective temperature $`T_{\mathrm{eff}}`$, surface gravity $`g`$, etc.) are interpolated over a large grid of stellar evolutionary tracks, based on Bertelli et al. (1994) for massive stars, Girardi et al. (2000) for low- and intermediate-mass stars, and complemented with the TP-AGB tracks from either Marigo et al. (1999) with $`\kappa _{\mathrm{fix}}`$, or new calculations with $`\kappa _{\mathrm{var}}`$ made on purpose (refer to Marigo 2002; see also Sect. 2).

The evolution along the TP-AGB phase is initially described in terms of the properties in the quiescent stages between thermal pulses. Actually, there are two additional factors that cause significant excursions in luminosity and effective temperature from their quiescent values, namely:

* occurrence of thermal pulses;
* instability against pulsation.

The former effect is included in our simulations. In practice, the luminosity variations driven by a thermal pulse (both the fast luminosity peak and long-lasting low-luminosity dip) are described according to the luminosity distributions over a pulse cycle as derived from Boothroyd & Sackmann (1988), and assuming that the envelope mass is the main factor determining the shape of the luminosity dip. Once a star is randomly scattered in $`\mathrm{log}L`$ according to this distribution, the corresponding excursion in effective temperature is determined from a grid of envelope integrations performed in precedence. In this way, the intrinsic dispersion of TP-AGB stars in the HR diagram is realistically simulated with just a modest computational cost. We recall that this dispersion is important mainly for low-mass AGB stars, with a maximum displacement of $`0.4`$ in $`\mathrm{log}L`$ (or $`+1`$ mag), affecting up to 30 percent of their TP-AGB lifetimes.

As for the latter point related to stellar pulsation, it is worth premising the following. We recall that both 2MASS and DENIS provide single-epoch observations. Hence, the sampling of variable stars (e.g. Miras, SR variables) at a single and random phase of their variability cycle causes a further dispersion of data points in the observed CMD. This effect has not been simulated here for two reasons. First, we do not know exactly how $`L`$ and $`T_{\mathrm{eff}}`$ vary during the pulsation cycle. Second, in any case the effect of this cyclical variation is easy to foresee: the bulk of TP-AGB stars will be scattered in the CMDs, to both higher and lower magnitudes and colours. The expected result is a sort of general blurring of the TP-AGB sequences in the CMD. Mira and SR variables have pulsation amplitudes ranging from 0.4 to 1.0 mag in the $`I`$-band (Hughes & Wood 1990). Therefore, to take this effect into account the simulated TP-AGB sequences should be blurred by just some tenths of magnitude.

Once the stellar intrinsic properties are singled out, the photometry is generated by applying the extended tables of bolometric corrections from Girardi et al. (2002) for O-rich stars (with C/O$`<1`$), and empirical relations – to be discussed in Sect. 4.1 below – for C-type stars (with C/O$`>1`$).

Our simulations are computed in the $`JHK`$ filter system defined by Bessell & Brett (1988). The expected differences from the 2MASS and DENIS systems – that use a “$`K`$-short” filter, $`K_\mathrm{s}`$ – are probably very low, of less than say 0.1 mag in $`(JK)`$. This is illustrated in the case of 2MASS by Cutri et al. (2002, section II,2d).

We simulate the photometric errors by means of a Monte Carlo approach, adopting the distribution of error values tabulated from 2MASS. We recall that $`\sigma `$ errors are lower than $`0.03`$ mag for $`K_\mathrm{s}13`$.

Finally, the population synthesis code generates simulations for both the LMC galaxy and the Galactic foreground, by using different distributions of reddening and distances. Results are shown in Fig. 5, for the specific case of TP-AGB models computed with solar-scaled molecular opacities ($`\kappa _{\mathrm{fix}}`$), and discussed in the following.

#### 3.2.1 The Galactic foreground

The field foreground stars have been included with the aid of a complete Galaxy model (Girardi et al., in preparation; Groenewegen et al. 2002). It consists of disk and halo components:

The disc is described by a double-exponential density law, with a radial scale length (in Galactocentric distance) of 2.8 kpc, and a scale height $`H`$ that increases with stellar age $`t`$ as

$$H=z_0(1+t/t_0)^\alpha .$$

We adopt the parameters $`z_0=95`$ pc, $`t_0=0.5`$ Gyr and $`\alpha =1.66`$, that well describe the disk components as derived by Ng et al. (1997). It is worth noticing that this prescription allows for a significant increase of $`H`$ at large ages, which is similar to assuming the presence of an old thick disk component. The distributions in mass, age and metallicity of the simulated disk stars are derived under the assumptions of: a constant star formation rate between 12 Gyr ago and now; the empirical age-metallicity relation from Rocha-Pinto et al. (2000) with the addition of a dispersion of 0.2 dex in $`[\mathrm{Fe}/\mathrm{H}]`$ at any age; and the log-normal IMF from Chabrier (2001).

The halo is described by an oblate spheroidal (cf. Gilmore 1984) of axial ratio $`q=0.65`$, and a core radius of 2.8 Kpc. Halo stars are assumed to be old (from 12 to 13 Gyr), and to follow a Gaussian distribution of metallicities of mean $`[\mathrm{Fe}/\mathrm{H}]=1.6`$ and standard deviation of 1.0 dex.

This set of prescriptions has been constrained by Groenewegen et al. (2002), by using a set of Deep Multicolor Survey (Osmer et al. 1998) and ESO Imaging Survey data (Prandoni et al. 1999). Interestingly, most of the calibration data were in visible pass-bands, whereas the model is now being applied to the interpretation of near-infrared data.

For the purpose of our study, we simulate a conic section of the Galaxy occupying the same sky area as in the selected 2MASS data. The field center has galactic coordinates $`\mathrm{}=280^{}.46`$, $`b=32^{}.89`$. We recall that, for this particular line-of-sight, other Galactic components (spiral arms, the bulge) are of no relevance.

An important aspect is that in the modelling we can distinguish different kinds of stars and hence tackle the origin of the several observed features in the CMD of 2MASS.

Our field simulation is presented in (Fig. 5a). In this plot, we can identify three marked vertical sequences of foreground disk stars (green dots), located at $`(JK)`$ colours of about 0.35, 0.65, and 0.9. Their origin is easily understood by looking at the corresponding plot in absolute magnitude (Fig. 6): The sequence to the left (at $`(JK)0.35`$) is defined mainly by the old disk turn-off, with typical masses of $`M0.9`$ $`M_{}`$; the central one (at $`(JK)0.65`$) hosts fainter RGB and red clump stars; and the right one (at $`(JK)0.9`$) corresponds to low-mass dwarfs with $`M0.6`$ $`M_{}`$. The vertical development of these strips in the CMD (Fig. 5a) simply reflects the large range of distances involved.

We notice that disk stars younger than $`4`$ Gyr ($`M1.3`$ $`M_{}`$) do not populate significantly this diagram, since they are found preferentially at small scale heights and hence at lower galactic latitudes and brighter apparent magnitudes (i.e. $`K10`$). The contribution of halo stars (black points) is also of minor importance, accounting for just a handful of points in the synthetic CMD, mostly at fainter apparent magnitudes (i.e. $`K13`$), due to their large mean distances.

Finally, an important point is that field stars are not expected to reach colours redder than $`(JK)1.0`$, so that they hardly contaminate the CMD features produced by the LMC population of AGB stars. Their intepretation is our next goal.

#### 3.2.2 The LMC population

The LMC population has been initially described by means of a very simplified model, under the assumptions of i) a fixed distance of 52.3 Kpc (Girardi & Salaris 2001), ii) constant SFR from 15 Gyr ago to now, and iii) a constant metallicity, $`Z=0.008`$, at any age. The latter choice is made for the sake of simplicity, since it requires the use of a single set of evolutionary tracks, hence relieving us, for the moment, from the effort to compute complete sets of TP-AGB models for several metallicities222Extended sets of TP-AGB tracks and isochrones for a wide range of metallicities are being computed and are soon to become available.. Anyway, the quality of our results should not depend much on these assumptions333In one of our test cases, we used TP-AGB tracks interpolated within an extended grid of 5 different metallicities, together with a more realistic AMR (from Pagel & Tautvaisiene 1998). No significant difference was noticed in the resulting CMDs..

In order to get a number of LMC giants comparable to the observed sample (Fig. 4), we have to simulate a total star formation of $`5\times 10^8`$ $`M_{}`$ during the entire LMC history. Of course, this is just a rough estimate of the mass consumed to form stars, its value depending on the adopted assumptions about SFR (constant), IMF (from Kroupa 2001), and size of the simulated region (i.e. within the innermost $`4^{}`$ of the LMC field). Anyway, it is interesting to notice that this number is just an order of magnitude lower than the total LMC mass (including dark matter; see van der Marel et al. 2002), and comparable to the present LMC H i mass (Staveley-Smith et al. 2003).

The simulated LMC population draws a few almost vertical features in the CMD (Fig. 5b), mainly populated by: young and intermediate-age main-sequence stars (with $`(JK)<0.5`$), intermediate-mass stars in the stage of core-He burning (the long diagonal sequence at $`(JK)0.8`$), intermediate-age and old RGB stars of larger luminosities (the prominent sequence at $`(JK)1`$ and with $`K>12`$), and intermediate-mass stars during the early-AGB phase (the weaker sequence departing slightly from the left of the RGB-tip, with the same inclination as the RGB). All these features (marked with blue dots) correspond to O-rich stars in evolutionary stages previous to the TP-AGB.

The O-rich TP-AGB stars (marked with cyan dots) trace two well-defined sequences: the first consists of intermediate-mass stars that somewhat extend the plume of early-AGB stars towards higher luminosities; the second defines a clear strip located directly above the RGB-tip, corresponding to low-mass stars. Finally, The C-rich TP-AGB stars (marked with red dots) distribute above and to the right of these latter sequences.

Combining both panels of Fig. 5, it turns out that our initial simulations explain quite well the several “fingers” displayed by the 2MASS data of Fig. 4, but for one major discrepancy: the location of C-stars, predicted by these particular models, simply align above the sequence of LMC O-rich giants, failing to form the red tail. As for the other CMD features, the agreement with data is indeed quite good, and our interpretations substantially agree with those advanced by Nikolaev & Weinberg (2000).

## 4 Probing the red tail of C-stars

In the following we will focus on the region of our simulated diagrams populated by TP-AGB stars, with the aim of casting light on the main cause that produces the observed red tail of C-stars in near-infrared colors. We investigate two possible factors separately, namely: the colour–$`T_{\mathrm{eff}}`$ relations (derived from empirical calibrations), and the evolution in effective temperatures (derived from stellar models).

### 4.1 Testing different colour–$`T_{\mathrm{eff}}`$ relations

Let us start by analysing the effect on our simulations of different $`(JK)T_{\mathrm{eff}}`$ transformations. It is important to remark that all computations discussed in this section are carried out by adopting TP-AGB tracks for solar-scaled $`\kappa _{\mathrm{fix}}`$ opacities. Results are shown in Fig. 7.

As a first step, in Sect. 4.1.1 we consider the sharp dichotomy between colour–$`T_{\mathrm{eff}}`$ relations derived from observed M-type stars (e.g. Fluks et al. 1994) and C-stars (e.g. Bergeat et al. 2001). Such differences essentially reflect the abrupt change in the dominant molecular species – producing the line blanketing – when passing from an O-rich spectrum (e.g. TiO, VO, H2O) to a C-rich one (e.g. C2, CN, SiC)

Then, within the class of C-stars, in Sect. 4.1.2 we attempt to account, in a simple way, for the possible composition dependence (related to the extent of carbon enrichment) of colour–$`T_{\mathrm{eff}}`$ relations.

#### 4.1.1 A $`T_{\mathrm{eff}}`$$`(JK)`$ relation valid for Galactic C-stars

Figure 7a shows our initial simulation: The TP-AGB models are calculated with $`\kappa `$-fix molecular opacities, and the conversion from $`T_{\mathrm{eff}}`$ to $`(JK)`$ is based on Fluks et al. (1994) empirical spectra and temperature scale for nearby M giants (see Girardi et al. 2002 for details). This means that a relation valid for O-rich stars has been applied also to C-stars. The discrepancy with observed data is obvious: C-stars just align above the sequence of O-rich stars, and do not develop any hint of red tail.

As a first attempt to include proper relations valid for C-stars, we apply the Bergeat et al. (2001) formulas, that give $`K`$-band bolometric corrections and $`(JK)`$ as a function of $`T_{\mathrm{eff}}`$. These fitting relations are derived from a sample of Galactic C-stars. In this context we would like to note that the Bergeat et al. (2001) relationship is essentially the same as that of Bessell et al. (1981). The resulting simulations are shown in Fig. 7b. With respect to the previous case, the sequence of C-stars is shifted to the blue, and not to the red as required to explain the data.

In this and other plots we show next, we have also adopted the empirical $`K`$-band bolometric corrections for C stars presented by Frogel et al. (1990). They differ little (by $`0.1`$ mag) from those independently derived by Bergeat et al. (2001), and yield closely similar simulated CMDs.

#### 4.1.2 The effect of molecular blanketing for C-stars

Molecular blanketing is responsible for the distinction between M- and C-spectral features and also affects the infrared colours of C-stars themselves (e.g. Cohen et al. 1981; Loidl et al. 2001). In this latter case, the extent of molecular blanketing is essentially related to the amount of carbon available to form molecules other than CO.

In their study of Galactic and LMC C-stars, Cohen et al. (1981) pointed out that the observed mean $`(JK)`$ colour increases steadily as the carbon abundance class – defined by the C2 band strength – increases (from say C,2 to C,5). Moreover, the authors showed that, for a given effective temperature, redder $`JHK`$ colours are predicted by model atmospheres with larger C2 and CN enhancements.

In our work we attempt to account, in a simple way, for the effect of molecules on the infrared colours of C-stars, by introducing a C/O-dependence into the empirical $`(JK)T_{\mathrm{eff}}`$ relation by Bergeat et al. (2001). The C/O data for Galactic C-stars are taken from Lambert et al. (1986). Figure 8 displays the observed trends of the $`(JK)`$ colour as a function of $`T_{\mathrm{eff}}`$ and C/O ratio. The linear bi-parametric fit relation to the data is:

$$(JK)=17.324.56\mathrm{log}T_{\mathrm{eff}}+0.052\mathrm{C}/\mathrm{O}$$
(1)
We see that the $`(JK)T_{\mathrm{eff}}`$ relation (top-left panel of Fig. 8) is quite narrow for Galactic C-stars, and it is well reproduced by our fitting formula for C/O$`=1.1`$, in agreement with the observational finding that most Galactic carbon stars have low C/O ratios, not much larger than one (Lambert et al. 1986). Moreover, despite some observed scatter and the small number of data points, we also note that there is a positive correlation between the $`(JK)`$ colour and the C/O ratio (bottom-left panel of Fig. 8). This trend is included in our fitting relation as it predicts a systematic shift towards redder $`(JK)`$ colours for e.g. C/O$`=3.0`$, compared to the case C/O$`=1.1`$ (top left panel).

Overall, Eq. (1) fits the data to better than 0.1 mag in $`(JK)`$, as indicated by the residuals in the right panels of Fig. 8.

Actually, in a recent study Matsuura et al. (2002) have pointed out that the molecular features in the spectra of LMC C-stars are consistent with C/O ratios larger than those of Galactic C-stars. This empirical finding completely supports the results of our TP-AGB calculations, as shown in Fig. 9. The same figure indicates that, in order to model the LMC population, we will have to apply Eq. (1) to C/O values as high as 4, that is beyond the maximum observed values of Galactic C-stars.

Thus, by using Eq. (1) in our simulations, we get the synthetic CMD of Fig. 7c. As it can be noticed, the dependence of $`(JK)`$ on the C/O ratio now causes the sequence of C-stars to be significantly redder, reaching $`(JK)=1.7`$. However, the C-stars are still largely distributed among and above the O-rich giants, and no clear red tail is drawn.

If we limit consideration to models with C/O$`<2`$, i.e. plotting only those stars for which Eq. (1) is strictly valid, we get Fig. 7d. In this case, we can notice that the low-C/O C-stars contaminate significantly the CMD region corresponding to O-rich AGB stars. On the contrary, such a high fraction of C-stars mixed up with M-stars is not present in observed samples. It is clear that this discrepant feature would remain for whatever hypothetical $`(JK)T_{\mathrm{eff}}`$ C/O that could push somehow the models to reach $`(JK)`$ colours as high as $`2.0`$.

In conclusion, our test calculations seem to indicate that the principal cause of the failure of models with solar-scaled molecular opacities in reproducing the red tail of carbon stars should not be ascribed to the $`(JK)T_{\mathrm{eff}}`$ relation.

### 4.2 Testing different TP-AGB tracks

If we assume that Eq. (1) reasonably well describes the $`(JK)`$ colours of C-stars, as supported by the low fitting residuals (of $`0.1`$ mag, see the right panels of Fig. 8), we are left with just one alternative: redder $`(JK)`$ colours could be attained only by assigning lower $`T_{\mathrm{eff}}`$ to C-star models. Moreover, in order to locate the C-star models onto a different branch from that of the O-rich models, a net $`T_{\mathrm{eff}}`$ segregation between the two classes turns out to be necessary. These requirements are indeed fulfilled by Marigo (2002) TP-AGB models with variable molecular opacities.

Figure 10 illustrates a set of $`Z=0.008`$ isochrones based on Marigo’s (2002) models. Similar isochrone sets have been computed for a large variety of dredge-up parameters and metallicities, and will be extensively described in a future paper. The conversion to $`(JK)`$ colours for C-star models is now based on Eq. (1). It is also interesting to notice that in this case we could not even have applied the same colour transformations derived for M-stars, since the effective temperatures of C-star models with $`\kappa _{\mathrm{var}}`$ are already much cooler than the validity lower-limits of such relations.

#### 4.2.1 The sole effect of variable molecular opacities

The effect of TP-AGB models with variable molecular opacities is illustrated in the simulations of Fig. 11. TP-AGB models are computed with different assumptions (see also Sect. 2) regarding the onset of the third dredge-up – controlled by the parameter $`T_\mathrm{b}^{\mathrm{dred}}`$ –, and the onset of the super-wind regime – depending on the pulsation mode, either fundamental ($`P=P_0`$), or first overtone ($`P=P_1`$). In all cases, the sequence of C-stars clearly departs from the branch of O-rich stars.

In these $`\kappa _{\mathrm{var}}`$ models the excursion towards redder $`(JK)`$ colours is mostly caused by the displacement towards cooler $`T_{\mathrm{eff}}`$. Instead, the increase of the C/O ratio during the TP-AGB evolution should play a smaller role. In fact, considering that the maximum C/O ratio reached by these models is of about 4 (Fig. 9), the C/O-dependence in the colour transformation of Eq. (1) accounts for at most $`0.2`$ mag of the total colour excursion at the end of the C-star evolution, the mean value being of about 0.1 mag.

In Fig. 11a, we adopt the same dredge-up parameters as calibrated in Marigo et al. (1999), i.e. $`\mathrm{log}T_\mathrm{b}^{\mathrm{dred}}=6.4`$ and $`\lambda =0.5`$. The simulated C-star red tail resembles very much the observed one, but for two aspects.

First, the red tail shows up at a too low luminosity. It starts at $`K=11.6`$ when $`(JK)=1.2`$, whereas the observed tail is at least 0.6 mag brighter. This problem can be easily solved with a re-calibration of the dredge-up parameters. To give an example, Fig. 11b shows a simulation that makes use of TP-AGB models recomputed with $`\mathrm{log}T_\mathrm{b}^{\mathrm{dred}}=6.45`$. This parameter choice has the effect to delay the onset of the third dredge-up to higher stellar core masses and hence higher luminosities. As we can see, now the red tail starts at $`K=11.0`$ when $`(JK)=1.2`$, which is about the right luminosity level to fit the data.

It is worth remarking that in both simulations presented in the upper panels (a and b) of Fig. 11, we have assumed that all TP-AGB stars are fundamental-mode pulsators, meaning that the Vassiliadis & Wood’s (1993) mass-loss formalism is applied with $`\dot{M}(P)=\dot{M}(P_0)`$. This assumption introduces the second problem of the models, which is discussed below.

#### 4.2.2 Extending the red tail with first-overtone pulsators

The mentioned second theoretical difficulty resides in the extension of the red tail, which in both simulations (Figs. 11a and b) clearly reaches $`(JK)=1.8`$, but hardly $`(JK)=2.0`$ as in the observations. Although this point could be considered as a minor discrepancy (at least compared with the $`\kappa _{\mathrm{fix}}`$ cases illustrated in Fig. 7), it deserves some attention. There are at least two possible ways out. The first is that we have not considered the effect of circumstellar reddening on the colours of C-stars. This should be ascribed to thick circumstellar material ejected by the coolest C-stars during an extreme mass-loss phase. In fact, a few LMC C-stars are known to extend up to $`(JK)6`$ (Cioni et al. 1999; Nikolaev & Weinberg 2000); they are among the so-called “obscured” or “dust-enshrouded” C-stars identified with IRAS sources (see van Loon 1999), which are found to undergo heavy mass-loss. The question, however, is if this sort of self-reddening is already effective in red tail C-stars with say $`(JK)_01.8`$, being able to shift their colours up to $`(JK)=2.0`$. This point has to be investigated with the aid of more detailed modelling.

The second explanation stands on the assumption that the dominant pulsation mode of AGB stars is the first overtone (with period $`P=P_1`$). This would allow us to get simulated colours further increased by 0.2 mag, while keeping the same $`T_{\mathrm{eff}}`$ vs. $`(JK)`$ relation. This is illustrated in the bottom panels of Fig. 11. For these simulations we use TP-AGB models computed with the same dredge-up parameters as in the upper panels, but assuming that Mira and SR variables are first overtone pulsators. Under this hypothesis, variable AGB stars should have shorter pulsation periods for given stellar $`L`$ and $`T_{\mathrm{eff}}`$, compared to the fundamental-mode counterparts.

As a consequence, by virtue of the empirical positive correlation between the mass-loss rates and the periods of Mira stars, we expect that for the shorter first overtone periods ($`P_0/P_12.2`$) the onset of the super-wind regime is postponed, thus implying that the end of the AGB evolution is delayed. We quantify this effect by adopting the Vassiliadis & Wood’s (1993) formalism with $`\dot{M}(P)=\dot{M}(P_1)`$ (Figs. 11c,d). In this way we find that C-star models are allowed to reach slightly higher $`L`$ and lower $`T_{\mathrm{eff}}`$, and their lifetimes increase accordingly. This latter effect is also relevant if we consider that longer lifetimes are needed to better reproduce the observed C-star counts in LMC clusters (Girardi & Marigo 2003, and work in preparation).

In Fig. 11c, we can notice that now the red tail actually extends up to $`(JK)=2.0`$. However, it still appears at too faint luminosities, due to the low value assumed for $`T_\mathrm{b}^{\mathrm{dred}}`$ (see remarks in Sect. 4.2.1). Instead, for a larger $`T_\mathrm{b}^{\mathrm{dred}}`$, the consequent later onset of the third dredge-up is able to shift the entire red tail to the right luminosities, as shown in Fig. 11d.

From a comparison with our reference plot of Fig. 4, we can conclude that the simulation of Fig. 11d best reproduces the observed morphology of the red tail. On the other hand, the introduction of the first overtone hypothesis produces an undesirable effect in the sequence of old O-rich stars: they now extend to much higher luminosities (up to $`K=10.4`$ and $`K=10.0`$ in Figs. 11c and d, respectively) than observed (up to $`K=11.2`$ in Fig. 4).

As a matter of fact, nowadays the observational scenario of long-period variables (LPVs) appears quite complex, mainly due to the observed mixture of different pulsation modes, and the lively controversy about their assignment (see e.g. Wood et al. 1999, Feast 1999, Whitelock & Feast 2000, Cioni et al. 2001, Noda et al. 2002). The situation gets even more intricate if one considers the evidence that some LPVs are found to be multi-period pulsators (Bergeat et al. 2002), and the possibility of an evolutionary path across different period-luminosity relations, so that LPVs may switch pulsation modes while evolving along the AGB (Cioni et al. 2001).

Without entering the open debate on the dominant pulsation mode, we perform a very simple test on the basis of purely theoretical considerations. Considering that the red tail of C-stars is somewhat better reproduced with the first-overtone assumption which, instead, does not seem to suitably describe the O-rich branch, we assume a mixed population of fundamental-mode and first-overtone pulsators, with the following composition:

1. All stars with mass lower than 1.3 $`M_{}`$, and larger than 3.5 $`M_{}`$, are fundamental-mode pulsators just because they do not evolve into a C-star phase;
2. The remaining stars are assigned either as fundamental-mode or as first-overtone pulsators. We assume a mix of 50 percent between the two types.

Assumption (2) may be not realistic enough, because the simulated stars do not switch between pulsation modes as they evolve, but they can just follow one between two different evolutionary paths. If this is an oversimplification of the problem – unavoidable because the mechanisms leading to change the pulsation mode cannot be presently modelled – at the same time it is a resonable and instructive experiment. In Fig. 12 we present the resulting simulation, compared with the selected 2MASS data. We can notice, in fact, that the resulting sequence of O-rich giants is noticeably shortened in luminosity, as expected.

We may therefore conclude that now our final simulation actually meets 2MASS and DENIS observations, as all CMD features are well reproduced by the models, including the C-star red tail and the plume of O-rich low-mass TP-AGB stars.

## 5 Concluding remarks

In this paper, we have shown that most features seen in the 2MASS and DENIS $`K_\mathrm{s}`$ vs. $`(JK_\mathrm{s})`$ diagram are quite well described by our present stellar models. In particular, the observed tail towards the reddest $`(JK_\mathrm{s})`$ colours provides a direct and stringent probe to the evolutionary theory of TP-AGB stars. In this regard, we have pointed out that:

1. The red tail of carbon stars is missed by models in which the TP-AGB evolution is followed using solar-scaled molecular opacities.
2. This failure can hardly be caused by some hypothetical inadequacy in the present-day $`T_{\mathrm{eff}}`$ vs. $`(JK)`$ relations for C-stars. Our tests indicate that in order to reach the observed $`(JK)`$ colours of C-stars with “solar-scaled opacity” models, quite high mean values of the C/O ratio should be exhibited by the bulk of carbon stars, which seems unrealistical.
3. Instead, the red tail comes out naturally of TP-AGB models which are calculated with variable molecular opacities (cf. Marigo 2002), i.e. opacities consistently computed for the current chemical composition of the stellar envelope. This can be understood considering that, as the third dredge-up makes the C/O ratio increase above 1, drastic changes in the main sources of molecular opacities cause both the transition from M- to C-spectral types, and a significant decrease of the effective temperature. The latter aspect has been so far largely ignored in AGB evolution models. Contrary to what claimed by Mouhcine & Lançon (2002), the present study demonstrates that in order to model properly the properties of AGB, in particular C-rich, stellar populations, one cannot ignore the large opacity effects have on the effective temperature of C-stars, and more generally on the AGB evolution.
4. Our simulations with “variable opacity” models provide a few additional useful indications about the TP-AGB evolution. For instance, models with $`\mathrm{log}T_\mathrm{b}^{\mathrm{dred}}=6.45`$ seem to reproduce well the $`K`$-band luminosity of red tail C-stars. This means that the onset of the third dredge up should occur at slightly larger core masses, as compared to previous models by Marigo et al. (1999).
5. Moreover, in the context of the present models, based on the use of Vassiliadis & Wood (1993) prescription for mass loss and an empirical $`T_{\mathrm{eff}}`$ vs. $`(JK)`$ relation, we may derive two indications about the pulsation modes of variable AGB stars. First, the colour extension of the red tail may imply a significant fraction of first-overtone pulsators among C-type stars. However, the same extension could be reproduced with fundamental-mode pulsators if circumstellar reddening starts to be important for C-stars with $`(JK)_01.8`$, a point that has still to be properly investigated. Second, the luminosity extension of the sequence populated by O-rich low-mass TP-AGB stars would suggest that most of M-type variables are fundamental-mode pulsators.

All these indications are being used to calibrate our TP-AGB models, from which extended sets of theoretical isochrones at varying metallicity will be derived and soon made available. Of course, new models will be computed for the variable-opacity case only.

###### Acknowledgements.

We would like to thank our referee Dr. M. Bessell, and P.R. Wood, for useful remarks that improved the final version of the paper. We thank M. Groenewegen for his fine calibration of the parameters used in the Galaxy model and for useful comments. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. The DENIS data is the result of an extended international collaboration with the use of ESO telescopes. We acknowledge financial support from the MIUR-COFIN 2001 project “The tridimensional structure of the Galaxy” (number 2001028112).


format: text

text:
# Power-law behavior in signal scattering process in vertical granular chain with light impurities

## 1 Introduction

In the last ten years physicists have been interested again in the study of dry granular materials. Granular materials are ubiquitous around us. They have many properties that are greatly different from those associated with common solids, liquids and gases. These unique properties lead to many applications. Many new ideas have emerged from this recent development. These systems have been studied from various different viewpoints experimentally, numerically, and theoretically. Nesterenkonoted that the propagation of a perturbation in loaded chain with Hertzian contacts possesses soliton-like features. The soliton-like behavior of the signal in two-dimensional granular beds was studied by Sinkovits and Sen. The soliton-like properties of the signal in the chain of iron balls were examined experimentally by Coste et al. Sen et al  studied the soliton-like pulse in perturbed and driven Hertzian chains and their possible applications in detecting buried impurities. For the vertical granular chain under gravity with power-law-type contact force, J. Hong et al  showed that there are two kinds of propagating modes, quasisolitary and oscillatory, depending on the strength of impulse. The type of dispersion and the oscillation frequency, wavelength, and period follow power laws in depth; J. Hong and Aiguo Xu found that the signal speed follows a power law in depth and the exponent varies smoothly from a saturant value to zero when the strength of the perturbation varies from very weak to very strong.

To explore the prospects for the application of impulse for detecting certain buried objects in granular beds, it is imperative that one must acquire a detailed understanding of the its properties. Up to now many problems on this topics are still openning. In present paper, we study the scattering process of impulse when it encounters some other kind of grains which we name impurities. The granular chain system can be described by the following motion equations

$`m_n\ddot{x}_n`$ $`=`$ $`V^{}(x_{n1},x_n)V^{}(x_n,x_{n+1})+m_ng,`$
$`n`$ $`=`$ $`1,2,\mathrm{},N.`$ (1)
where $`m_n`$ is the mass of the $`n`$th grain, $`V(x_n,x_{n+1})`$ is the inter-grain interaction between the $`n`$th and the$`(n+1)`$th grains. In this paper we study the impulses in a special case with $`V(x_{n1},x_n)=[\eta /(p+1)][\mathrm{\Delta }_0(x_nx_{n1})]^{p+1}`$, where $`\mathrm{\Delta }_0`$ is the diameter of the spherical grain, $`p`$ the exponent of the power-law type contact force, and $`\eta `$ the elastic constant of the grain under consideration. Thus equation (1) becomes

$`m_n\ddot{x}_n`$ $`=`$ $`\eta \{[\mathrm{\Delta }_0(x_nx_{n1})]^p[\mathrm{\Delta }_0(x_{n+1}x_n)]^p\}`$
$`+m_ng,n`$ $`=`$ $`1,2,\mathrm{},N.`$ (2)
For the Hertzian chain, i.e., the chain with Hertzian interaction between neighboring granular spheres, $`p=3/2`$.

## 2 Analytic treatment

In our system most of the grains are same and the second kind of grains only occur as impurities. Hence to analytically get some information of the system, we first fix the mass of all the grains to a constant, i.e., $`m_n=m`$ for all $`n`$. Then we introduce a new variable $`\psi _n`$, denoting the displacement of $`n`$th grain from equilibrium, defined by

$$\psi _n=x_nn\mathrm{\Delta }_0+\underset{l=1}{\overset{n}{}}\left(\frac{mgl}{\eta }\right)^{1/p},$$
(3)
where the last term in the right-hand side is the sum of grain overlaps up to $`n`$th contact and we set $`x_0=\psi _0=0`$. Eq. (2) is transformed into

$`m\ddot{\psi }_n`$ $`=`$ $`\eta \left[\left({\displaystyle \frac{mgn}{\eta }}\right)^{1/p}+(\psi _{n1}\psi _n)\right]^p`$ (4)
$``$ $`\eta \left[\left\{{\displaystyle \frac{mg(n+1)}{\eta }}\right\}^{1/p}+(\psi _n\psi _{n+1})\right]^p+mg`$
using Eq. (3).

For the weakly nonlinear regime, the condition

$$|\psi _{n1}\psi _n|\left(\frac{mgn}{\eta }\right)^{1/p}$$
(5)
is valid and the expansion of Eq. (4) under this condition reads

$$m\ddot{\psi }_n=\mu _n(\psi _n\psi _{n1})+\mu _{n+1}(\psi _{n+1}\psi _n)$$
(6)
where $`\mu _n=mpg\left(\frac{\eta }{mg}\right)^{1/p}n^{1\frac{1}{p}}`$ is the force constant of $`n`$th contact of the linear horizontal chain. We drop weak nonlinear terms of the expansion in Eq. (6) to make the system linear. Thus the vertical granular chain becomes a horizontal chain with varying force constants in which the gravity effect is contained. Both left and right side of Eq. (6) are linear in $`\psi _n`$. Therefore, the scaling analysis tells us that the equation of motion (6) has nothing to do with the initial impulse $`v_i`$. For this case we have obtained the following power-law,

$$v(h)h^{\frac{1}{4}\left(\frac{1}{3}+\frac{1}{p}\right)}\text{.}$$
(7)
The reducing of signal amplitude mainly due to the dispersion effect.

In our present studies two distances are important. The first one is the depth of the impurities $`d_{imp}`$. The second one is the distance between the impurity center and the grain where we get datas, $`d`$. Our main objectives are to get some knowledge on the reflection rate and transmission rate of the kinetic energy from the impurities. For fixed $`d`$, we expect the reflection rate follows a power-law and the exponent is independent of $`v_i`$.

For the strongly nonlinear regime, the condition

$$(\psi _{n1}\psi _n)\left(\frac{mgn}{\eta }\right)^{1/p}$$
(8)
is valid and the expansion under this condition leads Eq. (4) to

$`m\ddot{\psi }_n`$ $`=`$ $`\eta \left[(\psi _{n1}\psi _n)^p+pg_n(\psi _{n1}\psi _n)^{p1}\right]`$ (9)
$``$ $`\eta \left[(\psi _n\psi _{n+1})^p+pg_{n+1}(\psi _n\psi _{n+1})^{p1}\right]`$
where $`g_n=\left(mgn/\eta \right)^{1/p}`$ denotes grain overlap at $`n`$th contact. The gravity term can be neglected in the highly nonlinear regime, since the gravity effect appears in the coefficient $`g_n`$. Different order of $`\psi _n`$ in the left and right side of Eq. (9) implies that $`v_i`$-dependence must be appeared in the signal characteristics. Under the limitation $`v_i\mathrm{}`$, Eq. (9) recovers to

$$m\ddot{\psi }_n=\eta [(\psi _{n1}\psi _n)^p(\psi _n\psi _{n+1})^p]$$
(10)
That is to say, the vertical chain shows same behavior as that of the horizontal chain. For the horizontal chain the signal speed is a constant, so the exponent of the depth-dependence is $`0`$.

## 3 Simulation results

Our results are based on careful numerical integration of the coupled equations of motion for granular chains with $`N`$ spheres, where $`N`$ varies from $`1000`$ to $`5000`$ according to our need. The predictor-corrector method, Runge-Kutta method, Richardson extrapolation and the Bulirsch-Stoer method are used. The program units are just the same as those in Refs. and . That is to say, the units of distance, mass, and time are $`10^5m`$, $`2.36\times 10^5kg`$, and $`1.0102\times 10^3s`$, respectively. This kind of program units gives the gravitational acceleration $`g=1`$. We set the grain diameter $`100`$ and $`\eta =5657`$.

Figure 1 shows four snapshots for a scattering process, where the impurity depth is $`d_{imp}=500`$, which is described by the grain number; there are five impurities around the depth $`d_{imp}=500`$; the mass of the impurity grain is $`m_{imp}=0.1`$; the initial perturbation strength is $`v_i=0.1`$. Fig. 1 (a) shows a snapshot before the scattering process. Figs. 1 (b) and 1(c) show the scattering process. Fig. 1 (d) shows a snapshot where the scattering process nearly finished. We use “1st” and “2nd” to denote the leading and second leading peaks. The leading and second leading peaks of the incident impulse are pointed by an arrow, respectively, in Figs. 1 (a) and 1(b). The leading and second leading peaks of the scattered impulse are pointed by an arrow, respectively, in Figs. 1 (c) and 1 (d). It is clear that the leading peak of the scattering impulse is upwards, which is a characteristic property of the scattering process from light impurities. We use five impurities because that the scattered impulse is too weak if we use only one impurity. To understand the scattering process, we are interested in the variation of the ratios $`v_{ref}^{1st}/v_{inc}^{1st}`$ , $`v_{ref}^{2nd}/v_{inc}^{2nd}`$ and $`v_{ref}^{2nd}/v_{inc}^{1st}`$with the depth of the impurities $`d_{imp}`$, where $`v_{inc}^{1st}`$ and $`v_{ref}^{1st}`$ are maximum values of the leading peaks of the incident and reflected impulses, $`v_{inc}^{2nd}`$ and $`v_{ref}^{2nd}`$ are maximum values of the second leading peaks of the incident and reflected impulses. We obtain these datas from a gain with depth $`n=d_{imp}d`$. The three ratios describe the reflection rate from different sides. During the scattering process the shape of the impulse is not stable, so an appropriate distance $`d`$ is necessary. These ratios describe the reflection rate of the kinetic energy in the vertical granular chain. In order to get exact value for the leading and second leading peaks, we use the smooth velocity-time curve of the measured gain, instead of the rough snapshots.

Fig.2 shows the variations of $`v_{ref}^{1st}/v_{inc}^{1st}`$ , $`v_{ref}^{2nd}/v_{inc}^{2nd}`$ and $`v_{ref}^{2nd}/v_{inc}^{1st}`$ with the impurity depth $`d_{imp}`$. Fig. 2 (a) is for $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ and $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$, where $`d=200`$, $`v_i=1`$. The squares are numerical results for $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ and the circles are numerical results for $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$. The lines are corresponding fitting curves. Fig.2 (b) is for $`v_{ref}^{2nd}/v_{inc}^{1st}(d_{imp})`$, where $`d=200`$, $`v_i=5`$. The solid squares are for numerical results and the line is the fitting curve. We use the simplest power-law model, $`v_{ref}/v_{inc}d_{imp}^\beta `$ , to fit the numerical results. The fitting results suggest that the three ratios follow power-laws when $`d_{imp}`$ is large enough. From Refs.  we know that when the initial perturbation is weak the propagating mode is oscillatory, when the initial perturbation is strong the propagating mode is quasisolitary. It is interesting to check if the power-law model works for both the weak and strong perturbations. If it works, we have the following expected result: When $`v_i`$ is large enough, the behavior of the vertical chain is similar to that of the horizontal chain; hence there is no power-law. That means $`\beta 0`$ as $`v_i\mathrm{}`$.

Fig.3 shows the variation of the exponent for $`v_{ref}^{2nd}/v_{inc}^{2nd}`$ with the strength of initial perturbation, where the line with squares is for the case of $`d=100`$, the line with circles is for the case of $`d=200`$. The two curves show similar behaviors. It is clear that the exponent is independent of $`v_i`$ when $`v_i`$ is very small. For the case of $`d=100`$, $`\beta =0.19\pm 0.01`$. For the case of $`d=200`$, $`\beta =0.23\pm 0.01`$. In fact when $`v_i`$ is too small we can hardly get helpful information from the simulation datas due to numerical errors. The curve for $`d=100`$ has a peak around $`\stackrel{~}{v}_i=5`$. The curve for $`d=200`$ has a peak around $`\stackrel{~}{v}_i=10`$. When $`v_i>\stackrel{~}{v}`$ , $`\beta `$ decreases quickly with the increasing of $`v_i`$. We tried to add more points to the figure, but when $`v_i`$ is larger than $`20`$ we can hardly measure the second leading peak from the incident impulse. We can know the reason from Fig. 4.

Fig. 4 shows a $`v_n(t)`$ relation, where $`d_{imp}=300`$, $`d=200`$, $`n=100`$, $`v_i=15`$, and $`m_{imp}=0.1`$. The impulse near $`t=2`$ is the incident one. The impulse within $`9<t<10`$ is the reflected one (Only part of the reflected impulse is plotted). It is clear that for a strong perturbation the propagating impulse is quasisolitary and it has no evident the second leading peak; at the same time, the leading peak of the reflected impulse is very small and there exist some oscillations before the leading peak occurs.

When the perturbation is weak, the propagation of the wave is oscillatory; hence the two ratios $`v_{ref}^{1st}/v_{inc}^{1st}`$ and $`v_{ref}^{2nd}/v_{inc}^{2nd}`$ can describe well the reflection rate of the kinetic energy from the impurities. But with the increasing of $`v_i`$, the propagation of the wave is more and more quasisolitary and the ratio $`v_{ref}^{2nd}/v_{inc}^{1st}`$ describes the reflection rate better. Fig.5 shows the variation of the exponent $`\beta ^{}`$ for $`v_{ref}^{2nd}/v_{inc}^{1st}(d_{imp})`$ with $`v_i`$, where $`d=200`$. It shows similar behavior to the exponent $`\beta `$ for $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$. When $`v_i`$ is very small, $`\beta ^{}0.2`$.

Although we have found the power-law behavior for $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ for not very large $`v_i`$ (see Fig. 2 for an example), it is hard for us to get a figure which is similar to Fig.3 and Fig.5. The reasons will be given in the following part.

For a fixed distance $`d`$, when the depth of impurities $`d_{imp}`$ is large, we find the curve $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ oscillates greatly with $`d_{imp}`$. Fig. 6 shows two examples, where $`d=200`$; $`v_i=6`$ in (a) and $`v_i=4`$ in (b). We can understand the reason for this kind of behaviors through Fig. 7. Fig.7 shows several $`v_n(t)`$ curves from which we can study the scattering process of the impulse, where $`n=d_{imp}100`$, $`v_i=0.1`$, $`m_{imp}=0.1`$; $`d_{imp}=600`$ in (a), $`d_{imp}=800`$ in (b), $`d_{imp}=1000`$ in (c), and $`d_{imp}=1100`$ in (d), respectively. For the reflected impulse, only the leading peak and the second peak are shown. For a fixed value of $`d`$, when the depth $`d_{imp}`$ is small, the reflected impulse and the incident impulse are separated, i.e. there is a time during which the measured grain is static. For these cases, we can get exact values for $`v_{inc}^{1st}`$, $`v_{ref}^{1st}`$, $`v_{inc}^{2nd}`$, and $`v_{ref}^{2nd}`$. When the depth $`d_{imp}`$ is large enough, the tail of the incident impulse and the reflected impulse are connected. For these cases, we can hardly get exact values for $`v_{ref}^{1st}`$ and $`v_{ref}^{2nd}`$. From Fig.7 we also understand the fact these phenomena first influence the measurement of $`v_{ref}^{1st}`$ , then $`v_{ref}^{2nd}`$.

We name the increasing process of $`d_{imp}`$ before a power-law occurs a transient process. The simulation results show that the ratio $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ needs a much longer transient process than $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$. The initial impulse is stronger, we need a longer transient process. The curves for $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$ decreases with the increasing of $`d_{imp}`$, while the curves for $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ first increases, then decreases. The initial impulse is stronger, the curves decrease more slowly. We also notice that there are more fitting errors when $`v_i`$ is large.

## 4 Conclusion

In present paper we study the scattering process of impulse in vertical granular chain with light impurities. Most of the studies resort on numerical simulations and the analytic treatment helps to understand the behaviors. When the perturbation is weak, the quantities describing the reflection rate exhibit power law behavior with the impuity depth. The exponent is independent of $`v_i`$. When the perturbation is very strong, the vertical chain shows similar behavior to that of the horizontal chain, so the exponent is zero. Our numerical investigation begins from the weak perturbation region and extend to the nonlinear region and found a peak of the exponent. The difficulty in extending the numerical investigation to a stronger perturbation region is analyzed. More work needs to be done to completely understand the scattering process of the signal in vertical granular chain.

Acknowledgments

The author acknowledges Prof. Jongbae Hong for helpful suggestions and discussions. This work was supported by Korea Research Foundation Grant (KRF-2000-DP0106). The author also wishes to acknowledge the support by the BK21 project of the Ministry of Education, Korea.

Figure caption:

Figure 1: Snapshots for a scattering process, where the impurity depth is $`d_{imp}=500`$, which is described by the grain number; there are five impurities around the depth $`d_{imp}=500`$; the mass of the impurity grain is $`m_{imp}=0.1`$; the initial perturbation strength is $`v_i=0.1`$.

Figure 2: Variations of $`v_{ref}^{1st}/v_{inc}^{1st}`$ , $`v_{ref}^{2nd}/v_{inc}^{2nd}`$ and $`v_{ref}^{2nd}/v_{inc}^{1st}`$ with the depth of the impurities $`d_{imp}`$. Fig.2 (a) is for $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ and $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$ where $`d=200`$, $`v_i=1`$. The squares are numerical results for $`v_{ref}^{1st}/v_{inc}^{1st}(d_{imp})`$ and the circles are numerical results for $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$. The lines are corresponding fitting curves. Fig.2 (b) is for $`v_{ref}^{2nd}/v_{inc}^{1st}(d_{imp})`$ where $`d=200`$, $`v_i=5`$. The solid squares are for numerical results and the line is the fitting curve.

Figure 3: Variation of the exponent for $`v_{ref}^{2nd}/v_{inc}^{2nd}`$ with the strength of initial perturbation, where the line with squares is for the case of $`d=100`$, the line with circle is for the case of $`d=200`$.

Figure 4: The $`v_n(t)`$ relation, where $`d_{imp}=300`$, $`d=200`$, $`n=100`$, $`v_i=15`$, and $`m_{imp}=0.1`$. The impulse near $`t=2`$ is the incident one. The impulse within $`9<t<10`$ is the reflected one (Only part of the reflected impulse is plotted). The second leading peak of the incident impulse and the leading peak of the reflected impulse are very small. Some oscillations occur behind the main incident impulse.

Figure 5: Variation of the exponent $`\beta ^{}`$ for $`v_{ref}^{2nd}/v_{inc}^{1st}(d_{imp})`$ with $`v_i`$, where $`d=200`$. It shows similar behavior to the exponent $`\beta `$ for $`v_{ref}^{2nd}/v_{inc}^{2nd}(d_{imp})`$. When $`v_i`$ is very small, $`\beta ^{}0.2`$.

Figure 6: Reflection rate $`v_{ref}^{1st}/v_{inc}^{1st}`$ versus $`d_{imp}`$, where $`d=200`$; $`v_i=6`$ in (a) and $`v_i=4`$ in (b).

Figure 7: Velocity $`v_n(t)`$ curves, where $`n=d_{imp}100`$, $`v_i=0.1`$, $`m_{imp}=0.1`$; $`d_{imp}=600`$ in (a), $`d_{imp}=800`$ in (b), $`d_{imp}=1000`$ in (c), and $`d_{imp}=1100`$ in (d), respectively.


text:
# Hardy’s Second Axiom is insufficiently general

## Abstract

Hardy (quant-ph/0101012) conjectures in his Axiom 2 that $`K=K(N)`$, and that in classical probability $`K=N`$, while in quantum mechanics $`K=N^2`$. We offer an example in classical probability for which $`K=NV`$, $`V`$ the number of independent complete variables; with $`N=V`$ this classical example satisfies the purported quantal relation $`K=N^2`$.

Hardy’s Axiom 2 statesHardy01  that the informational degrees of freedom of a system, $`K`$, is a function of the number of values a variable is allowed, $`N`$: $`K=K(N)`$. He analyzes a classical probability system (balls in urns) to show that, in classical probability, $`K=N`$; in quantum mechanics, $`K=N^2`$. Hardy suggests that this difference characteristically distinguishes classical from quantal.

This is contradicted by the following simple example from classical probability (of a type I’ve discussed in Refs. Kirkpatrick:Quantal,  and Kirkpatrick:ThreeBox, ):

> The system has $`V`$ discrete-valued variables, call them Face, Suit, …, each taking on $`N`$ values, respectively $`\{\text{K},\text{Q},\mathrm{}\}`$, $`\{\text{S},\text{H},\mathrm{}\}`$, …. It consists of a box with a single display and pushbutton (momentary) switches marked “Face,” “Suit,” …. Inside the box is a deck of playing cards, each marked with a Face–value, a Suit–value, …(the cards are not necessarily unique). Some of the cards are segregated into a subdeck.
>
> When the switch labeled $`𝒱`$ (i.e., Face, Suit, …) is pressed, a card is selected from the subdeck by a nondeterministic process, and the value $`y`$ of the card’s variable $`𝒱`$ is shown momentarily in the display. A new subdeck is constructed, consisting of all the cards of the deck for which the value of the variable $`𝒱`$ is $`y`$.

In this system occurrences (“observations”) are repeatable: If, for a given occurrence, $`\text{Suit}=\text{H}`$ say, then a succeeding observation of Suit will return $`\text{Suit}=\text{H}`$ with certainty, while a succeeding observation of Face will randomly return K or Q or …. (These variables are incompatible: they do not have a joint probability distribution, nor are they simultaneously observable.)

It is necessary to determine the probability of each value of each variable in order to fully determine the state; thus, in this example, $`K=NV`$. Of course if $`V=1`$ (as is the case for his classical example, a ball in boxes), this corresponds to Hardy’s “classical” case, $`K=N`$. The maximal number of independent variables is equal to the dimension of the system: $`V=N`$, in which case $`K=N^2`$. This is the relation Hardy obtained for the quantum case—but this example is certainly not quantal.

In a quantum mechanical system of $`N`$ dimensions, there are maximally $`V=N`$ linearly independent variable operators, so $`K=N^2`$. However, if superselection rules are active, then there are fewer physical variables, $`V<N`$. Similarly, in the classical example presented here, the number of variables needn’t be maximal, so again it is possible that $`V<N`$.

Both in classical systems of the form exemplified here and in quantal systems we have $`K=VN`$, $`1<V<N`$. The functional form of $`K`$ cannot distinguish classical probabilistic systems from quantal systems.


text:
# Pattern formation inside bacteria: fluctuations due to low copy number of proteins

## Abstract

We examine fluctuation effects due to the low copy number of proteins involved in pattern–forming dynamics within a bacterium. We focus on a stochastic model of the oscillating MinCDE protein system regulating accurate cell division in E. coli. We find that, for some parameter regions, the protein concentrations are low enough that fluctuations are essential for the generation of patterns. We also examine the role of fluctuations in constraining protein concentration levels.

In recent years, dramatic experimental progress has been made in resolving the subcellular localization of bacterial proteins. Often the proteins form self–organized, spatially inhomogeneous patterns that can involve coherent spatiotemporal oscillations. These self-organized patterns are vital for accurate cell division Raskin ; Raskin1 ; Rothfield ; Hu ; Rowland ; Fu ; Hale . Understanding these patterns promises to reveal new mechanisms for the generation of subcellular bacterial structure. In this letter, we examine, for the first time, the impact of fluctuations on these patterns, focusing on the oscillatory MinCDE system in E. coli Raskin ; Raskin1 ; Hu ; Rowland ; Fu ; Rothfield ; Hale .

Each E. coli cell divides roughly every hour, first replicating its DNA into two nucleoids, and then dividing at midcell into two daughter cells. If division is not targeted accurately to midcell then DNA will not be distributed to both daughter cells, resulting in unviable anucleate “minicells”. Division is initiated by a polymeric ring of the protein FtsZ, which forms on the inner side of the cytoplasmic membrane. Precise positioning of the FtsZ ring to midcell is controlled both by the inhibiting effect of the nucleoids (“nucleoid occlusion”), and by the MinCDE system of proteins Yu . MinC inhibits the formation of the FtsZ ring, and is recruited to the membrane by MinD. MinE is also recruited to the membrane by MinD, where it forms a dynamical ring structure that drives pole to pole oscillations of MinC and MinD, with a period of about a minute. The oscillations lead to a time–averaged midcell MinC concentration minimum Hale . As MinC inhibits FtsZ ring formation, the FtsZ can only assemble near the cell midplane. The self–organized protein patterns are thus used by the cell to obtain positional information without stationary positional markers Howard ; MeinhardtBoer ; Kruse .

Very recently, the physical principles behind the MinCDE protein patterns have been explored through reaction–diffusion equations Howard ; MeinhardtBoer ; Kruse . The resulting equations represent protein diffusion, both in the cytoplasm and along the membrane, and also protein binding/unbinding from the cytoplasmic membrane. The slow membrane diffusion used in these models, relative to the cytoplasmic diffusion, results in Turing–like (Hopf) instabilities that spontaneously generate oscillatory patterns, in good agreement with experiment Howard ; MeinhardtBoer ; Kruse .

Bacterial proteins are, however, typically present in low numbers within the cell. This induces large fluctuations Mcadams ; Bialek  which have not been considered in previous pattern–forming models Howard ; MeinhardtBoer ; Kruse . In E. coli, a recent assay put the copy numbers for MinD and MinE at 2000 and 1400 respectively Shih . We have therefore investigated the role of fluctuations in a discrete particle model of the E. coli MinCDE system where each protein molecule is explicitly tracked. This allows for a full analysis of the fluctuations of both reactions and diffusion. Although the effects of noise have been studied in subcellular models without spatial dynamics Paulsson ; Vilar , and also in some spatially extended patterns Kessler ; Falcke ; Vilar1 ; Sancho ; Z , the effects of fluctuations on subcellular positional information, and the constraints on protein concentrations due to fluctuations, have not previously been studied. We find that:

(i) Even at surprisingly low concentrations, the noise does not destroy the oscillatory dynamics and indeed can be vital for generating patterns in regions of parameter space where the equivalent deterministic dynamics decays away. A bacterium can thus exploit low copy number fluctuations to produce stable, self–organized patterns. This result likely applies to any stochastic reaction–diffusion model of pattern formation involving sufficiently few protein copies (i.e. thousands or fewer).

(ii) We find evidence that the cell employs sufficient copy numbers of the MinCDE proteins to ensure reliable midcell division; using substantially fewer copies degrades accurate positioning of the midcell MinCD minimum, using more does not lead to significantly improved accuracy.

Stochastic model. We begin by introducing the stochastic model for the MinCDE dynamics, based on our deterministic model of Ref. Howard . We employ a $`1d`$ discrete particle model, where the particles hop between lattice sites and where the full fluctuation effects are intrinsically included by discrete particles. As in earlier models Howard ; MeinhardtBoer ; Kruse , the MinC dynamics is omitted, since it is known to closely follow the MinD dynamics. The occupancy at site $`i`$ is $`n_j^{\{i\}}`$, with $`j=\{D,d,E,e\}`$ representing cytoplasmic MinD, membrane MinD, cytoplasmic MinE, and membrane MinE, respectively. Each protein molecule is represented as a particle which, at each timestep, may hop with equal probability $`\stackrel{~}{D}_j\mathrm{\Delta }t/(\mathrm{\Delta }x)^2`$, where $`\mathrm{\Delta }t`$ is the time increment, to one of its neighboring sites at $`xx\pm \mathrm{\Delta }x`$ (except for the boundary sites at either end, where hard wall boundary conditions are imposed). At site $`i`$, the following reactions may occur, the first being for each $`D`$ particle, and then for each $`d,E,e`$ particle respectively:

$`\mathrm{Probability}:`$
$`n_D^{\{i\}}n_D^{\{i\}}1,n_d^{\{i\}}n_d^{\{i\}}+1\stackrel{~}{\sigma }_1\mathrm{\Delta }t/(1+\stackrel{~}{\sigma }_1^{}n_e^{\{i\}}),`$
$`n_D^{\{i\}}n_D^{\{i\}}+1,n_d^{\{i\}}n_d^{\{i\}}1\stackrel{~}{\sigma }_2\mathrm{\Delta }tn_e^{\{i\}},`$
$`n_E^{\{i\}}n_E^{\{i\}}1,n_e^{\{i\}}n_e^{\{i\}}+1\stackrel{~}{\sigma }_3\mathrm{\Delta }tn_D^{\{i\}},`$
$`n_E^{\{i\}}n_E^{\{i\}}+1,n_e^{\{i\}}n_e^{\{i\}}1\stackrel{~}{\sigma }_4\mathrm{\Delta }t/(1+\stackrel{~}{\sigma }_4^{}n_D^{\{i\}}).`$
These reactions are the stochastic analogs of the reaction processes used in our deterministic partial differential equation model Howard . The $`\stackrel{~}{\sigma }_1`$ term describes spontaneous membrane association of MinD, whereas the $`\stackrel{~}{\sigma }_2`$ term describes ejection of MinD from the membrane by membrane–bound MinE. Similarly, the $`\stackrel{~}{\sigma }_4`$ term describes spontaneous membrane disassociation of MinE, whereas the $`\stackrel{~}{\sigma }_3`$ term describes recruitment of MinE to the membrane by cytoplasmic MinD. The $`\stackrel{~}{\sigma }_1^{}`$, $`\stackrel{~}{\sigma }_4^{}`$ “suppression” terms correspond to membrane MinE suppressing the binding of MinD to the membrane, and to cytoplasmic MinD suppressing the unbinding of membrane MinE.

ATP dynamics. We next address the question of how ATP dynamics fits into the model. Experimentally, MinD is an ATPase deBoer  and it is the MinD-ATP complex that binds to the membrane, whereas the release of MinD back into the cytoplasm requires MinE–induced ATP hydrolysis Hu02 . Our model assumes that, following ATP hydrolysis and release of MinD into the cytoplasm, nucleotide exchange is sufficiently rapid to allow for membrane reattachment of MinD-ATP almost immediately. As a result, we only model MinD-ATP in the cytoplasm ($`n_D^{\{i\}}`$) and on the membrane ($`n_d^{\{i\}}`$). Like “actin treadmilling” in eukaryotic cells, the ATP driven binding and unbinding of MinD allows for a cyclic MinCDE pattern to be maintained with only low levels of protein synthesis: this makes the pattern–forming dynamics extremely energy efficient. This recycling is also consistent with experiments where protein synthesis was blocked, but where the MinCDE oscillations were observed to continue unaffected Raskin . Consequently, the above model (similar to that of Ref. Kruse ) does not include protein synthesis or degradation. These processes occur, but on longer time scales than the relatively rapid MinCDE oscillations.

Simulations. In our stochastic model simulations, we use time and spatial increments $`\mathrm{\Delta }x=0.02\mu m`$, $`\mathrm{\Delta }t=2\times 10^5s`$, so that $`100`$ lattice sites model a $`2\mu m`$ bacterium. We use: $`\stackrel{~}{D}_D=0.28\mu m^2s^1`$, $`\stackrel{~}{D}_d=0.003\mu m^2s^1`$, $`\stackrel{~}{D}_E=0.6\mu m^2s^1`$, $`\stackrel{~}{D}_e=0.006\mu m^2s^1`$, $`\stackrel{~}{\sigma }_1=20s^1`$, and $`\stackrel{~}{\sigma }_4=0.8s^1`$. Note that the membrane diffusion constants are much smaller than those in the cytoplasm; this agrees with recent data indicating that MinD may polymerize on the membrane Hu02 ; polymer . For the remaining variables of the model, we focus on 4 representative parameter sets shown in Table I, where we define $`N`$ as the total number of MinD proteins, equal to the total number of MinE proteins param . However, we emphasize that our results for the oscillatory behavior observed below are typical for large regions of parameter space. Initially, MinD and MinE particles are randomly distributed on the membrane and in the cytoplasm. Equal numbers of proteins are used since “wild–type” oscillations are observed when both proteins are equally expressed on plasmids (this is consistent, within experimental uncertainties, with the earlier quoted MinDE assay Shih ).

Fluctuation driven instability. Using the parameters in Table I, we find that the presence of noise is vital for the oscillations to persist. This is shown in Fig. 1, where the ratio of the average MinD density in the right–hand $`30\%`$ of the cell to that in the left–hand $`30\%`$ of the cell is plotted as a function of time for the stochastic model (at $`N=200`$ and $`N=1500`$) and for the deterministic model Howard  with equivalent parameters (i.e. with the above reaction probabilities directly transformed into deterministic reaction rates). In both cases, for the deterministic model, the protein concentrations rapidly decay away to the homogeneous steady state (in agreement with linear stability analysis Howard ), whereas regular oscillations continue for the stochastic model. Hence the average behavior of the stochastic model is clearly not describable using the naively equivalent deterministic model. To investigate this issue in more detail, we have examined how steady–state/oscillation bifurcations in the deterministic model are altered in the stochastic model. As a representative example, at $`N=1500`$, using the above parameter set in the stochastic model, but with $`\stackrel{~}{D}_d=0.001\mu m^2s^1`$, $`\stackrel{~}{D}_e=0.003\mu m^2s^1`$, $`\stackrel{~}{\sigma }_1^{}=0.2`$ and varying $`\stackrel{~}{\sigma }_4`$, we find that the transition from oscillatory to steady–state behavior is reduced by around $`40\%`$ from $`\stackrel{~}{\sigma }_4=0.63s^1`$ (deterministic) to $`\stackrel{~}{\sigma }_40.39s^1`$ (stochastic), as determined using the end to end MinD ratio. The noise does smear out the transition somewhat in the stochastic model (the transition at $`\stackrel{~}{\sigma }_40.39s^1`$ has width $`\pm 0.05s^1`$), but this effect is rather small. Hence this smearing out cannot account for the large regions of parameter space where oscillations occur in the stochastic, but not the deterministic, model. Rather, we have a fluctuation driven instability, where the noise has shifted the location of the transition, thereby promoting oscillations in large regions of parameter space where it would be forbidden in the equivalent deterministic model (see Fig. 1). Cells can in principle exploit low copy number fluctuations to generate pattern–forming dynamics. The oscillations continue down to very low concentrations ($`N=200`$), underlining the robustness of the dynamics.

The above examples show that the fluctuations are often essential for pattern formation. However, it also possible for the stochastic and deterministic models with equivalent parameters Howard  both to generate oscillations. Hence we cannot definitively conclude that fluctuations are essential for the MinCDE oscillations. Nevertheless, our analysis does show conclusively that cells can exploit low copy number fluctuations for the generation of dynamical subcellular structure.

Effect of fluctuations on the midcell MinD minimum. In Fig. 2, we plot the MinD and MinE concentration profiles for $`N=200`$ and $`N=1500`$, showing their averages over $`160`$ and $`110`$ successive cycles, respectively, and also for $`4`$ individual data sets each, where each set is averaged over only an individual oscillation cycle. For the long time data sets, we find that the midcell MinD concentration minimum (and a MinE concentration maximum) are still robustly reproduced even in the presence of noise. However, as can be seen from the data averaged over individual cycles, the fluctuations around this average can be very large for small $`N`$. In Fig. 3 we show histograms of the position of the MinD concentration minimum, where each minimum is determined over a single oscillation cycle (the use of a single cycle here is explained below). For $`N=1500`$, the histogram is sharply peaked around the cell center at $`1.0\pm 0.07\mu m`$ ($`1`$ standard deviation). As expected, with decreasing $`N`$ the width increases: $`0.09\mu m`$ at $`N=800`$, $`0.16\mu m`$ at $`N=400`$, and $`0.27\mu m`$ at $`N=200`$ param . The width of the midcell localization is large, particularly at protein counts ($`N=200`$ and $`N=400`$) that are significantly below those seen naturally. Hence, using significantly fewer protein copies degrades the accuracy of midcell division.

Nucleoids, when present, also affect the positioning of the FtsZ ring through the poorly understood phenomenon of “nucleoid occlusion” Yu , where FtsZ rings do not nucleate over nucleoids and are restricted to either near the midcell or at the cell poles. Segregated nucleoids (at the $`1/4`$ and $`3/4`$ positions along the cell) will truncate the tails of the distributions shown in Fig. 3, further enhancing the accuracy of midcell division (in agreement with experiment sun1998 ; Yu ). In normal cells with nucleoids, it is particularly important that the MinCDE system block polar FtsZ rings, since the nucleoids themselves will inhibit FtsZ rings elsewhere away from midcell. Assuming that FtsZ nucleation occurs at a single cycle MinD minimum (see below), then from Fig. 3 we see that $`N=1500`$ is a high enough concentration to reduce the probability of polar division to considerably less than $`0.01`$ per oscillation cycle. Given that about $`50`$ complete oscillation cycles normally occur between successive divisions, we see that attaining this level of accuracy is important. Significantly lower concentrations than $`N=1500`$ will lead to an unacceptable probability of polar division, while higher concentrations will lead to only marginally increased accuracy, but at the cost of manufacturing many additional protein copies. From these simple arguments based on fluctuation effects, we see that E. coli may be using an optimal number of Min proteins, trading off midpoint precision against the cost of protein synthesis. There will be other constraints on the protein copy numbers (e.g., sufficient MinC to successfully inhibit off–center FtsZ ring formation), but fluctuations set useful bounds on the concentration levels.

Comparison with experiment. Experimentally, the precision of the MinCDE system can be probed in anucleate cells by measuring the position of the FtsZ ring. In these cells the only positional guide for division is the MinCDE system, which functions even in the absence of the nucleoids. Indeed in these anucleate cells the FtsZ ring position is placed at midcell with a width of $`0.12\mu m`$ (scaled to a $`2\mu m`$ length) Yu , somewhat larger than the MinD distribution width $`0.07\mu m`$ at $`N=1500`$. We would expect that FtsZ ring nucleation would not precisely track the MinCD minimum, which could account for the difference. In order to make this comparison we need to know how many cycles the FtsZ nucleation averages over to identify the location of the MinD minimum. Experimentally, this timescale is on the order of a minute, since, from Ref. Raskin , oscillation cycles of between $`30s`$ and $`120s`$ give normal division. As referred to earlier, this justifies a rough comparison of the width of the single cycle ($`90s`$) MinD density minimum distribution at $`N=1500`$ with that of the FtsZ distribution in anucleate cells.

In conclusion, we have studied fluctuations in the E. coli MinCDE system. In some regimes, we have found that fluctuations are essential for pattern generation. We have also found evidence that the MinCDE concentrations may be optimal for reliable midcell division. Based on the MinCDE system, we see that $`O(10^3)`$ copies of pattern forming proteins are required in bacteria to obtain positional information accurate to within a few percent.

We acknowledge useful conversations with P. de Boer, M. van Hecke, K. Kruse, J. Lutkenhaus, and S. de Vet. M.H. acknowledges financial support from Stichting FOM and from The Royal Society. A.D.R acknowledges financial support from NSERC Canada.

Present address: Department of Mathematics, Imperial College London, South Kensington Campus, London SW7 2AZ, U.K.