text:
# Evolutionary models for low-mass stars and brown dwarfs: uncertainties and limits at very young ages

## 1 Introduction

The development of a new generation of stellar evolution models based on the accurate coupling between interior and atmosphere models yielded a major advance in the description of very-low-mass stars (VLMS; Baraffe et al. 1995, 1997, 1998, BCAH98) and substellar objects (brown dwarfs BD and giant planets GP; Burrows et al. 1997 for objects with $`T_{\mathrm{eff}}<2000`$K; Chabrier et al. 2000a, CBAH00). One of the main advantages of such models over previous generation models is the direct comparison between theory and observation in colour-colour and colour-magnitude diagrams (CMD). Several observational tests now assess the validity of the theory devoted to the description of low-mass ($`<\mathrm{\hspace{0.17em}1}M_{}`$) astrophysical objects (see Allard et al. 1997; Chabrier and Baraffe 2000, for recent reviews). Although some discrepancies between models and observations still remain, uncertainties due to the input physics are now significantly reduced.

In previous papers, we have developed evolutionary models of VLMS and BDs (BCAH98; CBAH00), taking into account the most recent improvements of the physics describing the interior (equation of state for dense plasmas, screening factors) and the atmosphere (molecular opacities, formation of dust) devoted to the analysis of objects with an age $`t>`$ 100 Myr. Comparison between observations and models for very young objects ($`t<`$ 100 Myr) are more uncertain, in particular from the observational viewpoint, mainly for two reasons : (1) extinction due to the surrounding dust modifies both the intrinsic magnitude and the colours of the objects and (2) spectra of very young objects ($`t<`$ 1 Myr) may still be affected by the presence of an accretion disk or circumstellar material residual from the protostellar stage. On the theoretical side, an important source of uncertainty at the very early stages of evolution is the choice of the initial conditions, resulting from prior protostellar collapse and accretion phases. The simple picture of non-accreting objects contracting from large initial radii, as used in our models, is clearly an idealized description of reality (Stahler, Shu & Taam 1980; Stahler 1983, 1988; Hartmann, Cassen & Kenyon 1997 and references therein).

In spite of these uncertainties, numerous surveys are devoted to very young clusters because of the potential detection of substellar objects down to planetary masses (Zapatero et al. 2000; Lucas et al. 2001), and the determination of the very (sub)stellar initial mass functions (IMF), free of dynamical evolution effects that affect older clusters. Consequently, a wealth of data for low mass objects with ages spanning $``$ 1-10 Myr is available and provides the basis for a better understanding of the early phases of stellar evolution. Given the reliability of the present theory for VLMS and BD, a confrontation of the models with the observational data of such very young objects provides important insight into the early epoch of protostellar collapse and accretion phases. For this purpose, however, a correct appreciation of the model uncertainties and limitations is required.

For this reason, this paper is specifically devoted to the analysis of evolutionary models by BCAH98 and CBAH00 at early phases of evolution ($`t<`$ 100 Myr). The input physics is briefly recalled in section 2, followed by a discussion of uncertainties specific to evolutionary tracks at early ages (§3). Our models are compared to other evolutionary tracks available in the literature in §4 and the confrontation to observations is presented in §5.

## 2 Model description

### 2.1 Evolutionary tracks

The models analysed in the present paper are based on the input physics already described in BCAH98 and CBAH00. Both sets of models use the same ingredients describing the stellar interior but use different sets of atmosphere models, which provide the outer boundary conditions and the synthetic spectra. The BCAH98 evolutionary tracks are based on the non-grey atmosphere models by Hauschildt, Allard and Baron (1999a). These models are dust-free and are appropriate to the description of objects with effective temperatures $`T_{\mathrm{eff}}>\mathrm{\hspace{0.17em}\hspace{0.17em}2300}`$ K. The CBAH00 models are based on atmospheres including the formation and opacity of dust (Allard et al. 2001, hereafter DUSTY models ). As illustrated in CBAH00, dust must be taken into account in order to explain the near-IR colors of late M-dwarfs and L-dwarfs. The latter models are thus more appropriate to the description of objects with $`T_{\mathrm{eff}}<\mathrm{\hspace{0.17em}\hspace{0.17em}2300}`$ K. As emphasized in CBAH00, the DUSTY models are not appropriate for the description of spectral and photometric properties of methane dwarfs ($`T_{\mathrm{eff}}<`$ 1600K), which require a different treatment of dust (the so-called COND models in CBAH00 and Allard et al. 2001). Evolutionary models for such cool objects will be presented elsewhere (Chabrier et al. 2001, in preparation).

The BCAH98 grid covered a mass range from 0.02 $`M_{}`$ to 1.4 $`M_{}`$, for ages $``$ 1 Myr up to the Main Sequence for stars. Originally, the CBAH00 grid covered masses from 0.01 $`M_{}`$ to 0.1 $`M_{}`$ for ages $``$ 10 Myr. In the present work, we extend it down to 1 $`M_\mathrm{J}`$ (10$`{}_{}{}^{3}M_{}^{}`$) and ages $``$ 1 Myr. Figure 1 presents the complete grid of models in a Hertzsprung-Russell diagram (HRD) from 0.001 $`M_{}`$ to 1.4 $`M_{}`$. The small variation of radius with mass and age in the substellar regime (see e.g. Chabrier & Baraffe, 2000) yields the merging towards very similar tracks below the hydrogen-burning limit. Figure 2 displays the time evolution of the effective temperature and luminosity for selected masses. Objects below 2 $`M_\mathrm{J}`$ evolve essentially with $`T_{\mathrm{eff}}<`$ 1600K (see Figs. 1 and 2), even at very early ages. Their atmospheric properties are thus better described by the COND models.

The initial conditions of the models are described in the next section. The stellar/substellar transition is located at $`m`$ 0.075 $`M_{}`$. Below this limit, objects become partially degenerate and their nuclear energy production cannot compensate the energy lost by radiation (e.g $`L_{\mathrm{nuc}}<L`$), which is required to reach the Main Sequence. The deuterium burning minimum mass is $`m`$ 0.012 $`M_{}`$ (Saumon et al. 1996; Chabrier et al. 2000b). The initial D-burning phase lasts less than 1 Myr for $`m>\mathrm{\hspace{0.17em}0.2}M_{}`$, between 1 and 5 Myr for $`0.05<m<\mathrm{\hspace{0.17em}0.2}M_{}`$ and almost 20 Myr for a 0.02 $`M_{}`$ brown dwarf (see Chabrier et al. 2000b, for details).

As shown in CBAH00, grains have little effects on the evolution of $`L(t)`$ and $`T_{\mathrm{eff}}(t)`$, because of the reduced dependence of evolution upon opacity, $`L(t)\kappa _R^{1/3}`$, $`T_{\mathrm{eff}}(t)\kappa _R^{1/10}`$ (Burrows & Liebert, 1993). We verified that the difference between the different TiO and H2O molecular line lists used in BCAH98 and CBAH00 models (see §3 in CBAH00), respectively, affect essentially the outer atmospheric layers, and thus the synthetic spectra and colors, but not the deeper atmospheric layers, and thus the outer boundary conditions. The effect of these different molecular linelists on the evolution of the effective temperature $`T_{\mathrm{eff}}(t)`$ and the bolometric luminosity $`L(t)`$ is small, less than 100 K in $`T_{\mathrm{eff}}`$ and 10% in $`L`$ at a given age. This is illustrated in Figure 2, where tracks based on the BCAH98 input physics (solid lines) are compared to the CBAH00 tracks (dashed lines) for the same masses (0.075, 0.02 and 0.012 $`M_{}`$). As stressed in CBAH00 and Baraffe et al. 2001, the computation of more reliable H2O and, to a lesser extend, TiO linelists is badly needed to solve this shortcoming in the present theory.

## 3 Main uncertainties: initial model and convection

### 3.1 Initial conditions

