A&A 424, 919-925 (2004)
DOI: 10.1051/0004-6361:20040470

New grids of stellar models including tidal-evolution constants up to carbon burning[*]

I. From 0.8 to 125 $M_{\odot }$ at Z = 0.02

A. Claret

Instituto de Astrofísica de Andalucía, CSIC, Apartado 3004, 18080 Granada, Spain

Received 18 March 2004 / Accepted 11 May 2004

Abstract
We present new stellar models based on updated physics (opacities, expanded nuclear network and mass loss rates). We compute stellar models suitable for the mean solar neighborhood, i.e. for Z=0.02 and X=0.70. The covered mass range is from 0.8 up to 125 $M_{\odot }$ and the models are followed until the exhaustion of carbon in the core, for the more massive ones. In addition, the effective temperatures of the more massive models are corrected for the effects of stellar winds, while models with lower effective temperatures are computed using a special treatment of the equation of state (CEFF). Convective core overshooting is assumed to be moderate and is modelled with $\alpha_{\rm ov} = 0.20$.

Besides the classical ingredients of stellar models, we also provide the internal structure constants needed to investigate apsidal motion and tidal evolution in close binaries. The latter constants are made public for the first time. According to the current theories of tidal evolution, the time scales for synchronization and circularization for cool stars depend - apart from the mass, radius and effective temperature - on the depth of the convective envelope  $x_{\rm bf}$ and on the radius of gyration $\beta$. For stars with higher effective temperatures, these dependencies are mainly incorporated in the tidal torque constant E2. All these parameters are steep functions of mass and time, and thus require a special numerical treatment. The new mass loss formalism produces more mass concentrated configurations than previously, especially for more massive and more evolved stellar models. As the present grid is designed mainly for the study of double-lined eclipsing binaries, the gravity-darkening exponents necessary to calculate the surface brightness distribution in rotationally and/or tidally distorted stars are computed following the method described recently by Claret (1998), and made available for each point of every evolutionary track.

Key words: stars: binaries: close - stars: evolution - stars: interiors - stars: fundamental parameters - stars: abundances - stars: rotation

1 Introduction

Since the first investigations on stellar structure one of the most exciting goals was to look into the interior of stars. Oscillating stars, in principle, offer such an opportunity, but generally do not yet have sufficiently accurate mass determinations. An accurate view inside the interior of a star is therefore not yet possible at this time. Double-lined eclipsing binaries (DLEB), on the other hand, have average errors in masses and radii of the order of 1%, which allow us to investigate stellar interiors using two different techniques: apsidal motion and tidal evolution.

Until very recently, systematic investigations of DLEB were restricted to the comparison between the observed radii, masses and effective temperatures with their theoretical counter-parts derived from evolutionary models. The only serious constraint was the coevity of the components. This situation has now changed drastically: in addition to these primary conditions, present-day studies are also oriented towards the properties of the stellar interior. The series of recent papers by Lacy, Torres and collaborators (Torres et al. 1999, 2000; Lacy et al. 2000, 2002, 2003) is a good example of the new systematics. On the other hand, extragalactic DLEB are now being investigated in different chemical environments such as the Magellanic Clouds or M 31 (for example, Fitzpatrick et al. 2003). Concerning direct investigations of the interior of extragalactic stars, North & Zahn (2003) performed an analysis of the critical fractional radius for the circularization of detached binaries in the SMC and LMC using data from OGLE and MACHO.

On the other hand, the study of rapidly rotating stars requires knowledge of how the flux is distributed over the stellar surface. The brightness distribution depends on the gravity-darkening exponent $\beta _1$. The first investigations are due to von Zeipel (1924) who determined that, for stars with radiative envelopes, the flux is directly proportional to the local gravity, i.e. $T_{\rm eff}^4 \propto g^{\beta_1}$, with $\beta_1= 1.0$. Kippenhahn (1977) and Maeder (1999) later generalized the functional relationship between the flux and the local gravity. More than 40 years after the pioneering work by von Zeipel, Lucy (1967) investigated the case of stars with convective external layers and found that $\beta _1$ was, on average, of the order of 0.3. In 1998, Claret made use of a new numerical method based on the triangle strategy of Kippenhahn to derive $\beta _1$ as a function of mass, age, chemical composition, and effective temperature. These calculations superseded the traditional values of 1 and 0.32 for radiative and convective envelopes, and provided a smooth transition in $\beta _1$ for the two mechanisms of transport of energy. Claret (2000) later also extended the method to be applicable to brown dwarfs. In the latter paper the effect of differential rotation was also considered for more complicated systems, such as those presenting an accretion disk. Since investigations of rapidly rotating stars acquire more and more importance, gravity-darkening exponents are also presented for all models in this paper.