As mentioned in the previous section, although shortcomings still remain in current molecular opacities, the resulting uncertainty on the evolution is small. One of the main source of uncertainty for models at early stages of evolution is the choice of the initial conditions. Most of low mass pre-Main Sequence (PMS) models available in the literature (D’Antona & Mazzitelli 1994, 1997; Burrows et al. 1997; BCAH98; Siess et al. 2000) start from arbitrary initial conditions, totally independent on the outcome of the prior proto-stellar collapse and accretion phases. The initial configuration is that of a fully convective object starting its contraction along the Hayashi line from arbitrary large radii. Evolution starts prior to or at central deuterium ignition, with initial central temperature $`\mathrm{\hspace{0.17em}5\hspace{0.17em}10}^5`$ K. According to studies of low-mass protostellar collapse and accretion phases, such initial conditions are oversimplified, and low mass objects should rather form with relatively small radii (Hartmann et al. 1997, and references therein). Based on spherical accretion protostellar models, Stahler (1983, 1988) defined a birthline in the Hertzsprung-Russell diagram where young objects become visible. Evolutionary tracks should then start from this birthline, which fixes the age $`t`$ = 0. Ages determined from models based on the above-mentioned oversimplified conditions should then be corrected accordingly, with substantial corrections for systems younger than a few Myr (see Palla and Stahler 1999). Furthermore, collapse and accretion are unlikely to proceed spherically. Spherical collapse does not consider angular momentum transport, an important issue of the early phases, which affects the subsequent cooling and formation of the protostar. Hartmann et al. (1997) recently illustrate the sensitivity of the birthline locus assuming that accretion proceeds through a disk rather than spherically. This work stresses again the high uncertainty of assigning ages from HRD positions for the youngest objects.

This analysis demonstrates convincingly that assigning an age to objects younger than a few Myr is totally meaningless when the age is based on models using oversimplified initial conditions. As shown in the next section, for $`t<\mathrm{\hspace{0.17em}\hspace{0.17em}1}`$ Myr, the evolutionary tracks themselves are sensitive to the initial conditions, whereas after a few Myr, the models converge toward the same track. This is the main reason why we provide confidently evolutionary models for ages $`t\mathrm{\hspace{0.17em}1}`$ Myr, considering that below such ages, models are too uncertain. To solve this substantial uncertainty requires the consistent evolution between the 3D collapse of the protostellar phase and the subsequent PMS evolution.

#### 3.1.1 Effect of the initial radius and the mixing length parameter

In this section we examine the effect of the variation of the initial radius (or initial gravity) in the initial conditions. For this purpose we have extended our grid of atmosphere models to gravities $`\mathrm{log}g<\mathrm{\hspace{0.17em}3.5}`$. We have checked the validity of the plane-parallel approximation used in these models by comparison with similar atmosphere models which take into account the effects of spherical geometry. These effects are found to become important only for surface gravities $`\mathrm{log}g<\mathrm{\hspace{0.17em}\hspace{0.17em}2}`$ (Hauschildt et al. 1999b; Allard, Hauschildt & Schweitzer 2000), well below the range of gravities characteristic of the evolution of young low mass objects.

We first calculate a set of models, labeled (A), with initial radii fixed to obtain initial surface gravities $`\mathrm{log}g\mathrm{\hspace{0.17em}3}3.5`$ and initial thermal time-scales $`t_{th}`$ a few Myr. For masses above the deuterium burning minimum mass (DBMM), $`m>\mathrm{\hspace{0.17em}\hspace{0.17em}0.02}M_{}`$, such tracks start at deuterium ignition, with central temperature $`\mathrm{5\hspace{0.17em}10}^5\mathrm{\hspace{0.17em}10}^6`$K, defining a Zero-Age-Deuterium Burning Sequence. Such initial conditions are similar to that used in BCAH98 and CBAH00. A second set of models, labeled (B), starts with larger radii such that the initial surface gravity $`\mathrm{log}g\mathrm{\hspace{0.17em}2.5}`$. These initial models are more luminous than the previous ones, with central temperatures below 5 105K and initial thermal time-scales $`t_{th}\mathrm{\hspace{0.17em}10}^5`$ yr.

We also analyse the sensitivity of the models to the mixing length $`l_{\mathrm{mix}}`$, characteristic of the mixing length formalism (MLT) used to describe convection. We used two different values of the mixing length parameter $`\alpha _{\mathrm{mix}}`$111$`\alpha _{\mathrm{mix}}=l_{\mathrm{mix}}/H_\mathrm{P}`$, with $`H_\mathrm{P}`$ the pressure scale-height = 1 and 2 (see §3.2 for justifications).

The time evolution of $`L`$ and $`T_{\mathrm{eff}}`$ for models (A) (solid lines) and (B) (dashed lines) is displayed in Figure 3 for several masses and for $`\alpha _{\mathrm{mix}}`$ = 1. Figure 3 shows the importance of the initial radius on $`T_{\mathrm{eff}}`$, and thus on $`L`$, during the first Myr of evolution. After a few Myr, however, the tracks based on different initial radii merge, for a given mass. For masses ranging from $``$ 0.01 to 0.5 $`M_{}`$, the two sets of models start at very different $`T_{\mathrm{eff}}`$. Models starting with the lowest gravity are cooler by up to several hundreds K compared to initially denser, less luminous models with the same mass. Note also that for models (B), $`T_{\mathrm{eff}}`$ increases during the early evolution, under contraction, in contrast to the first set of models with initial $`\mathrm{log}g>\mathrm{\hspace{0.17em}\hspace{0.17em}3.0}`$. This is the consequence of the different surface gravities, which strongly affect the atmosphere profiles for $`T_{\mathrm{eff}}`$ = 2200-3500K, as illustrated in Figure 4. As shown in this figure, for this range of effective temperatures, the atmosphere profiles also show an extreme sensitivity to the mixing length $`l_{\mathrm{mix}}`$ at low gravity ($`\mathrm{log}g<\mathrm{\hspace{0.17em}3.5}`$) but become rather insensitive to a variation of this parameter for $`\mathrm{log}g\mathrm{\hspace{0.17em}3.5}`$ (see Fig. 4b).

Such a sensitivity of the atmospheric profiles to gravity and mixing length in this range of $`T_{\mathrm{eff}}`$ stems from the less efficient formation of molecular H2 as gravity decreases. Figure 5 displays the fraction of H2 along the atmosphere profiles shown in Fig. 4 (below the photosphere). For $`T_{\mathrm{eff}}`$ = 2500K and $`\mathrm{log}g=2.5`$, the fraction $`n_{\mathrm{H}_2}`$ of H2 decreases rapidly along the inner profile. For $`\alpha _{\mathrm{mix}}`$ = 1, it becomes negligible at optical depth $`\tau `$ = 100, where the boundary condition for the inner structure is defined (see Chabrier & Baraffe 1997, their §2.5). For $`\mathrm{log}g\mathrm{\hspace{0.17em}3.0}`$ (see Fig. 5), $`n_{\mathrm{H}_2}`$ exceeds 20%, even in the deepest atmospheric layers.

As described in Chabrier & Baraffe (2000, §2.2.2), hydrogen atom recombination yields a rapid increase of the H2 Collision Induced Absorption opacity ($`\kappa _{\mathrm{CIA}}n_{\mathrm{H}_2}^2`$) and a decrease of the adiabatic gradient . Both effects favor convection, yielding a flatter atmosphere profile, i.e a decrease of $`T`$ at fixed $`P`$. This in turn favors the formation of H2 and increases $`n_{\mathrm{H}_2}`$ at fixed $`P`$, increasing again the opacity and the efficiency of convection. This illustrates the non-linear response of the atmospheric structure to the formation of H2, as soon as its fraction becomes significant. It explains the huge effect of gravity when this latter increases from $`\mathrm{log}g=`$ 2.5 to 3.0 on the atmospheric profile, as illustrated in Figure 4a. The decrease of $`_{ad}`$ in regions of H2 formation favors adiabatic convection, limiting the extension of super-adiabatic layers and thus the sensitivity of the thermal profile to the mixing length parameter. This explains the strong dependence of the atmosphere profile on $`\alpha _{\mathrm{mix}}`$ for $`\mathrm{log}g`$ = 2.5 (see Fig. 4b). At higher gravities convection described by the MLT is essentially adiabatic and almost insensitive to the choice of $`\alpha _{\mathrm{mix}}`$ (see §3.2).

A similar behavior is found in the range of effective temperatures 2200 K $`<T_{\mathrm{eff}}<`$ 4000K. Above $``$ 4000K, the structure is too hot for H2 to form in significant fraction whereas below $`T_{\mathrm{eff}}\mathrm{\hspace{0.17em}2200}K`$, the outer layers are dense and cool enough for H2 to form efficiently, even at $`\mathrm{log}g=\mathrm{\hspace{0.17em}2.5}`$, .

The drastic modification of the atmospheric profile, for $`\alpha _{\mathrm{mix}}`$ = 1, as gravity increases from $`\mathrm{log}g`$ = 2.5 to $`\mathrm{log}g`$ 3.0 explains the different $`T_{\mathrm{eff}}`$’s for models starting from different initial radii (Fig. 3). For models (B), starting with $`\mathrm{log}g`$ = 2.5, gravity increases as contraction proceeds and favors H2 molecular formation, yielding significantly flatter atmospheric profiles, as mentioned above. This yields a significant increase of $`T_{\mathrm{eff}}`$ at fixed mass $`m`$ as contraction proceeds (Chabrier & Baraffe, 2000, §3.2). For $`\alpha _{\mathrm{mix}}`$ = 2, the sensitivity of atmosphere profiles to $`g`$ is less pronounced, but still yields up to 200 K differences between models (A) and (B) at a given age.

The effect on evolutionary tracks in a HRD is displayed in Figure 6. Note that evolution along the Hayashi line does not necessarily proceed at constant $`T_{\mathrm{eff}}`$, because of the effects described above. The common picture of vertical (constant $`T_{\mathrm{eff}}`$) Hayashi tracks is therefore an oversimplified picture of PMS evolution. As demonstrated in the pioneering papers of Hayashi (1961) and Hayashi & Nakano (1963), fully convective and adiabatic objects which have H2 dissociation zones contract almost vertically in the HRD. The very low value of $`_{ad}`$ ($`<0.1`$) in such zones is responsible for such evolution at constant $`T_{\mathrm{eff}}`$. Taking into account super-adiabaticity in convective layers displace the Hayashi line toward lower effective temperature. This is equivalent to a decrease of $`\alpha _{\mathrm{mix}}`$ in our calculations. Hayashi & Nakano (1963) have also shown that omitting the presence of H2 molecules yields an evolution proceeding at decreasing $`T_{\mathrm{eff}}`$ (from the left to the right in a HRD), rather than at constant $`T_{\mathrm{eff}}`$. Although based on a simplified treatment of the atmospheric properties, opacities and equation of state, these basic works already illustrated the extreme sensitivity of the shape of PMS tracks for low mass objects to super-adiabaticity and to the degree of H2 formation/dissociation.

Models (A) and (B) follow the same track for a given mass and $`\alpha _{\mathrm{mix}}`$, but do not reach the same position at the same age. Significant differences appear at ages $`<`$$``$ 1 Myr but vanish after a few Myr. We thus consider 1 Myr as the characteristic time required to forget our arbitrary initial conditions and below which models are too sensitive to input physics and thus too uncertain.

#### 3.1.2 Effect of the initial deuterium abundance

Deuterium burning plays a key role during the protostellar collapse phase (Stahler 1988) and the following $``$ 10 Myr of evolution. In a previous paper, we focused on the initial deuterium burning phase (Chabrier et al. 2000b), adopting an initial mass fraction \[$`{}_{}{}^{2}D_{0}^{}`$\] = 2 10-5, characteristic of the local interstellar medium (hereafter LISM, Linsky 1998). In the present section we examine the effect of a variation of \[$`{}_{}{}^{2}D_{0}^{}`$\] by a factor of 2 on the early stages of evolution. Such a variation is motivated by a recent deuterium abundance determination towards quasars, which suggests a primordial abundance only slightly larger (by less than a factor of 2) than the LISM value (O’Meara et al. 2001).

The main effect of a variation of the initial \[$`{}_{}{}^{2}D_{0}^{}`$\] is the modification of the age at a given $`L`$ and $`T_{\mathrm{eff}}`$, with the most important effect for masses $`m<\mathrm{\hspace{0.17em}0.07}M_{}`$ (Figure 7). For a standard \[$`{}_{}{}^{2}D_{0}^{}`$\], a 0.07 $`M_{}`$ depletes by a factor of 100 its initial deuterium content within $`3`$ Myr and a 0.02 $`M_{}`$ needs $`17`$ Myr (see Chabrier et al. 2000b). If \[$`{}_{}{}^{2}D_{0}^{}`$\] is increased by a factor of 2, the 0.07 $`M_{}`$ brown dwarf then needs $`5`$ Myr and the 0.02 $`M_{}`$ requires $`26`$ Myr for the same depletion factor. A 0.015 $`M_{}`$ brown dwarf requires 50 Myr in the standard case and 70 Myr if \[$`{}_{}{}^{2}D_{0}^{}`$\] is twice as large to reach the 99% depletion limit. Conversely, if \[$`{}_{}{}^{2}D_{0}^{}`$\] is smaller by a factor of 2, the 99% depletion limit is reached in 1.7, 12 and 40 Myr for respectively 0.07, 0.02 and 0.015 $`M_{}`$. If the initial \[$`{}_{}{}^{2}D_{0}^{}`$\] is as small as 2$`\times `$10-6, e.g a factor of 10 smaller than the standard LISM value, all objects with masses $`m>\mathrm{\hspace{0.17em}0.04}M_{}`$ depletes their initial deuterium within 2 Myr.

The effect of a variation of \[$`{}_{}{}^{2}D_{0}^{}`$\] on isochrones in CMDs, however, remains small, compared to other sources of uncertainties at young ages (extinction, $`T_{\mathrm{eff}}`$ calibration, initial models) and can be ignored for this present. If the improvement of observable techniques in the future, however, allows the detection of deuterium in the atmosphere of young objects (Chabrier et al. 2000b), this effect needs to be taken into account for a correct age estimate of very young clusters.

### 3.2 Convection

One of the major uncertainties affecting the evolution of stars with masses $`m>\mathrm{\hspace{0.17em}0.6}M_{}`$ is due to the treatment of convection. These stars show relatively extended super-adiabatic outer layers during PMS and MS evolution, of which description is extremely sensitive to the adopted treatment of convection, i.e., within the MLT formalism, to the adopted mixing length (at any gravity). The effect of a variation of $`l_{\mathrm{mix}}`$ on evolutionary tracks for solar-type stars is well known and is illustrated e.g in Figure 2b of Baraffe et al. (2001). Ludwig, Freytag & Steffen (1999) calibrated the mixing length parameter $`\alpha _{\mathrm{mix}}`$ for these stars with 2D hydrodynamical models performed in the parameter space 4300 K $`T_{\mathrm{eff}}`$ 7100 K and gravities 2.54 $`\mathrm{log}g`$ 4.74. They found a moderate variation of the mixing length parameter around typically $`\alpha =`$1.5. Note that for the present models and input physics 222As stressed in BCAH98, the value of $`l_{\mathrm{mix}}`$ required to fit the Sun depends on the input physics (equation of state, opacities, outer boundary condition), $`l_{\mathrm{mix}}`$ = 1.9 $`H_\mathrm{P}`$ is the value required to fit the Sun at its present age. An increase of $`\alpha _{\mathrm{mix}}`$ from 1 to 2 yields an increase of $`T_{\mathrm{eff}}`$ up to 500K for the highest masses during their PMS evolution.

For masses $`m<\mathrm{\hspace{0.17em}0.6}M_{}`$, the extension of the super-adiabatic layers retracts appreciably and the transition from convective to radiative outer layers is characterized by an abrupt transition from a fully adiabatic to a radiative structure with a very small entropy jump. This means that during most of the evolution, except at early ages, as discussed in §3.1.1 and below, the sensitivity of the evolutionary models to $`l_{\mathrm{mix}}`$ is small. Multi-dimensional hydrodynamical simulations for conditions characteristic of M-dwarf atmospheres, $`T_{\mathrm{eff}}`$ 3000 K, $`\mathrm{log}g=5`$, have been recently conducted by Ludwig, Allard & Hauschildt (2001). These simulations confirm the afore-mentioned small entropy jump found in the 1-D models described by the MLT, illustrating the large efficiency of atmospheric convection for these objects, a direct consequence of the formation of molecules, as mentioned earlier. Under these circumstances, the 3D simulations show that the MLT does indeed provide a correct thermal profile, providing a value of $`\alpha _{\mathrm{mix}}>`$ 1, at least for high gravities and relatively old objects ($`t>>`$ 10 Myr).