These new challenges support the necessity of ongoing updates of stellar evolution codes to provide new and more accurate theoretical material to be compared with observations. The main objective of the grid of stellar models presented in this paper is therefore to serve as a theoretical tool to explore the commented possibilities.

2 A summary of the input physics of the code

The four main differences between the present code and the 1995 version are: the opacities, the nuclear network, the mass loss rates and the mass range for which the models are computed. For radiative opacities we adopted the calculations by Iglesias & Rogers (1996) completed with Alexander & Ferguson (1994) results for lower temperatures. This set of opacities forced a revision in the calibration of the mixing-length parameter $\alpha$ for the Sun. For X = 0.684 and Y = 0.296, the new calibration yields $\alpha = 1.68$. Core overshooting is considered through the parametrization $d_{\rm ov}= \alpha_{\rm ov}H_{\rm p}$, where $H_{\rm p}$ is the pressure scale height and $\alpha_{\rm ov}$ is the overshooting parameter. For the latter, we adopt the value 0.20 which seems to fit with acceptable accuracy the DLEB (Ribas 1999; Lastennet & Valls-Gabaud 2002; Claret & Willems 2002). Other mixing processes are not considered, like semiconvection or that due to rotation. In particular, semiconvection poses an interesting problem: in an intermediate region above the convective core, the mean molecular weight decreases with the distance and depending on the conditions, the zone should be mixed due to the differences between the radiative and the adiabatic gradients. On the other hand, if this zone is mixed, the chemical profile may change yielding a reverse situation. The equilibrium is achieved by means of a slow mixing (semiconvection), which is really a diffusion process. Semiconvection is important in massive stars but its effects are more notable in later evolutionary phases. We plan to implement semiconvection to investigate the effects of partial mixing in the stellar internal structure, specially the influence on the size of the convective core, as pointed out by Langer (1991).

In order to extend our investigation on stellar evolution of massive binary stars, we have implemented a more extensive nuclear network than before. We follow the evolution of 13 isotopes: 1H, 4He, 12C, 13C, 14N, 16O, 17O, 18O, 20Ne, 22Ne, 24Mg, 25Mg, and 26Mg. A more elaborate network is available as an option in the code in case the production/destruction of light elements as Li and Be or when more advanced stages of nucleosynthesis are required. This optional version is an adaptation of the routine written by Timmes (Timmes, private communication) for a general network with 47 isotopes. The calculations performed with our network and with the adaptation of the TORCH routine are compatible, at least for the evolutionary stages we have explored in this paper. However, this routine is very CPU intensive so that we only use it in special cases.

In our 13 isotope nuclear network we have used the usual reactions of the pp chains and the CNO tri-cycle. For helium burning the following reactions were introduced:

4He( $\alpha,\gamma$)8Be( $\alpha,\gamma$)12C( $\alpha,\gamma$)16O( $\alpha,\gamma$)20 Ne

17O($\alpha$, n)20Ne

20Ne( $\alpha,\gamma$)24Mg

13C($\alpha$, n)16O

22Ne($\alpha$, n)25Mg

22Ne( $\alpha \gamma$)26Mg