In §3.1.1, however, we have shown that even below 0.6 $`M_{}`$, very young models with gravities $`\mathrm{log}g<`$ 4 can be affected by a variation of $`l_{\mathrm{mix}}`$. To minimize such uncertainties, a correct calibration of $`l_{\mathrm{mix}}`$ requires the extension of the Ludwig et al. (2001) calculations to lower gravities. This work is under progress and represents the most promising method for an accurate description of convection in optically-thin media, through an accurate calibration of $`l_{\mathrm{mix}}`$. In contrast, given all the already mentioned uncertainties inherent to either observation or models for very young objects, a calibration of $`l_{\mathrm{mix}}`$ based on a comparison of PMS tracks with observations is at best highly speculative.

## 4 Comparison with other work

A detailed comparison between various evolutionary tracks available in the literature has already been done by Siess, Dufour & Forestini (2000). In the present section, we compare the evolutionary tracks most widely used by the community to describe the observational properties of objects in young clusters (Burrows et al. 1997, B97; D’Antona & Mazzitelli 1994, 1997, DM94, DM97; Palla & Stahler 1999, PS99). The comparison between these different models in a theoretical HR diagram is illustrated in Figure 8 for a few low-mass star and BD masses.

Let us first summarize the main differences between these models. One crucial difference is due to the outer boundary condition: DM94, DM97 and PS99 use approximate boundary conditions based on $`T(\tau )`$ relationships assuming gray approximation and radiative equilibrium. Outer boundary conditions based on such approximations are wrong as soon as molecules form near the photosphere and convection reaches optically thin layers, i.e below $`T_{\mathrm{eff}}4000K`$ (see Chabrier and Baraffe 1997, 2000; and references therein) They usually yield hotter effective temperatures for a given mass. Above $`4000K`$, the choice of the outer boundary condition is less consequential, as shown by the good agreement between the 0.8 $`M_{}`$ tracks of BCAH98 and of PS99 (see Figure 8), all the other physical inputs being similar for such mass. Note that, although the PS99 models start from a more realistic birthline, this does not affect the results after 1 Myr, as expected from the tests on initial conditions performed in §3.1.1.

The B97 models for low mass stars and hot brown dwarfs ($`T_{\mathrm{eff}}>\mathrm{\hspace{0.17em}2000}K`$) are based on gray atmosphere models obtained by solving the radiative transfer equation, as described in Burrows et al. (1993). Such an approximation, although it represents an improvement over the previous $`T(\tau )`$ relationships, still overestimates $`T_{\mathrm{eff}}`$ at a given mass compared to models based on full non-gray atmosphere models. This is certainly the reason why for VLMS on the Main Sequence (for $`m>0.08M_{}`$), the B97 models are about 100-200K hotter than the BCAH98 ones. We note, however, that the very early evolution of the B97 models proceeds at much hotter $`T_{\mathrm{eff}}`$ than any other model. Although also true for masses below 0.2 $`M_{}`$, this does not appear in Fig. 8 since the tracks start at 1 Myr and the 0.06 $`M_{}`$ of B97 displayed in Fig. 8 is indeed much hotter but for $`t<`$ 1 Myr. The reasons for such differences along the Hayashi line is not clear. We only notice that similar shape of the Hayashi line (i.e. a strong decrease of $`T_{\mathrm{eff}}`$) is obtained if one assumes a fully adiabatic initial structure and adiabatic convection (i.e $`\alpha _{\mathrm{mix}}\mathrm{}`$). Based on test cases, we find that such initial models are much hotter for a given $`L`$. If $`T_{\mathrm{eff}}`$ is high enough for H2 formation to be negligible, such models start to evolve toward cooler $`T_{\mathrm{eff}}`$. Once the fraction of H2 becomes significant and the adiabatic gradient small enough near the photosphere, evolution proceeds almost vertically in the HRD, as also expected from the analysis of Hayashi & Nakano (1963) (see §3.1.1). Such an evolution resembles the shape of the B97 tracks. It may however be a pure coincidence, since B97 do not describe their initial conditions.

Convection in all afore-mentioned models is treated within the framework of the mixing length theory, with $`1<\alpha _{\mathrm{mix}}<2`$, except DM97, who use the Canuto-Mazzitelli formalism. As demonstrated in §3, the treatment of convection affects significantly evolutionary tracks and a variation of $`\alpha _{\mathrm{mix}}`$ modifies the shape of the Hayashi lines (as long as super-adiabatic layers are present). Figure 8 displays quite similar shapes for a given mass between BCAH98, PS99 and DM94, although not at the same $`T_{\mathrm{eff}}`$ because of different outer boundary conditions. The DM97 Hayashi lines behave differently with respect to models of other groups, as already noticed by Siess et al. (2000). This is certainly due to their different treatment of convection. Unexpectedly, the DM97 models predict a MS for VLMS close to the BCAH98 MS, which, as mentioned by DM97, is purely coincidental and stems from unexpected canceling effects between different treatments of convection and outer boundary conditions. This is illustrated for the 0.2 $`M_{}`$ in Figure 8. Since the grey-like outer boundary condition used in DM97 is expected to yield hotter effective temperatures than the BCAH98 models, the close agreement on the MS suggests that the Canuto-Mazzitelli convection treatment yields larger super-adiabatic layers in the atmosphere. This is indeed required to decrease $`T_{\mathrm{eff}}`$ for a given mass. Such a behaviour, however, is not found by the recent simulations by Ludwig et al. (2001) (see §3.2). Indeed, the Canuto-Mazzitelli treatment of convection yields results at odds with 3D hydro simulations for the outer thermal profile of the Sun (see e.g. Nordlund & Stein 1999), and does not provide an accurate treatment of convection in optically-thin media, at least for solar-type stars and low-mass stars ( see e.g. §2.1.3 of Chabrier & Baraffe, 2000 for a discussion of this topic).

## 5 Observational tests

In the previous sections, we emphasize the large uncertainty of evolutionary models at early ages and low gravities. Testing different initial conditions, we find out that after a few Myr these uncertainties become inconsequential and all our models converge toward the same position in a HRD at a given age, for a given mass. Such a result does not necessarily imply that the models are reliable at ages of a few Myr, since we have only tested simple cases. More sophisticated initial conditions are beyond the scope of the present paper and were already explored by Hartman et al. (1997), illustrating the sensitivity of tracks and birthline positions to the (poorly known) details of the protostellar accretion process (geometry, rate, temperature of added matter).

A better knowledge of initial conditions may come from the determination of the minimum age below which present models start to depart significantly from observations. Estimation of this age can constrain the characteristic time-scales and accretion rates of the protostellar collapse phase. Unfortunately, direct comparisons of observations with models directly in colour - magnitude diagrams are extremely uncertain due to the large extinction in star formation regions, which affects the observed energy distribution and thus the spectra and the colors. Only very few exceptions, such as $`\sigma `$ Orionis, exhibit low extinction. Recently, Béjar et al. (1999) and Zapatero et al. (1999, 2000) obtained optical and near-IR photometry for low mass objects in this cluster. In a $`(IJ)`$ vs $`M_\mathrm{I}`$ CMD, the data lie between the 1 and 10 Myr isochrones, respectively, for masses down to $`0.01M_{}`$, using the BCAH98 and CBAH00 models (Zapatero et al. 2000; Béjar et al. 2001). If statistics is improved and if the membership of the objects to the cluster is confirmed, such observations provide an unique opportunity to test directly the validity of young theoretical isochrones. They also offer the best chance to determine the Initial Mass Function (IMF) down to the substellar regime and the minimum mass formed by a collapse process (see Béjar et al. 2001).

Young multiple systems provide also excellent tests for PMS models at young ages, because of the assumed coevality of their different components. In addition, another strong constraint is supplied by the estimate of dynamical masses deduced either from binary systems (Covino et al. 2000; Steffen et al. 2001) or determined from the orbital motion of circumstellar/circumbinary disks (Simon, Dutrey & Guilloteau 2000). An example is provided by the quadruple system GG TAU (White et al. 1999), with components covering the whole mass-range of VLMS and BDs from 1 $`M_{}`$ to $``$ 0.02 $`M_{}`$. Orbital velocity measurements of the circumbinary disk surrounding the two most massive components imply a constraint on their combined stellar mass (Dutrey, Guilloteau & Simon 1994; Guilloteau, Dutrey & Simon 1999). This mass constraint and the hypothesis of coevality provides a stringent test for PMS models. The BCAH98 models are the most consistent with GG Tau (for details see White et al 1999; Luhman 1999) and provide the closest agreement with derived masses of other young systems (see Fig.9).