14N(( $\alpha,\gamma$)18F($\beta\nu$)18O( $\alpha,\gamma$)22Ne.
Carbon burning was followed using the reactions:

12C(12C, $\alpha$)20Ne

12C(12C, p)23Na

23Na(p, $\alpha$)20Ne

12C( $\alpha,\gamma$)16O

16O( $\alpha,\gamma$)20Ne

20Ne( $\alpha,\gamma$)24Mg.
The basic nuclear reactions rates are from Caughlan & Fowler (1988), while the screening factors are taken from Graboske et al. (1973). The loss of energy due to neutrinos was considered using the work by Itoh et al. (1989). For cooler models, the adopted equation of state was CEFF (Däppen 2000, private communication), while for hotter ones partial ionization was taken into account for the same elements as in CEFF (H, He, C, N, O, Ne, Na, Mg, Al, Si, Ar, Fe) by interacting the Saha equations. The Fermi-Dirac integrals are used where necessary.

For the more massive models, the mass loss rates are the main determinants acting on the evolutionary path of a given model. We adopted the formalism by Nieuwenhuijzen & de Jager (1990) for all models except in two situations: Wolf-Rayet stars and red giants with masses smaller than 4 $M_{\odot }$. In the latter case, we used the work by Reimers (1977). The parametrization by Langer (1989) was used for WNE and WC Wolf-Rayet stars while for WNL Wolf-Rayet stars, we use the mass-loss rates given by Conti (1988).

If one compares our stellar models of 1995 with the present ones, the differences in the tracks are evident, mainly for more massive stars. In part, this is due to the fact that in the previous work the calculations were interrupted during the first stages of helium burning which renders a direct comparison impossible. In addition, more advanced evolutionary stages such as the WR stages were not computed with the 1995 code. The principal differences, however, are the mass range and the implications of the updated physics for the internal structure (see next sections).

Another parameter that influences the calculation of the apsidal-motion constants, the tidal-torque constant, and the depth of the convective envelope is obviously the radius of the stellar configuration. In stars with high mass-loss rates, the optical thickness of the winds are finite. In order to introduce this effect in the code, we follow the recipe given in Schaller et al. (1992), and references therein. The corrections are mainly important for WR stars, and are almost negligible during early evolutive phases. Although such a scheme is consistent (the opacity flux weighted, for example), there are still large uncertainties in the theory of stellar winds. We therefore consider these models to be less reliable than those located near the Main-Sequence.

The theoretical HR diagram corresponding to our grid of evolutionary tracks is represented in Fig. 1. The basic differences between our calculation and those for example by Schaller et al. (1992) are associated with the morphology of the tracks of more massive stars during later evolutionary phases: our tracks predict, generally, cooler configurations. Recent comparisons between several grids of stellar models (Lastennet & Valls-Gabaud 2002) indicate that the models from Geneve, Padova and Granada give very similar results. However, the sample of known DLEB is still limited in mass range as it contains only few systems with components more massive than 25 $M_{\odot }$. Unfortunately, there are no sufficient massive DLEB to help to elucidate the question. Hence, as mentioned above, we should take these advanced phases with caution: only when a systematic comparison with observational data becomes possible (mainly using DLEB) will we be able draw definitive conclusions.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig1.eps}
\end{figure} Figure 1: The theoretical HR diagram. The numbers below the tracks indicate the stellar masses in solar units.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig2.eps}
\end{figure} Figure 2: The evolution of the apsidal-motion constant k2 for 5 and 40 $M_{\odot }$ models. The numbers on the tracks indicate the initial masses (in solar units). The solid lines represent the present calculations, while the dashed ones represent the models of 1995.
Open with DEXTER

3 The apsidal-motion constants

The well known apsidal-motion phenomenon has been used for decades to test evolutionary models (see Claret & Giménez 1993 for a historic review). Briefly, the distortions in the shape of the components of an eccentric binary due to rotation and tides lead to secular changes in the position of the periastron. These distortions can be described in terms of the mass distribution of the stars, and can thus be directly derived from the theoretical stellar models. In addition to the contribution of the stellar shape distortions, there is also a general relativistic contribution which makes some systems (mainly those with high orbital eccentricity) particularly interesting as tests to the theory of General Relativity. The historic discrepancies between theory and observations seems to be reduced to acceptable levels after the introduction of modern stellar models which take into account the effects of dynamic tides (Claret & Willems 2002). However, some systems are yet to be explained (DI Her, for example) and require further refinements in the calculation of the theoretical internal structure constants (k2, k3 and k4).