Most of the observed systems displayed in Figure 9 are better reproduced by tracks using a large value of $`\alpha _{\mathrm{mix}}`$ (=1.9). However, for some systems, such as 1 (Covino et al. 2000), 2 (Steffen et al. 2001) and 4 (BP Tau from Simon et al. 2000), a better agreement is obtained with $`\alpha _{\mathrm{mix}}`$ = 1. Although a variation of $`\alpha _{\mathrm{mix}}`$ with effective temperature and gravity is possible, as suggested by the simulations of Ludwig et al. (1999), none of these three systems occupies a peculiar position in ($`T_{\mathrm{eff}}`$, $`g`$) to suggest a different value of $`\alpha _{\mathrm{mix}}`$. This puzzle may reflect the uncertainties of PMS models based on arbitrary initial conditions. It may also be due to the large uncertainties of observationally-derived spectral type classifications, luminosity estimates and $`T_{\mathrm{eff}}`$ calibrations for such very young objects. Since they display spectral features between that of giants and dwarfs (see Luhman 1999), a better representation of their spectral properties may require new indices more appropriate to these intermediate surface gravities, in the same vein as the pseudo-continuum ratios used by Martín et al. (1996) for Pleiades objects. The transformation of the inferred spectral type into $`T_{\mathrm{eff}}`$ is even more difficult, because of the lack of reliable $`T_{\mathrm{eff}}`$ \- scales for such young T-Tauri like objects. Significant efforts were devoted within the past recent years to the elaboration of improved $`T_{\mathrm{eff}}`$ \- scales for M-dwarfs (Leggett et al. 1996) and M-giants (Perrin et al. 1998; van Belle et al. 1999). However, work remains to be done for T-Tauri like objects. Luhman (1999) defined a $`T_{\mathrm{eff}}`$ \- scale intermediate between giants and dwarfs and based on the isochrone of BCAH98 which goes through the four components of GG Tau. Interestingly enough, applying this $`T_{\mathrm{eff}}`$ \- scale to young clusters such as IC348 (Luhman 1999) and star forming regions like Chamaeleon I (Comerón, Neuhäuser & Kaas 2000), the cluster members show a small scatter in age and no obvious correlation between age and mass. As mentioned by Comerón et al. (2000), this suggests in Chamaeleon I an almost coeval population which formed within less than 1 Myr. This short timescale supports the suggestion that star formation is controlled primarily by large-scale turbulent flows rather than by magnetic processes such as ambipolar diffusion (Hartmann 2000).

Given the assumption of coevality and the uncertainties of our initial conditions, the $`T_{\mathrm{eff}}`$ \- scale for young objects suggested by Luhman (1999) needs to be confirmed in order to confirm these exciting results about star formation process. Although still very preliminary, the comparison of observed and synthetic spectra, as recently done by Lucas et al. (2001) for Orion objects, provides a promising way to define such a $`T_{\mathrm{eff}}`$ \- scale.

## 6 Conclusion

The very good agreement of models based on improved physics with observations for relatively old (t $`>`$$``$ 100 Myr) low-mass objects yields confidence in the underlying theory. Such evolutionary models can now be confronted to the complex realm of very young objects, providing important information on star formation processes and initial conditions for PMS models. Although based on extremely simple initial conditions (no accretion phase, no account of protostellar collapse phase and time scale, spherical symmetry), these models provide the most accurate comparison with present observations of very young objects (dynamical masses, tests of coevality in multiple systems, CMDs). Given the combining effects of large observational and theoretical uncertainties at very young ages, however, one must remain cautious. It is probably too premature to conclude on the validity of the present models at early phases of evolution. We have examined in the present paper the uncertainties on evolutionary models of very-young low-mass objects arising from initial conditions, in particular the initial radius of the object, the efficiency of convection in the outermost layers and the initial abundance of deuterium. We have shown that at least the two first afore-mentioned uncertainties can affect drastically the fundamental properties, luminosity and effective temperature of objects younger than about 1 Myr. Therefore, any attempt to infer an age or a mass from observable quantities for these objects, in particular the initial mass function of very young clusters, must be considered with highly limited - if any - validity !

Realistic initial conditions can only be provided by multi-dimensional protostar collapse simulations, not by spherically-symmetric models for PMS initial conditions involving many free, ill- or unconstrained parameters. Because of numerical subtleties and complex physical processes (accretion fronts, turbulent time-dependent convection, hydrodynamical radiative transfer, magnetic field etc…), the construction of star formation models is a harsh task, which very likely will necessitate several years of efforts. Besides these theoretical difficulties, observations of very young objects can provide only limited guidance to such simulations, since most phases involved during the collapse are embedded in dusty cocoons. Only the final product can be observationally tested, when the protostar becomes visible. This stage marks essentially the beginning of PMS evolutionary tracks. PMS tracks tested against observations thus provide a precious link to gather insight about star formation models from subsequent evolution.

To progress in the field and in parallel with the development of star formation models, efforts can be directed toward: (1) a reduction of the main theoretical uncertainties affecting PMS models at low gravities, which involves in particular a better determination of $`l_{\mathrm{mix}}`$ through multi-D hydrodynamical simulations and (2) the elaboration of a reliable $`T_{\mathrm{eff}}`$-scale for young objects.

Note:  Tracks and isochrones for $`t`$ 1 Myr of the BCAH98 models (from 0.02 $`M_{}`$ to 1.4 $`M_{}`$) and of the CBAH00 models including dusty atmospheres (from 0.001 $`M_{}`$ to 0.1 $`M_{}`$) are available by anonymous ftp:

ftp ftp.ens-lyon.fr

username: anonymous

ftp $`>`$ cd /pub/users/CRAL/ibaraffe

ftp $`>`$ get README

ftp $`>`$ get BCAH98\_models.\*

ftp $`>`$ get BCAH98\_iso.\*

ftp $`>`$ get DUSTY00\_models

ftp $`>`$ quit

###### Acknowledgements.

We are indebted to Lee Hartmann for valuable discussions. I.B thanks the Max-Planck Institut for Astrophysik in Garching for hospitality during elaboration of part of this work. The calculations were performed using facilities at Centre d’Etudes Nucléaires de Grenoble.


format: text

text:
# I Introduction

## I Introduction

Recently Jeon, Koch  and Asakawa, Heinz and Mueller  argued that the study of event-by-event fluctuations of electric charge in high energy nucleus-nucleus collisions may provide information on the state of matter in an early stage of the collision. The authors calculated the magnitude of the charge fluctuations in the Quark-Gluon Plasma (QGP) and in a hadron gas. These calculations show that the fluctuations in the QGP should be significantly smaller than in a hadron gas. Thus they can be used as a signal for the creation of the deconfined phase providing that the initial fluctuations survive hadronization and their relaxation time is significantly longer than the time of the hadronic stage of the collisions .

The charge fluctuations should also be sensitive to other effects like the number of resonances at chemical freeze-out  and fluctuations occurring in the initial stage . Consequently their analysis is interesting also beyond the QGP hypothesis.

First experimental results on charge fluctuations in central collisions of heavy nucleus at SPS  and RHIC  are already available and don’t show the expected reduction of fluctuations due to QGP creation.

The analysis of charge fluctuations is relatively easy to perform experimentally due to the typically good resolution of the measurement of electric charge in tracking detectors positioned in a magnetic field.

The magnitude of the measured charge fluctuations is, however, also dependent on trivial, uninteresting effects, which may shadow the studied physics. The two most important of these effects are:

\- the fluctuations in the event multiplicity, introduced mostly by the variation of the impact parameter;

\- changes in the mean multiplicity due to changes of the colliding system (e.g. p+p, Pb+Pb), changes of collision energy and due to changes of the size of the acceptance in which fluctuations are studied.

In order to minimize the sensitivity to these effects two measures of fluctuations were proposed: $`\stackrel{~}{D}`$  and $`\mathrm{\Phi }`$ .

The aim of this paper is to study the properties of these measures using simple toy models of particle production in nuclear collisions.

In Section II the fluctuation measures will be introduced. Their properties will be tested in Section III. Summary and conclusions are given in Section IV.

## II MEASURES OF CHARGE FLUCTUATIONS: $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$

The measure $`\stackrel{~}{D}`$ is defined as follows  :

$$\stackrel{~}{D}=\frac{\delta R^2N_{ch}}{C_yC_\mu },$$
(1)
where:

$$R=\frac{N_{pos}}{N_{neg}},$$
(2)
$$\delta R=RR,$$
(3)
$$N_{ch}=N_{pos}+N_{neg},$$
(4)
$$C_y=1\frac{N_{ch}}{N_{ch}_{tot}},C_\mu =\frac{N_{pos}^2}{N_{neg}^2}.$$
(5)
$`N_{pos}`$ and $`N_{neg}`$ are measured multiplicities of positively and negatively charged particles within the analysed acceptance. The symbol $`\mathrm{}`$ represents averaging over events. $`N_{ch}`$ is the mean number of charged particles within the acceptance and $`N_{ch}_{tot}`$ denotes the mean total event multiplicity of charged particles in full phase space. The factors $`C_y`$ and $`C_\mu `$ were introduced to remove the influence of global charge conservation ($`C_y`$) and effect of non-zero net charge ($`C_\mu `$). Under the assumption of Poissonian distribution of $`N_{pos}`$ and $`N_{neg}`$ (providing the number of accepted particles is much smaller than the total number of particles $`C_y1`$) and assuming zero net charge ($`Q=N_{pos,tot}N_{neg,tot}=0C_\mu =1`$) one obtains $`\stackrel{~}{D}=4`$ . Many hadrons measured in the final state originate from the decay of resonances. Therefore those are correlated in kinematic quantities like rapidity. Consequently if we consider an acceptance window in a kinematic variable which is large compared to the mean separation of the decay products in this variable the probability that all decay products will fall into it is high and therefore the charge fluctuations should be reduced significantly . In the QGP phase the unit of charge is $`\frac{1}{3}`$ instead of $`1`$ as in the hadron phase. Smaller charge units cause smaller charge fluctuations. Assuming that the QGP fluctuations are frozen one expects $`\stackrel{~}{D}1`$ .

A well established measure of event-by-event fluctuations is the variable $`\mathrm{\Phi }`$  which is defined in the following way:

$$\mathrm{\Phi }=\sqrt{\frac{Z^2}{N}}\sqrt{\overline{z^2}},$$
(6)
where:

$$z=x\overline{x},Z=\underset{i=1}{\overset{N}{}}(x_i\overline{x}).$$
(7)
$`x`$ is a single particle variable, $`N`$ is the event multiplicity within acceptance and overline denotes averaging over a single particle inclusive distribution. By construction, for a system which is an independent sum of identical particle sources the value of $`\mathrm{\Phi }`$ is equal to the value of $`\mathrm{\Phi }`$ for a single particle source independent of the number of superimposed particle sources and its distribution . In the original proposal  $`\mathrm{\Phi }`$ was introduced as a measure of fluctuations of a kinematic variable like transverse momentum. Furthermore it was proposed  to use $`\mathrm{\Phi }`$ as a measure of fluctuations of discrete quantities. Following this suggestion in this publication $`x`$ is identified with the particle charge $`q`$. Some properties of the $`\mathrm{\Phi }`$ measure for discrete quantities were derived in Ref. . For a scenario in which particles are correlated only by global charge conservation ($`gcc`$) and providing $`Q=0`$ we calculated the value of $`\mathrm{\Phi }_q`$:

$$\mathrm{\Phi }_{q,gcc}=\sqrt{1P}1,$$
(8)
where

$$P=\frac{N_{ch}}{N_{ch}_{tot}}.$$
(9)
Thus the value of $`\mathrm{\Phi }_{q,gcc}`$ depends only on the fraction of accepted particles. In order to remove the influence of $`gcc`$ we suggest therefore to use a difference:

$$\mathrm{\Delta }\mathrm{\Phi }_q=\mathrm{\Phi }_q\mathrm{\Phi }_{q,gcc}.$$
(10)
Consequently the value of $`\mathrm{\Delta }\mathrm{\Phi }_q`$ is zero if the particles are correlated only by global charge conservation. It is negative in case of an additional correlation between positively and negatively charged particles and positive if there are anticorrelations. One possible source of additional correlations could be the creation of a QGP. The first part in the definition of $`\mathrm{\Phi }`$ takes into account all, statistical and dynamical fluctuations. The second part quantifies only the statistical ones. Therefore the value of $`\mathrm{\Phi }`$ calculated for different systems is independent of their degrees of freedom (quarks and gluons or hadrons) of the considered system if its fluctuations are only statistical. However, experimentally we measure the fluctuations of particles only after hadronization. A hadronization which conserves charge and entropy produces correlations. Consequently the measured value of $`\mathrm{\Phi }_q`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ should be smaller, if a QGP was created in the early stage of the collision.

## III Toy Models of Nuclear Collisions

In this section properties of the fluctuation measures, $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$, are discussed. First we will assume that the only source of correlations is the global charge conservation. Under this assumption we will study the behaviour of the measures at small $`N_{ch}`$ and their response to non-zero net charge and fluctuations of event multiplicity. Furthermore we also calculate their values in two toy models, a ’$`\rho `$’–gas model and a QGP model in which additional correlations are present. In the study we use a simple Monte Carlo (MC) procedure, which simulates the particle distribution in momentum space. The values of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ are calculated for particles falling into rapidity windows centred around mid-rapidity. In the MC simulations we used values of total multiplicity ($`N_{ch}_{tot}`$) and net charge ($`Q`$) which are typical for central Pb + Pb collisions at 40 A$``$GeV as measured by the NA49 collaboration . In the figures the measures are plotted as a function of $`N_{ch}`$, which is a positive monotonic function of the size of the rapidity window. The ranges of the measures $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ plotted on the vertical scales are chosen in order to make both measures similar sensitive to resonance decay kinematics and QGP creation. Statistical errors of the MC simulations are plotted but in most of the cases they are smaller than the symbol size.

### A Model A

In the first model we analyse the dependence of the measures $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ on $`N_{ch}`$. It is assumed that the total multiplicity is fixed, $`N_{ch,tot}=850`$, and the net charge in each event is zero ($`Q=0`$). The dependence of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ on $`N_{ch}`$ is shown in Fig. 1.

As expected the value of $`\mathrm{\Delta }\mathrm{\Phi }_q`$ is zero and the value of $`\stackrel{~}{D}`$ is close to 4 for large $`N_{ch}`$. However for $`N_{ch}`$ smaller than 200 a significant deviation of $`\stackrel{~}{D}`$ from 4 is observed. It is because the assumptions used to derive the value $`\stackrel{~}{D}=4`$  are not fulfilled in the limit $`N_{ch}0`$.

### B Model B

Both measures are constructed in order to be independent of the fluctuations of $`N_{ch,tot}`$. To test it we varied the total number of charged particles from event to event according to a Gaussian distribution with a mean value 850 and a dispersion of 60 and 120. As in the model A, the global net charge was assumed to be zero. The results of this simulation are shown in Fig. 2.

The behaviour of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ is similar to this obtained in Model A. Thus the influence of multiplicity fluctuations is small, especially on $`\mathrm{\Delta }\mathrm{\Phi }_q`$.

### C Model C

After we have observed that the measures are insensitive to fluctuations of $`N_{ch,tot}`$ we turn to check the influence of non-zero net charge. We assume $`N_{ch,tot}=850`$ and $`Q=60`$. Fig. 3 shows the results.

The measure $`\stackrel{~}{D}`$ is explicitly corrected for non-zero net charge by the factor $`C_\mu `$, but it still shows a significant dependence on it. Much weaker dependence on the net charge is observed for $`\mathrm{\Delta }\mathrm{\Phi }_q`$. For an acceptance larger than 50% its value is somewhat higher than the calculated one assuming $`Q=0`$.

### D Model D

In nuclear collisions the number of charged particles fluctuates and the net charge is approximately proportional to $`N_{ch,tot}`$. To see how the measures behave under these assumptions we used a model in which $`N_{ch,tot}`$ was generated from a Gaussian distribution ($`N_{ch,tot}=850`$, $`\sigma (N_{ch,tot})=80`$) and $`N_{pos,tot}=\frac{7}{12}N_{ch,tot}`$, $`N_{neg,tot}=\frac{5}{12}N_{ch,tot}`$ . The results of this simulation are presented in Fig. 4.