As usual, the apsidal-motion constants are derived by integrating the Radau equation. Figure 2 shows the comparison between the present calculations (solid lines) and those by Claret (1995) (dashed ones) for 5 and 40 $M_{\odot }$ models. The mass distributions are very similar for both 5 $M_{\odot }$ models, at least up to exhaustion of hydrogen in the core. The present model was computed until later stages, which is reflected by the peak in log k2 around $\log~ g = 1.7$. The situation for the 40 $M_{\odot }$ model is very different. The mass distributions are only comparable during the early evolutionary stages, while during the later ones the present models are much more mass concentrated than the previous ones. Although the new opacities and the expansion of the nuclear network play some role, the main cause of these differences is due to the higher mass-loss rates adopted for the present models (Claret & Giménez 1991 showed that mass loss increases the grade of mass concentration). The differences may achieve 0.5 dex near core-hydrogen exhaustion for the 40 $M_{\odot }$ model. These differences further improve the comparison between theory and observations.

4 The tidal evolution

4.1 The turbulent breaking mechanism

Concerning tidal evolution, the situation is not so clear but this is compensated for the presence of an additional scenario besides the case of isolated binaries: the critical time scales for the circularization of binaries belonging to clusters. Spectroscopic observations (radial velocity curves) of binary systems in clusters indicate the presence of a cut-off period ( $P_{\rm cut}$) below which the binaries present circular orbits (Mathieu et al. 1992, for example). In addition, there seems to be a relationship between the age of the cluster and the cut-off periods. Dealing with cut-off periods of binaries in clusters is a very interesting and iterative exercise since the age of the cluster (which can be used as an upper limit for the critical time scales) needs to be computed,and, simultaneously, the differential equations governing the tidal evolution for the most representative close binaries need to be integrated in order to compare the theoretical  $P_{\rm cut}$ with the observed  $P_{\rm cut}$-age relationship (Claret, in preparation).

The stellar interior/envelope conditions dictate what mechanism will be responsible for tidal braking. For stars with convective envelopes this agent is turbulent dissipation. The associated time scale is of the order of R2/$\nu_{t}$ and the convective friction time is, following Zahn (1966),

\begin{displaymath}{t_{\rm f}}= {\left[{MR^{2}\over{L}}\right]^{1/3}}
\end{displaymath} (1)

where M is the stellar mass, R the radius, and L the luminosity. The differential equations for synchronization and circularization can be written as (Zahn 1984):

                     $\displaystyle {\left(\tau_{\rm sync}\right)_{\rm turb}}{^{-1}}$ = $\displaystyle -\frac{1}{\Omega-\omega}\frac{\rm d\omega}{{\rm d}t}$  
  = $\displaystyle {3.95\times10^{2}}{\beta^2}{M^{7/3}}{{(1+q)^2}\over{q^2}}
{L^{-1/3}}{{\lambda_{2}}^{-1}}{P^{4}\over{R^{16/3}}}$ (2)


                     $\displaystyle {\left(\tau_{\rm circ}\right)_{\rm turb}}{^{-1}}$ = $\displaystyle -\frac{1}{e}\frac{{\rm d}e}{{\rm d}t}$  
  = $\displaystyle {1.99\times10^{3}}{M^{3}}{{(1+q)^{5/3}}\over{q}}{L^{-1/3}}
{{\lambda_{2}}^{-1}}{P^{16/3}\over{R^{22/3}}},$ (3)

where $\Omega$ is the rotational angular frequency, $\omega$ the orbital angular frequency, $\beta$ the radius of gyration, e the orbital eccentricity, q the mass ratio, P the orbital period, and the time scales are expressed in years. The constant $\lambda_2$ depends on the depth of the convective envelope  $x_{\rm bf}$ in normalized units and, when the convective turnover time is small in comparison to the tidal period, is given by the following expression (Zahn 1989):

\begin{displaymath}{\lambda_{2}} = {{0.607}{\alpha^{4/3}}{E^{2/3}}{\int_{x_{\rm bf}}^{1}}
{x^{22/3}}{(1-x)^2}{\rm d}x}.
\end{displaymath} (4)