The value of $`\stackrel{~}{D}`$ is significantly higher than 4 in the full range of $`N_{ch}`$. $`\mathrm{\Delta }\mathrm{\Phi }_q`$ is much less dependent on the studied effects. It deviates slightly from 0 for $`N_{ch}>0.5N_{ch,tot}`$.

### E QGP Model

Following Ref.  we construct a simple toy model of the QGP. Under the assumption of zero baryo-chemical potential and providing zero quark mass we calculate the ratio of up-, down-, antiup- and antidown- quarks and gluons in equilibrium. Assuming entropy and net charge conservation in each rapidity window during the evolution from the QGP to the final hadron state the number of pions ($`N`$) and their net charge is calculated. The number of charged pions is taken to be $`N_{ch}=\frac{2}{3}N`$ based on isospin symmetry. Using this model we obtain the results shown in Fig. 5.

Our simulation confirms results published in Ref. , the $`\stackrel{~}{D}`$ value in our QGP toy model is close to $`1`$. The new result is the calculation of $`\mathrm{\Delta }\mathrm{\Phi }_q`$ in this model. Its value is mostly significantly smaller than zero. It increases from $`0.5`$ for small $`N_{ch}`$ to $`0.1`$ for the largest multiplicity. In case all particles are used, there is no net charge fluctuations due to charge conservation and therefore the value of $`\mathrm{\Delta }\mathrm{\Phi }_q`$ should be $`0`$. We note that the previously discussed trivial effects were increasing $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ to values above $`4`$ and $`0`$, respectively. Thus they can not mimic suppression of fluctuations observed in the QGP toy model.

### F $`\rho `$’–gas Model

In this model we study the influence of the resonance decay kinematics on $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$. We use a MC simulation which generates 400 neutral resonances (’$`\rho `$’) in each event. The mass of the ’$`\rho `$’ is fixed to 770 MeV/$`c^2`$. The ’$`\rho `$’ rapidity distribution is Gaussian with $`\sigma (y)=0.8,1,2`$ or $`3`$ and the transverse mass distribution is ’thermal’: $`\frac{d^2n}{dm_Tdy}=Cm_Te^{m_T/T}`$ with $`T=170`$ MeV. All ’$`\rho `$’s decay into two charged pions which are then used to calculate values of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$. In Fig. 6 we show the results of a simulation for $`\sigma (\mathrm{y})=1.0`$.

For low $`N_{ch}`$ the values of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ are higher then 4 and 0, respectively. With increasing $`N_{ch}`$ these values decrease. Resonance decay kinematics influences charge fluctuations in two different ways. It dilutes the effect of global charge conservation if only one of the decay products falls into the acceptance. This will increase the values of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$. If both decay products fall into the acceptance $`N_{ch}`$ increases but the net charge doesn’t change. $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ measure charge fluctuations normalised by the number of charged particles. Therefore their values will drop in this case. Which effect dominates depends on the size of acceptance window ($`\mathrm{\Delta }\mathrm{y}`$, upper scale on the plots), mass of the decaying resonances and their rapidity distribution. The sensitivity on the last one is shown in Fig. 7.

At SPS energies $`\sigma (\mathrm{y})`$ of the $`\rho `$ distribution is expected to be close to $`1`$ . In Figs. 6 and 7 (upper plots) we see that small changes of $`\sigma (\mathrm{y})`$ around $`1`$ have strong influence on the values of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$. For $`\sigma (\mathrm{y})=1`$ we observe a transition from fluctuation enhancement to suppression with an increasing rapidity window. For $`\sigma (\mathrm{y})=0.8`$ the values of $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ are above or equal to their values expected for independent particle production with global charge conservation for all $`\mathrm{\Delta }\mathrm{y}`$. We would like to stress here that in our model the ’$`\rho `$’s are neutral and therefore there are no charge fluctuations before the ’$`\rho `$’s decay. In Fig. 6 and Fig. 7 (upper plots) we see that for $`\sigma (\mathrm{y})`$ values expected at CERN SPS the resonance decay kinematics increases the initial zero charge fluctuations to values which are close to 4 or 0, respectively. The fluctuations are significantly lower for larger values of $`\sigma (\mathrm{y})`$, e.g. $`\sigma (\mathrm{y})=2`$ and $`3`$, as are expected at RHIC and LHC energies. In the case of QGP creation the measured fluctuations will be composed of the initial QGP fluctuations and fluctuations produced by the decay of resonances. The latter ones are, as our model shows, still larger than the predicted QGP fluctuations even for $`\sigma (\mathrm{y})=2`$ and $`3`$. Therefore we conclude that resonance decay kinematics may significantly shadow charge fluctuations developed at the early stage of the collision.

## IV Summary and Conclusions

The properties of two measures of charge fluctuations $`\stackrel{~}{D}`$ and $`\mathrm{\Delta }\mathrm{\Phi }_q`$ were studied within several models. We have shown the dependence of both measures on trivial effects, in particular, on number of particles used per event, on the net electric charge and on fluctuations in the total event multiplicity. Furthermore their response to resonances decay kinematics and to the reduction of fluctuations due to QGP creation were analysed, too. From our simulations we learn that the measure $`\mathrm{\Delta }\mathrm{\Phi }_q`$ is significantly less dependent on trivial effects than $`\stackrel{~}{D}`$ and that both measures are sensitive to frozen QGP fluctuations and to resonance decay kinematics. Therefore we suggest to use $`\mathrm{\Delta }\mathrm{\Phi }_q`$ as a proper measure of charge fluctuations. Using the ’$`\rho `$’ gas model we have shown that the influence of resonance decay strongly depends on their rapidity distribution and on the acceptance window. Given the typical widths of the rapidity distributions of $`\rho `$’s at CERN SPS energies they can easily increase charge fluctuations from QGP values to values expected for independent particle production with global charge conservation. A further development of methods which allow to distinguish various sources of charge fluctuations is obviously necessary.

Acknowledgements

I am very indebted to Marek Gaździcki for numerous fruitful discussions and ideas. I’m also grateful to S. Mrówczyński, V. Koch, H. Ströbele, R. Stock, P. Seyboth, R. Renfordt and J. Berschin for helpful comments.


text:
# The energy states of quantum dots through the variational method and SQM*footnote **footnote *Talk given at the VII International Wigner Symposium, University of Maryland, College Park, USA, August 2001

## I Introduction

Twenty years ago Supersymmetric Quantum Mechanics (SQM) appeared inside the study of the SUSY breaking mechanism of higher dimensional quantum field theories, . Since then it has been considered as a new field of research, providing not only a supersymmetric interpretation of the Schr$`\ddot{o}`$dinger equation, but interesting answers to all sorts of non-relativistic quantum mechanical systems. Particular points to be mentioned include the better understanding it brought of the exactly solvable, -, the partially solvable, , , the isospectral, , , the periodic,  and the non-exactly solvable potentials, -. The latter were studied by the association of the variational method with SQM formalism. The scheme employed in them is based on an ansatz made to the superpotential. Through the superalgebra it is possible to evaluate the wave function, the trial wave function, which contains the free parameters that minimize the energy expectation value.

We use this methodology to study the confined hydrogen atom, a system described by the confined Coulomb potential. The confining approach is the same as the one used to study the problem of a hydrogenic donor located at the centre of a spherical $`GaAs(Ga,Al)As`$ quantum dot, a semiconductor device that confines electrons, , . Here the energy of the $`1s`$, $`2p`$ and $`3d`$ states of such confined system is evaluated and the results are compared to other variational and exactly numerical results, .

## II The variational method associated to SQM

Consider a system described by a given potential $`V_1`$. The associated Hamiltonian $`H_1`$ can be factorized in terms of bosonic operators, in $`\mathrm{}=c=1`$ units, -.

$$H_1=\frac{1}{2}\frac{d^2}{dr^2}+V_1(r)=A_1^+A_1^{}+E_0^{(1)}$$
(1)
where $`E_0^{(1)}`$ is the lowest eigenvalue. Notice that the function $`V_1(r)`$ includes the barrier potential term. The bosonic operators are defined in terms of the so called superpotential $`W_1(r)`$,

$$A_1^\pm =\frac{1}{\sqrt{2}}\left(\frac{d}{dr}+W_1(r)\right).$$
(2)
As a consequence of the factorization of the Hamiltonian $`H_1`$, the Riccati equation must be satisfied,

$$W_1^2W_1^{}=2V_1(r)2E_0^{(1)}.$$
(3)
Through the superalgebra, the eigenfunction for the lowest state is related to the superpotential $`W_1`$ by

$$\mathrm{\Psi }_0^{(1)}(r)=Nexp(_0^rW_1(\overline{r})𝑑\overline{r}).$$
(4)
At this point we address to the variational method. It was conceived to be an approximative method to evaluate the energy spectra of a given Hamiltonian $`H`$ and, in particular, its ground state. Its central point is the search for an optimum wave-function $`\mathrm{\Psi }(r)`$ that depends of a number of parameters. This is called the trial wave-function. The approach consists in varying these parameters in the expression for the expectation value of the energy

$$E=\frac{\mathrm{\Psi }^{}H\mathrm{\Psi }𝑑r}{\mathrm{\Psi }^2𝑑r}$$
(5)
until it reaches its minimum value. This value is an upper limit of the energy level pursued.

Thus, using physical arguments an ansatz is made to the superpotential and, through the superalgebra, the wave function is evaluated. This is our trial wave function which depends on a set of the free variational parameters, which were introduced by the ansatz to the superpotential.

Notice that since the potential is non-exactly solvable, the Hamiltonian is not exactly factorizable, in other words, there is no superpotential that satisfies the Riccati equation exactly. However, the superpotential we have in hands through the ansatz does satisfy the Riccati equation by an effective potential, $`V_{eff}`$

$$V_{eff}(y)=\frac{\overline{W}_1^2\overline{W}_1^{}}{2}+E(\overline{\mu })$$
(6)
where $`\overline{W}_1=W_1(\overline{\mu })`$ is the superpotential that satisfies (3) for $`\mu =\overline{\mu }`$, the set of parameters that minimise the energy of eq.(5).

## III The Confined Coulomb Potential

The Hamiltonian of an on-center impurity in a spherical quantum dot can be written in the effective-mass approximation as

$$H=\frac{\mathrm{}^2}{2m^{}}^2\frac{e^2}{ϵr}$$
(7)
where $`m^{}`$ is the effective mass and $`ϵ`$ is the dielectric constant of the material of the quantum dot. The donor is assumed to be at the centre of the quantum dot of radius $`R`$ with an infinite barrier height. This means that the wave function vanishes at $`r=R`$.

In atomic units, the radial Hamiltonian equation for the Coulomb Potential is given by

$$H=\frac{d^2}{dr^2}+\frac{l(l+1)}{r^2}\frac{2}{r}.$$
(8)
We use the variational method associated with the SQM in order to get the energy states of the confined atom.

We make the following ansatz to the superpotential

$$W(r)=\frac{c}{r}+b+\frac{a}{Rr}$$
(9)
which depends on three free parameters, the variational parameters. From the superalgebra, the eigenfunction obtained from equation (4) is given by

$$\mathrm{\Psi }(a,b,c,r)e^{br}r^c(Rr)^a.$$
(10)
This expression is a trial wavefunction for the variational method having three parameters, $`a`$, $`b`$ and $`c`$ and vanishing at $`r=R`$.

The energy is obtained by minimisation of the energy expectation value with respect to the free parameters. The equation to be minimised is given by the following evaluation

$$E(a,b,c)=\frac{_0^{Rr}\mathrm{\Psi }(a,b,c,r)[\frac{d^2}{dr^2}\frac{2}{r}+\frac{l(l+1)}{r^2}]\mathrm{\Psi }(a,b,c,r)𝑑r}{_0^{Rr}\mathrm{\Psi }(a,b,c,r)𝑑r}.$$
(11)
In support of the choice made in (9) for the superpotential are the physical arguments: the first two terms of (9) are related to the exact Coulomb case; the last term is the confining term. This is remarked by the form of the effective potential $`V_{eff}`$ that satisfies the Riccati equation, (eq. 6), evaluated at the values of $`a`$, $`b`$, and $`c`$ that minimize the energy expectation value,

$`V_{eff}(r,a,b,c)={\displaystyle \frac{c(c1)}{r^2}}+{\displaystyle \frac{a(a1)}{(Rr)^2}}{\displaystyle \frac{2bc}{r}}{\displaystyle \frac{2ac}{r(Rr)}}+{\displaystyle \frac{2ab}{Rr}}+b^2+\overline{E}(a,b,c)`$ (12)
which is infinite at $`r=R`$.

## IV Results

Thus, minimizing the equation (11) with respect to the parameters $`a`$, $`b`$ and $`c`$ we find the minimum value of the energy for different values of $`R`$ and $`l`$. These results are shown in the tables bellow. Comparison is made with results from , which are based on the standard variational method in which a variational wavefunction is proposed. Ref.  also contains exact numerical results.

Table 1. Energy eigenvalues (in Rydbergs) for different values of $`R`$ for $`l=0`$ with a, b and c as variational parameters. Comparison is made with results of Ref. .
R $`E_{EXACT}`$ $`E_V`$ $`\underset{¯}{|E_{EXACT}E_V|}\%`$ $`E_{SQM}`$ $`\underset{¯}{|E_{EXACT}E_{SQM}|}\%`$   Ref.  Ref.  $`E_{EXACT}`$  $`E_{EXACT}`$  0.1 937.986 937.999 0.00 940.688 0.29  0.5 29.496 29.497 0.00 29.571 0.25  1 4.7480 4.7484 0.01 4.7565 0.18  2 -0.25000 -0.24990 0.04 -0.25000 0.00  3 -0.84793 -0.84523 0.32 -0.84706 0.10  4 -0.96653 -0.95518 1.17 -0.96509 0.15

Table 2. Energy eigenvalues (in Rydbergs) for different values of $`R`$ for $`l=1`$ with a, b and c as variational parameters. Comparison is made with results of Ref. .
R $`E_{EXACT}`$ $`E_V`$ $`\underset{¯}{|E_{EXACT}E_V|}\%`$ $`E_{SQM}`$ $`\underset{¯}{|E_{EXACT}E_{SQM}|}\%`$   Ref.  Ref.  $`E_{EXACT}`$  $`E_{EXACT}`$  0.4 116.896 116.925 0.02 117.038 0.12  2.0 3.1520 3.1530 0.03 3.1555 0.11  4.0 0.28705 0.28706 0.00 0.28732 0.09  6.0 -0.111111 -0.11069 0.38 -0.111111 0.00  7.0 -0.17496 -0.17392 0.59 -0.17490 0.03  8.0 -0.20890 -0.20691 0.95 -0.20882 0.04

Table 3. Energy eigenvalues (in Rydbergs) for different values of $`R`$ for $`l=2`$ with a, b and c as variational parameters. Comparison is made with results of Ref. .
R $`E_{EXACT}`$ $`E_V`$ $`\underset{¯}{|E_{EXACT}E_V|}\%`$ $`E_{SQM}`$ $`\underset{¯}{|E_{EXACT}E_{SQM}|}\%`$   Ref.  Ref.  $`E_{EXACT}`$  $`E_{EXACT}`$  1.0 29.935 29.950 0.05 29.958 0.08  3.0 2.5856 2.5867 0.04 2.5876 0.08  4.0 1.2427 1.2431 0.03 1.2437 0.08  5.0 0.65823 0.65836 0.02 0.65873 0.08  7.0 0.19318 0.19318 0.00 0.19333 0.08  8.0 0.09212 0.09216 0.05 0.09220 0.08  12.0 -0.06250 -0.06181 1.11 -0.06250 0.00  14.0 -0.08623 -0.08479 1.66 -0.08622 0.01

## V Comments and conclusions

The variational method associated with SQM was applied to get the energy states of a confined hydrogen atom. The problem is similar to the confinement of electrons in a quantum dot. We used a confining effective potential that depends on tree parameters (a, b, c). Minimizing the energy expectation value with respect to these three parameters we found good results, even for increasing values of the radius $`R`$, when compared to results coming from other approximative variational method and exact numerical results. The term $`\frac{1}{(Rr)^2}`$ in the effective potential leads to border effects, since it is infinite at $`R=r`$. These effects do not appear in the original problem, eq. (8) and become more perceptive for smaller values of $`R`$. For increasing values of $`R`$, smaller will be the border effects so that the variational results get better.

In conclusion, we remark that the results presented here once again showed that the association of the superalgebra of SQM with the variational method provides good information about atomic systems.

## VI Acknowledgements

RMR and EDF would like to acknowledge partial financial support from FAPESP and CNPq respectively.