Here $\alpha$ denotes the mixing-length parameter and E describes the characteristics of a polytropic envelope. In this paper we evaluate the latter constant for each evolutionary configuration by means of the approximation

\begin{displaymath}{E} = {{Q_{\rm conv}}\over{{\int_{x_{b}}^{1}}{{\left({2}{(1-x)}\over{5x}\right)}^
{3/2}}{x^{2}}{{\rm d}x}}}
\end{displaymath} (5)

where Q $_{\rm conv}$ is the mass contained in the convective envelope.

The parameter $\lambda_2$ is a steep function of the effective temperature. In fact, the precision in the determination of the depth of the convective envelopes depends strongly on the convergence criterions used to compute the stellar models. Moreover, the triangle used to define an envelope in the HR diagram must be as small as possible to avoid numerical oscillations in  $x_{\rm bf}$. We typically use $\Delta$log  $T_{\rm eff} = 0.001$ and $\Delta$log L = 0.004 or smaller, depending on the numerical conditions. This numerical consistency is, however, only possible at the expense of more CPU time. Even taking into account the mentioned precaution, we still detected small numerical oscillations but these do not have a significant influence on the time scales. Figure 3 shows the behavior of the depth of the convective envelope (0 means a fully convective model) as a function of the effective temperature for a 1.41 $M_{\odot }$ model. The "zig-zag'' around log  $T_{\rm eff} \approx 3.8$ is due to the red hook. The constant E is represented in Fig. 4 for the same model and presents similar features: the smaller the depth of the convective envelope, the smaller the corresponding values of E. The integral which appears in Eq. (4) is a simple one and users can easily compute the time scales and thus integrate the differential equations for e and $\omega$ directly from the files providing $\beta, E$ and  $x_{\rm bf}$. As a final note on the turbulent dissipation mechanism, it is worth to mention that $\lambda_2$, as expected, also depends on the chemical composition. Such a dependency may have some implications for the levels of synchronization and circularization of some problematic systems.

  \begin{figure}
\par\includegraphics[height=8cm,width=8.8cm,clip]{fig3.eps}
\end{figure} Figure 3: The depth of the convective envelope (normalized radius) as a function of the effective temperature for a 1.41 $M_{\odot }$ model.
Open with DEXTER


  \begin{figure}
\par\includegraphics[height=8cm,width=8.8cm,clip]{fig4.eps}
\end{figure} Figure 4: The constant E characterizing a polytropic envelope (Eq. (5)) as a function of the effective temperature. Same models as in Fig. 3.
Open with DEXTER

4.2 The radiative damping and the calculation of the tidal torque constant

The case of more massive stars with convective cores and radiative envelopes is somewhat different. The dominant mechanism for tidal braking in this kind of stars is radiative damping, characterized by the tidal torque constants En. The time scales for synchronization and circularization are:

\begin{displaymath}{\left(\tau_{\rm sync}\right)_{\rm rad}}{^{-1}} = -\frac{1}{\...
...3}}{{(1+q)^{2}}\over{q^2}}{{E_2}^
{-1}}{P^{17/3}\over{R^{7}}}
\end{displaymath} (6)


\begin{displaymath}{\left(\tau_{\rm circ}\right)_{\rm rad}}{^{-1}} = -\frac{1}{e...
...}{{(1+q)^{5/3}}\over{q}}
{{E_2}^{-1}}{P^{7}\over{R^{9}}}\cdot
\end{displaymath} (7)

Studies on tidally driven oscillations did not consider this dissipation until the work by Zahn (1975) who investigated the effects of the radiative damping on dynamical tides. The gravity waves due to the potential of the companion are damped when they reach the external layers. This damping occurs because the radiative cooling time is of the order of the tidal period. Therefore, part of the angular momentum of rotation is transfered to the orbital motion. The associated time scales are characterized by the tidal torque constants En. The complete set of equations we have used to compute En are:

\begin{displaymath}{E_n} =
{3^{8/3}{\left(\Gamma(4/3)\right)^2}\over{(2n+1)\left...
...B\over
{x^{2}}}\right)^
{\prime}_{f}\right]^{-1/3}} {H_{n}^2}
\end{displaymath} (8)

where $\Gamma$ is the gamma function, x is the normalized radius of the configuration, the symbol f denotes the border of the convective core, s indicates surface values, the prime denotes the derivative with respect to x, R is the radius, M the mass, g the gravity, -gB is the square of the Brunt-Väisälä frequency, and B is given by

\begin{displaymath}{B} = {{{\rm d}\over{{\rm d}r}}\ln\rho} - {{1\over{\Gamma_{1}}}{{{\rm d}\over{{\rm d}r}}\ln P}}.
\end{displaymath} (9)

The equation for Hn is given by

 \begin{displaymath}{H_n} = {1\over{X(x_{\rm f})Y(1)}} {\int_{0}^{x_{\rm f}}} {\left[Y^{''} -
{n(n+1)Y\over{x^{2}}}\right]}{X{\rm d}x},
\end{displaymath} (10)

while X is the solution of the differential equation


\begin{displaymath}{X^{''}} - {\rho^{'}\over \rho}{X^{'}} - {n(n+1)\over{x^2}}{X} = 0.
\end{displaymath} (11)

In Eq. (10), Y(1) is the solution of the Clairaut equation. The differential equations were solved for each configuration using a fourth order Runge-Kutta method. We restrict the calculations to n=2 since for larger n the contribution to the dynamic tide is not significant. The above parameters are more complicated to derive than the depth of a convective envelope mentioned before: among other calculations, they depend on the derivative of the Brunt-Väisäla frequency at the boundary of the convective core. For less evolved models, the derivatives are rather simple but as the star evolves, the convective core begins to recede, increasing the chemical composition gradient. The situation is worst near the border of the convective core. In these conditions, additional shell points have to be added by changing the convergence criterions. In some cases up to 6000-7000 mesh points were used.

  \begin{figure}
\par\includegraphics[height=8cm,width=8.8cm,clip]{fig5.eps}
\end{figure} Figure 5: The behavior of the tidal torque constant E2 as a function of the logarithm of the age. The models shown are for 2 and 10 $M_{\odot }$.
Open with DEXTER

The behavior of E2 is illustrated in Fig. 5 for 2 and 10 $M_{\odot }$ stellar models. The influence of the size of the convective core is obvious: the tidal torque constant for the more massive model is about 2 orders of magnitude larger than that for the 2 $M_{\odot }$ model. The tidal torque constant E2 is furthermore more sensitive to stellar evolution than, for example, the apsidal-motion constants. Indeed, from the Main-Sequence up to the last stages of carbon burning, log E2 changes from -5.5 to -16, which is almost 10 orders of magnitude. In the same interval, k2 changes by less than 2 orders of magnitude.

There is also a strong discussion on the validity of the hydrodynamical mechanism proposed by Tassoul (Rieutord 1992; Tassoul & Tassoul 1996, for example). According to this formalism, the distortions due to the proximity effects in close binaries cause hydrodynamical currents which tend to circularize their orbits. A full discussion of this controversy is outside the scope of this paper, but we would like to remark that our tables provide all the necessary material to compute the critical time scales of circularization and synchronization for the hydrodynamical mechanism as well.

By examining Eqs. (2), (3), (6) and (7) we note that the time scales for synchronization and circularization also depend on the moment of inertia, through the radius of gyration $\beta$. The moment of inertia and the potential energy have additional important applications in stellar physics, as for example in investigations of Jacobi dynamics (Ferronsky et al. 1987). The later authors have shown that the product $\alpha _{\rm P}\beta $ ( $\Omega_{\rm P} =
-\alpha_{\rm P} GM^2/R$, where  $\Omega_{\rm P}$ is the potential energy) is almost constant for several laws of mass distribution. However, we have found that this product is not constant for the physical conditions in stellar interiors. In Fig. 6, we show how $\alpha _{\rm P}\beta $ varies with the evolution for 3 selected stellar models: 1, 20 and 125 $M_{\odot }$. The general aspect of Fig. 6 is independent of the mass of the models. As predicted by Ferronsky et al., the product is around 0.4, but only up to core-hydrogen exhaustion, which is indicated by an asterisk. After this evolutionary phase, it increases steeply, breaking the constant tendency.

  \begin{figure}
\par\includegraphics[height=8cm,width=8.8cm,clip]{fig6.eps}
\end{figure} Figure 6: The product $\alpha _{\rm P}\beta $ as a function of the logarithm of the age for models of 1, 20 and 125 $M_{\odot }$. The asterisks indicated the end of hydrogen burning.
Open with DEXTER

5 Rotationally and tidally distorted configurations: The surface brightness distribution

As commented in the Introduction, the rapid rotators and/or components of close binaries are subject to distortions which will alter the brightness surface distribution. The flux depends on the local gravity following the formula by von Zeipel. However, the gravity-darkening exponent $\beta _1$ depends on which mechanism of transport of energy is dominant and on the evolutionary status of the configuration. Therefore, adopting the value 1 for envelopes in radiative equilibrium and 0.32 for envelopes in convective equilibrium is not a fully correct procedure.

The numerical method used to derive $\beta _1$ is based on the triangles strategy (Claret 1998). If the external boundary conditions are unchanged during the evolution, the integrations in the outer layers will be the same, provided that the corresponding triangles in the HR diagram are small. By generalizing this concept to several points, we can simulate distorted configurations with different effective temperature distributions over the surface. By numerically differentiating these envelopes, we finally obtain $\beta _1$. It is important to emphasize that the triangles must be as small as possible in order to guarantee that the external integrations be really unchanged and also to guarantee the validity of the definition $ {\beta_1} = -{\Delta{\log ~ K}\over{\Delta{\log ~g}}}/{{\Delta{\log~ K}\over{\Delta{\log~ T_{\rm eff}}}}}$, where K is the equation of an adiabate.

  \begin{figure}
\par\includegraphics[height=8cm,width=8.8cm,clip]{fig7.eps}
\end{figure} Figure 7: The evolution of the gravity-darkening exponent $\beta _1$ as a function of the effective temperature for a 10 $M_{\odot }$ model. Note how $\beta _1$ decreases/increases during the loops in the HR diagram.
Open with DEXTER

Figure 7 shows how $\beta _1$ depends on the effective temperature for a fixed mass (10 $M_{\odot }$). We selected this model because it crosses the radiative-convective boundary three times, thus presenting envelopes with different physical conditions. For higher effective temperatures the value of $\beta _1$ is 1.0, as predicted by von Zeipel. However, when convection begins to dominate the situation changes drastically and $\beta _1$ decreases rapidly. During the loops in the HR-diagram the gravity-darkening exponent increases-decreases and finally, during the last stages of Carbon burning, it recovers the typical value for convective envelopes (around 0.3). The auxiliary figure within Fig. 7 helps us to connect the chemical-$\beta _1$ evolution.

The first five columns of the tables contain the age (in years), log L (in solar units), log g, log  $T_{\rm eff}$ and the mass (in solar units). Columns 6-9 contain the logarithm of the mass-loss rate ($M_{\odot }$/year), the logarithm of the central temperature, of the density and the size of the convective core $q_{\rm c}$. Columns 10-35 contain the central and surface chemical composition of X, Y, 12C, 13C, 14N, 16O, 17O, 18O, 20Ne, 22Ne, 24Mg, 25Mg, and 26Mg. In the next three columns we have log k2, log k3 and log k4. The depth of the convective envelope $x_{\rm bf}$ is given in Col. 39, while log E2 is given in Col. 40. The constant E, the form factor  $\alpha_{\rm P}$ and $\beta$ are given in positions 41, 42 and 43, and finally, Col. 44 is reserved for the gravity-darkening exponent $\beta _1$. We furthermore provided two kinds of tables: more complete ones with our computed models and more restricted ones containing only equivalent evolutionary points (usually more adequate for the computation of isochrones).

Acknowledgements
I am grateful to B. Willems, who helped to improve the paper. The Spanish MCYT (AYA2003-04651) is gratefully acknowledged for its support during the development of this work.

References

 

Copyright ESO 2004