Issue 
A&A
Volume 546, October 2012



Article Number  A100  
Number of page(s)  12  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201219610  
Published online  11 October 2012 
The new ToulouseGeneva stellar evolution code including radiative accelerations of heavy elements
^{1}
Université de Toulouse, UPSOMP, IRAP,
Toulouse,
France
email: sylvie.theado@ast.obsmip.fr
^{2}
CNRS, IRAP, 14
avenue Édouard Belin, 31400
Toulouse,
France
^{3}
LUTH, Observatoire de Paris, CNRS, Université Paris
Diderot, 5 place Jules
Janssen, 92190
Meudon,
France
email: georges.alecian@obspm.fr
^{4}
Département de Physique et d’Astronomie, Université de
Moncton, Moncton,
NB, E1A 3E9, Canada
email: francis.leblanc@umoncton.ca
Received: 15 May 2012
Accepted: 16 August 2012
Context. Atomic diffusion has been recognized as an important process that has to be considered in any computations of stellar models. In solartype and cooler stars, this process is dominated by gravitational settling, which is now included in most stellar evolution codes. In hotter stars, radiative accelerations compete with gravity and become the dominant ingredient in the diffusion flux for most heavy elements. Introducing radiative accelerations into the computations of stellar models modifies the internal element distribution and may have major consequences on the stellar structure. Coupling these processes with hydrodynamical stellar motions has important consequences that need to be investigated in detail.
Aims. We aim to include the computations of radiative accelerations in a stellar evolution code (here the TGEC code) using a simplified method (SVP) so that it may be coupled with sophisticated macroscopic motions. We also compare the results with those of the Montreal code in specific cases for validation and study the consequences of these coupled processes on accurate models of A and earlytype stars.
Methods. We implemented radiative accelerations computations into the ToulouseGeneva stellar evolution code following the semianalytical prescription proposed by Alecian and LeBlanc. This allows more rapid computations than the full description used in the Montreal code.
Results. We present results for Atype stellar models computed with this updated version of TGEC and compare them with similar published models obtained with the Montreal evolution code. We discuss the consequences for the coupling with macroscopic motions, including thermohaline convection.
Key words: stars: evolution / stars: chemically peculiar / diffusion
© ESO, 2012
1. Introduction
Accurate stellar modeling has recently been given a new boost with the advent of asteroseismology. The observations of oscillating stars and the analysis of the stellar oscillation properties brought new and powerful constraints on stellar models and allowed major progress in our understanding of stellar internal structure. Furthermore, since the discovery of the first extrasolar planets (Wolszczan & Frail 1992; Mayor & Queloz 1995), the spectacular development of the exoplanet research field has also sparked renewed interest in stellar physics, the accurate knowledge of the host star being a necessary condition for characterizing the surrounding planets.
One of the most important successes of astrophysics was the understanding of the basics of stellar internal structure and evolution. This progress was supported by the computations of numerical models. Growing computational resources contributed to refining the description of the physics of the stellar medium (equation of state, opacities, nuclear reactions, etc.) and allowed building a simplified but efficient and widely used “standard model”. However, this standard model does not take the effects of rotation and magnetic fields or the occurrence of accretion or mass loss into account, and it considers convection as the only chemical transport process.
Observations of chemical abundance anomalies in stars and unexpected stellar seismic behaviors have proved the necessity of including “non standard processes” in stellar evolutionary computations. In this framework, the main challenges encountered today by stellar physicists are to better determine the effects of rotation and magnetic fields and understand the chemical transport processes better. This last point, and more specifically modeling of atomic diffusion including the radiative accelerations on individual elements, represents a key ingredient for accurate stellar modeling.
The importance of atomic diffusion inside stars is well established: not only can it modify the atmospheric abundances (Michaud 1970), as observed in the socalled “chemically peculiar stars”, but it can also have strong implications for the stellar internal structure (e.g. Richard et al. 2001). In mainsequence solartype and cooler stars (below about 1.2 M_{⊙}), the radiative accelerations on the heavy elements are generally slower than gravity in absolute value, owing to the small radiation flux compared to hotter stars (Michaud et al. 1976). The elements heavier than hydrogen sink, even if the efficiency of this sinking may be modulated by the radiative accelerations (Turcotte et al. 1998).
One of the biggest success of helioseismology was to prove the importance of atomic diffusion in the Sun. Its introduction in solar evolutionary models has significantly improved the agreement between the sound speed profile inside solar models and that deduced from helioseismic inversion techniques (e.g. ChristensenDalsgaard et al. 1996; Richard et al. 1996; Brun et al. 1998; Turcotte et al. 1998; Schlattl 2002). (This agreement has however been spoiled by the new abundances proposed by Asplund et al. 2005, 2009.) Gravitational settling is now introduced in most stellar evolution codes.
The effects of atomic diffusion on the seismic frequency of mainsequence, solartype star models have been studied by Théado et al. (2005). They have shown that element segregation significantly alters the internal structure of the models and their oscillations frequencies. The frequency differences between models with and without diffusion reach several microHertz for stars with masses greater than 1.3 M_{⊙}.
In hotter stars, the radiative accelerations may become significantly stronger than gravity for many metals, which are pushed up. The variations with depth of the radiative accelerations of specific elements combined with the selective effects of gravitational settling can lead to element accumulation or depletion in various stellar layers. In A and Btype stars, ironpeak element accumulation appear in the Zopacity bump (located at ≃ 200 000 K). The induced opacity increase may lead to local convection (Richer et al. 2000; Richard et al. 2001).
The ironpeak element accumulation in the opacity bump region can help trigger stellar pulsations, therefore improving the agreement between seismic observations and theoretical frequency spectra in many stars: e.g. in γ Doradus, Am, SPB, β Cephei or sdB stars (see Théado et al. 2009, for a detailed discussion of this subject).
The effects of the radiative accelerations on the oscillation frequencies have been tested by Escobar et al. (2012). Very weak effects are observed for models with masses up to 1.28 M_{⊙}, but significant effects are expected for more massive stars.
Progress in the understanding of the physics of earlytype stars is limited by the radiative accelerations not being computed in most stellar evolutionary codes. Up to now, the Montreal stellar evolution code is the only one in which a complete, accurate, and consistent treatment of radiative accelerations has been introduced (Richer et al. 2000; Richard et al. 2001; Turcotte et al. 2000). The price to pay for this accuracy is that the heavy and CPUtime consuming computations of atomic processes do not allow additional sophisticated treatments of macroscopic motions. In the Montreal code, turbulence is only treated as an extremely simplified process with a turbulent diffusion coefficient proportional to a parameterized power of the density.
In this context, we introduced into the ToulouseGeneva evolution code (hereafter TGEC) new computations of the radiative accelerations on heavy elements, following the semianalytical prescription proposed by Alecian & LeBlanc (2002) and LeBlanc & Alecian (2004). This prescription leads to fast but reasonably accurate computations, which represent a good compromise between accuracy and CPUtime consumption and allows coupling with macroscopic motions. Such a treatment of abundance variations inside earlytype stars is necessary for a good understanding of these stars.
In the present paper, we do not introduce thermohaline convection as described in Théado et al. (2009) because we want to present a detailed comparison with the Montreal code in which this specific physical process is not included. Very precise tests of the results obtained by the two codes for iron accumulation inside nearly identical models with similar physics are still underway. In particular, the computations of iron fluxes seem to lead to differences in some cases, which still have to be understood. These tests are beyond the scope of the present paper and will be presented elsewhere in the near future. Here we present a first step in the comparisons, namely the detailed study of the radiative accelerations obtained with the two methods. Some abundance profiles are only shown as indicators of the results that the TGEC code is presently able to obtain.
In the following section, we explain in detail the major improvements implemented in the TGEC code. In Sect. 3, we present a comparison between the Montreal and TGEC computations for two similar stellar models. In Sect. 4, additional models are presented to illustrate the capabilities of the updated TGEC version and its application fields. Our conclusions are given in Sect. 5.
2. New opacities and atomic diffusion computations in TGEC
The ToulouseGeneva stellar evolution code is described in detail in HuiBonHoa (2008); however, the code has recently undergone major improvements that are reported in the following sections. The code can follow the timedependent abundance variations of 21 species (12 elements and their main isotopes: H, ^{3}He, ^{4}He, ^{6}Li, ^{7}Li, ^{9}Be, ^{10}B, ^{12}C, ^{13}C, ^{14}N, ^{15}N, ^{16}O, ^{17}O, ^{18}O, ^{20}Ne, ^{22}Ne, ^{24}Mg, ^{25}Mg, ^{26}Mg, ^{40}Ca, and ^{56}Fe) in detail. The remaining metals are collected into an average species Z.
2.1. Opacities
In the TGEC code, the opacities are computed using the OPCD v3.3 codes and data^{1} (Seaton 2005). They allow computating selfconsistent Rosseland opacities taking the detailed composition of the chemical mixture into account. The opacities are recalculated by considering the abundance variations at each time step and at each mesh point.
2.2. Atomic diffusion
2.2.1. Computational methods
The diffusion computations are based on the Boltzmann equation for dilute collisiondominated plasma. At equilibrium, the solution of the equation is the Maxwellian distribution function. Transport properties in stars are computed considering small deviations from the Maxwellian distribution. In this framework, two different formalism have been proposed to obtain approximate solutions to the Boltzmann equations.
The first method lies on the ChapmanEnskog theory (Chapman & Cowling 1970), which assumes that the total distribution function of a given species can be written as a convergent series, each term of the series representing successive approximations to the distribution function. The formalism is first applied to a test element diffusing in the stellar plasma, taking the diffusion of electrons and the induced electric field into account. The computations lead to a statistical “diffusion velocity” of the testelement with respect to the main component of the plasma. In stars, this basic component is generally hydrogen, so that the diffusion of every element is computed with respect to neutral hydrogen and/or to protons. Meanwhile, the hydrogen abundance is renormalized to recover the abundance consistency and local stellar equilibrium. For the case of helium, which has a nonnegligible abundance, corrections on the diffusion velocity as proposed by Montmerle & Michaud (1976) are added.
The second method has been developed by Burgers (1969). It is based on the Grad 13 moment approximation and the use of a FokkerPlanck collision term in the Boltzmann equation. In this approach, separate flow and heat equations for each component of a multicomponent mixture are solved simultaneously. This method provides a more convenient way for handling multicomponent gases than the ChapmanEnskog one but is heavier to apply. In both methods, the diffusion coefficients can be written as functions of the collisions integrals that depend on the exact nature of the interaction between colliding particules.
In TGEC, the atomic diffusion is computed following the Chapman & Cowling (1970) formalism. We use the gravitational and thermal diffusion coefficients as derived by Paquette et al. (1986), who performed a detailed computation of the collision integrals for ionized elements diffusing in an ionized medium. In the case of the diffusion of neutral atoms in ions, or reverse, the polarization of the neutrals is taken into account as proposed by Michaud et al. (1978). For neutrals atoms diffusing in a neutral medium, the rigid sphere approximation is used. The average diffusing charge Z_{i} of the particules is computed by solving the SahaBoltzmann equations. The partition functions are limited to the statistical weight of fundamental levels, and no correction for high density or degenerate electrons is introduced. The set of SahaBoltzmann equations are solved by using a NewtonRaphson iterative scheme. When the computations of radiative accelerations are added to the diffusion computations, the ChapmanEnskog method is easier to handle than Burgers’. This situation helps treating cases with short time scales and rapid variations, as shown in Sect. 4.
2.2.2. Radiative accelerations
The radiative accelerations on C, N, O, Ne, Mg, Ca, and Fe have been included in the TGEC code following the improved version of LeBlanc & Alecian (2004) of the semianalytical prescription proposed by Alecian & LeBlanc (2002). This method allows very fast computation of radiative accelerations with a reasonable accuracy. Radiative accelerations due to boundbound (Alecian 1985; Alecian & Artru 1990) and boundfree (Alecian 1994) transitions are obtained using a parametric form of the radiative acceleration equation. The basic idea of this parametric method is to derive formula where the terms depending explicitly on atomic data (such as gf values for instance) are separated from those depending on the stellar plasma and the abundances of the considered ion (with the aim of accounting for the saturation effects). In this framework, the radiative accelerations may be approximated by calculating a single value for each parameter found in the related equations. This is the socalled singlevalued parameter (or SVP) approximation (LeBlanc & Alecian 2004).
The interface of the SVP method with the models consists of a set of six parameters per ion, which allows to estimate the radiative acceleration of each element through simple algebraic expressions (LeBlanc & Alecian 2004), for each time step of the run of the evolution code. These parameters are determined only at the beginning of the computation through interpolation as a function of the stellar mass, inside a preestablished grid. The computation of the total radiative acceleration for a given element with SVP also requires computing the ions’ relative populations. This is included in the set of the SVP numerical routines added to TGEC.
The SVP approximation may be implemented in existing codes. There are much less data to process than for detailed radiative acceleration calculations because complete and detailed monochromatic opacities for each element are not needed. A new grid of SVPparameters, well fitted to the stellar mass range considered in this work, was computed following the procedure described in LeBlanc & Alecian (2004). An implementation of the SVP method was first used in Théado et al. (2009).
3. Comparison with the Montreal computations
In this section, we compare the TGEC results to those obtained with the Montreal code. For this purpose we compare the internal structures and the results of the diffusion computations for two models: one computed with the Montreal code, the other computed with TGEC. The Montreal model chosen for this comparison is the 5.3D503, 1.7 M_{⊙}model presented in Richard et al. (2001). A similar model was computed with TGEC, including input physics and initial parameters as close as possible to the Montreal model. As a first step, we present comparisons of seven elements (listed in Table 1), which are especially important for A type stars.
The detailed physics of the Montreal model can be found in Richard et al. (2001), Richer et al. (2000), and Turcotte et al. (1998). The key ingredients are reported in the following sections and compared to those introduced in the TGEC model.
3.1. Input physics
3.1.1. Basic input physics
Initial chemical composition (in mass fraction) for models presented in Sect. 3.
Here is a list that compares the basic input physics used in the TGEC and Montreal models:
Equation of state:
the Montreal model uses the CEFF equation of state (ChristensenDalsgaard & Daeppen 1992). The TGEC model is computed using the OPAL2001 (Rogers & Nayfonov 2002) equation.
Opacities:
in both models, the opacities are computed at every point in the star and at each evolution time step. The Montreal computations take the local abundance of 21 chemical elements into account (see Table 1 of Turcotte et al. 1998) using the OPAL monochromatic data. In TGEC, the opacities are computed as described in Sect. 2.1 considering the detailed abundance of the elements listed in Table 1.
Nuclear reactions:
the Montreal code uses the nuclear energy generation routine of Bahcall & Pinsonneault (1992). The nuclear reaction rates used in TGEC are from the analytical formulae of the NACRE compilation (Angulo 1999).
Convection:
in both Montreal and TGEC models, the delimitation of the convective zones is based on the Schwarzschild criterion. The energy flux in the convection zones is computed following the BöhmVitense mixing length formalism. In the Montreal model, convective zone mixing is modeled as a diffusion process, described by a diffusion coefficient D_{mix} (cf. Sect. 3.1.3). A mixing length approximation is used for D_{mix}, which always leads to high D_{mix} values and very homogeneous convective zones. In TGEC, convective zones are assumed to be instantaneously homogenized. The mixing length parameter α is taken to be equal to 1.68 in both the TGEC and Montreal models.
Initial chemical mixture:
the initial mixture used in the Montreal model is given in Table 1 of Turcotte et al. (1998). The initial composition introduced in the TGEC model is given in Table 1 of this paper: the chemical elements followed in both models are introduced with the same initial abundance.
3.1.2. Atomic diffusion
In the Montreal code, atomic diffusion is computed as described in Richard et al. (2001) and Turcotte et al. (1998). The Montreal calculations take the timedependent variations of 28 species (including isotopes) into account: H, ^{3}He, ^{4}He,^{6}Li, ^{7}Li, ^{9}Be, ^{10}B, ^{11}B, ^{12}C, ^{13}C, N, O, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, K, Ca, Ti, Cr, Mn, Fe, and Ni. Monochromatic OPAL data are used to calculate the Rosseland opacities and the radiative accelerations at each evolution time step. The abundances are updated at every iteration over the stellar structure. The diffusion coefficients and velocities are determined by solving the Burgers’s flow equations for ionized gases (Burgers 1969) for all diffusing elements. The collisions integrals are from Paquette et al. (1986).
As a first step, we include, in the TGEC model, the atomic diffusion (including radiative accelerations) of seven elements (eight species including isotopes), especially important for A type stars: H, ^{3}He, ^{4}He, ^{12}C, ^{14}N, ^{16}O, ^{40}Ca, and ^{56}Fe. The atomic diffusion is computed as described in Sect. 2.2. The diffusion velocities are computed following the Chapman & Cowling (1970) formalism in the testatom approximation, and the diffusion coefficients are from Paquette et al. (1986).
3.1.3. Macroscopic transport of chemical elements
Macroscopic transport in the Montreal code In the Montreal code, macroscopic transport processes, including (semi)convection and turbulent mixing, are modeled as diffusion processes by adding a pure diffusion term to each element’s diffusion velocity V_{pi}: (1)where D_{T} and D_{mix} are turbulent diffusion coefficients with D_{mix} representing the effects of convective and semiconvective motions and D_{T} parametrizing turbulent transport. The variable X_{i} represents the mass fraction of species i.
In radiative layers D_{mix} = 0, while in convective zones it is computed using the mixing length theory (as stated in Sect. 3.1.1). In semiconvection zones, D_{mix} is assumed proportional to (∇−∇_{ad})/(∇_{L} − ∇), where ∇ and ∇_{ad} stands for the local and adiabatic logarithmic temperature gradients and ∇_{L} is a function of ∇_{μ} = dlnμ/dlnP. We refer the reader to Sect. 2.2 of Richard et al. (2001) for a more detailed description of the treatment of convection.
In the 5.3D503model, the turbulent transport is chosen strong enough to guarantee, during the whole evolution period considered, a complete mixing throughout the region between the surface and the point where log T_{0} = 5.3 (see Fig. 2 of Richard et al. 2001). This mixing mimics the effects of a Fe convection zone except that it is imposed throughout evolution. Below log T_{0} = 5.3, the diffusion coefficient D_{T} obeys the following algebraic dependence on density (Eq. (1) of Richer et al. 2000): (2)where ρ_{0} = ρ(T_{0}), ω = 50, n = 3 and (3)Macroscopic transport in the TGEC code In the TGEC code, the macroscopic motions of a given species i are parametrized, as in the Montreal code, by including a turbulent diffusion term to the diffusion velocity of the considered species: (4)The TGEC code does not include a D_{mix} term, since convection is treated as an instantaneous dilution. The diffusion coefficient D_{T} is chosen to mimic the mixing (D_{T} + D_{mix}) introduced in the 5.3D503 Montreal model. Figure 1 represents this diffusion coefficient inside our model, it may be compared to Fig. 2 of Richard et al. (2001). From the surface down to log T = 5.3, D_{T} is chosen to homogenize the stellar material, below this region D_{turb} is computed using the same expression as introduced in the Montreal code (i.e. following Eq. (2)).
Fig. 1 Turbulent diffusion coefficient in the 170 Myr, 1.7 M_{⊙}model computed with TGEC. D_{turb} is constant from the surface down to log T = 5.3. Just below this region, D_{turb} is by definition 50 times larger than the He atomic diffusion coefficient (see Eq. (2)). The vertical dashed line locates the region where log T = 5.3. 
3.2. Results
The Montreal and TGEC codes compute stellar evolution from premain sequence up to the subgiant branch. They were assumed to be homogeneous on the premain sequence, atomic diffusion, and turbulent transport were introduced at the beginning of the mainsequence phase.
Figure 2 compares the evolutionary tracks of the two models (for clarity the premainsequence evolution is not shown). In spite of different treatments of some physical ingredients (e.g. equation of state, opacities, nuclear reactions, diffusion), the two tracks appear relatively close to each other.
Fig. 2 Evolutionary tracks of two 1.7 M_{⊙} sequences computed with TGEC or with the Montreal code. The dashed curve represents the 5.3D503 model presented in Richard et al. (2001) (and computed with the Montreal code). The solid curve represents a model computed with TGEC and including similar input physic. The black dots and crosses respectively represent the 30 and 400 Myrmodels. 
We first give a detailed comparison between the Montreal and TGEC calculations for the two 30 Myr models. Figure 3 compares the internal structure of these two models. The relative differences in pressure, temperature, and density never exceed 6% (and even less for P and ρ). The maximum relative differences are observed near log (ΔM/M) ≃ −6.3, which corresponds to the base of the potential iron convective zone.
Similarly to Fig. 7 of Richard et al. (2001), Fig. 4 displays the radiative accelerations and the local gravity of the two models. The discrepancy between g_{MTR} and g_{TGEC} remains smaller than around 0.3 dex except for O in the upper layers (in the convection zone). We notice that the average accuracy of the SVP method is estimated to ± 0.1 dex (Alecian & LeBlanc 2002) with respect to the detailed computations of Seaton (1997). However, the SVP approximation is less accurate for light elements (lighter than O) than for heavy ones because this method uses parameters determined at the position of the maximum of each ion relative population. Because there are fewer ionization states for light elements, the gap in temperature between these maxima is larger for them. The comparison of Fig. 4 is fairly satisfactory, since radiative accelerations are computed in completely different ways and use different atomic and opacity databases.
Fig. 3 Comparison of the internal structures (pressure, temperature and density) of the two 1.7 M_{⊙}models represented in Fig. 2. The comparison is shown at 30 Myr. MTR stands for the model computed with the Montreal stellar evolution code while TGEC stands for the model computed with the ToulouseGeneva stellar evolution code. 
Fig. 4 Radiative accelerations (in cgs unit) at 30 Myr in the TGEC and Montrealmodels represented in Fig. 2. The solid curves represents the accelerations in the TGEC model, while the dashed curve represent those in the Montreal model. The longdashed curve represents the local gravity in the TGEC model. (The local gravity in the Montreal model is not represented as it is indistinguishable from that of the TGEC model.) 
Figure 5 presents the diffusion velocities of several elements for the same models. Despite the different atomic diffusion prescriptions used, the diffusion velocities are quite close in most of the star.
Fig. 5 Diffusion velocities (in cgs units) below the surface convective zone in the Montreal (dashed curves) and TGEC (solid curves) 30 Myrmodels represented in Fig. 4. The vertical lines locate the bottom of the completely mixed region. The peaks observed for Ca and Fe correspond to sign changes in the diffusion velocities: the Ca diffusion velocity is positive for log (ΔM/M) ≲ −6 and −3.5 ≲ log (ΔM/M) ≲ −2.5. The Fe diffusion velocity is positive for −7 ≲ log (ΔM/M) ≲ −4.5. 
Figure 6 shows the abundance profiles in the two 30 and 400 Myrmodels. At 30 Myr the profiles are quite close for all the species followed in TGEC. However, at 400 Myr discrepancies in the abundance profiles are observed. These differences are small for most elements (He, C, N, O, and Fe) except for Ca. As shown in Fig. 2, the positions of the two 400 Myr models in the HR diagram are slightly different, and are in particular more distant from each other than the 30 Myr models. As a consequence, their internal structure is also expected to be different. In this context, it seems difficult to disentangle the discrepancies in the models due to structure variations from those caused specifically by the difference in diffusion calculations.
Fig. 6 Abundance profiles in the TGEC (solid curves) and Montreal (dashed curves) models shown in Fig. 2 as a function of the outer mass fraction. log (X/X_{0}) represents the ratio between the current element mass fraction and its initial values (at t = 0). The blue and orange curves present the abundances at 30 Myr and 400 Myr after the ZAMS, respectively. 
4. Application field of the TGEC code
In this section, we propose to demonstrate the abilities of the TGEC code and its application field. Models of AF stars including atomic diffusion and minimal mixing (i.e. only a mild mixing below the convective zones to avoid steep and unrealistic composition gradients at the transition between radiative and convective regions) have already been presented in Théado et al. (2009). These models, which were evolved until the end of the mainsequence phase, have shown the ability of the code to compute complete evolutionary tracks including radiative levitation effects in the presence of minimal mixing.
In this paper, we propose another test of the abilities of the code to manage rapid variations in the chemical composition and the thermal structure. For this purpose we present a set of models computed with TGEC with masses ranging from 1.5 to 2.5 M_{⊙}. Like the models presented in Théado et al. (2009), they include atomic diffusion and a mild mixing below the convective zones. In the present models, the iron convective zone, when it appears, is assumed to be connected and completely mixed with the H/He convective envelope. These models were evolved from premain sequence up to hydrogencore exhaustion. They were assumed to be chemically homogeneous on the premain sequence. Atomic diffusion was introduced at the beginning of the mainsequence phase. Because the present computations are done to test the TGEC code, we did not include here the thermohaline convection that must be added for comparison with real stars (Théado et al. 2009).
4.1. Input physics
As previously described, our models were computed using the OPAL2001 (Rogers & Nayfonov 2002) equation of state. The nuclear reactions were from the NACRE compilation (Angulo 1999), and the opacities were computed as described in Sect. 2.1. Atomic diffusion (including the radiative forces) was introduced as described in Sect. 2.2, for the following elements: H, ^{3}He, ^{4}He, ^{12}C, ^{14}N,^{16} O, ^{40}Ca, and ^{56} Fe. The initial metal mixture was the solar mixture of Grevesse & Noels (1993). The initial mass fractions for the diffusing species are given in Table 3. Convection was computed using the mixing length theory with a mixing length parameter α = 1.8. The HI and HeII convective zones were assumed to be connected by overshooting and mixed together. As stated previously, the iron convective region that may appear during mainsequence evolution was also assumed connected and completely mixed with the surface convective region.
To avoid the appearance of sharp and nonphysical abundance gradients at the transition between radiative and convective regions, we introduced mild mixing below the surface convective zone. This mixing was modeled as a diffusion process (Schatzman 1969) as described in Eq. (4) with (5)where D_{czb} and r_{czb} are the value of D_{T} and the value of the radius at the base of the convective zone, respectively and Δ is the half width of the mixing region, D_{czb} and Δ are free parameters chosen to produce a mild mixing to a small extent. The value of D_{czb} is taken as equal to 2 × 10^{5} cm^{2} s^{1}, and Δ is fixed to 0.5% of the radius below the surface convective region. Figure 7 presents this diffusion coefficient below the surface convective zone of 1.7 M_{⊙} models at two evolutionary stages. At 170 Myr, the surface convective region includes the H and He convective zones (there is no Fe convective zone), and at 588 Myr the external convective zone is composed of the H, He, and Fe convective zones assumed connected and mixed.
Characteristics of the models.
The physics introduced in these models is close to that of the r303M, 1.5 M_{⊙}model presented in Richard et al. (2001). Like our models, the r303M model includes mild mixing below the He convective zone. The turbulent diffusion coefficient is chosen as large as 10^{4} cm s^{1} at the base of the convective zone. It then rapidly decreases with depth, varying like ρ^{3}. Turbulent transport is chosen to be large enough to guarantee complete mixing between the Fe and He convective zones whenever the Fe convective zone appears. The initial metal mixture and the mixing length parameter are, however, different in the Montreal and TGEC models. For CPUtime consuming reasons, the evolution of the Montreal model was stopped at 89 Myr, and the results for a 1.5 M_{⊙} model are the only ones presented. TGEC models with masses up to 2.5 M_{⊙} are computed along the whole mainsequence phase. Since the physics of our models and the r303M Montreal model are close but not similar, no detailed comparison between them can be carried out. However, qualitative similarities are underlined.
4.2. Results
Figure 8 displays the evolutionary tracks of our models. (For clarity premainsequence evolution is not represented.) In the following, we describe, as a representative example, the results obtained for a 2.1 M_{⊙} model. We focus on irondiffusion related features.
Initial chemical composition (in mass fraction) for models presented in Sect. 4.
Fig. 7 Turbulent diffusion coefficient (according to Eq. (5)) in a 1.7 M_{⊙}model computed with TGEC at two evolutionary stages. The solid line represents the diffusion coefficient below the surface convective zone (due to the ionization of H and He) at 170 Myr, the dashed line represents the diffusion coefficient below the iron convective zone (assumed connected with the H and He ones) at 588 Myr, the vertical dotted line represents the base of the iron convective zone. 
Fig. 8 Evolutionary tracks of 1.5 to 2.5 M_{⊙} models computed with atomic diffusion (including radiative acceleration effects) and a mild mixing below the convective enveloppe. 
Figure 9 presents the profiles of various quantities inside the 2.1 M_{⊙} model along the mainsequence evolution. The left and middle columns illustrate early mainsequence evolution, the right hand column displays later evolutionary stages.
Fig. 9 Internal structure of the 2.1 M_{⊙}model represented in Fig. 8. The profiles are represented as a function of the outer mass fraction and at various evolutionary steps. Upper panel: iron abundance profiles, ^{56}Fe/^{56}Fe_{0} represents the ratio between the current Fe abundance and its initial value (cf. Table 3). The various curves are defined in the upper panels and are valid for all the plots of the respective column. Middle panel: Rosseland opacity profiles. Lower panel: difference between the radiative and adiabatic gradients. 
Fig. 10 Variation with time of the bottom of the surface convective zone in the 2.1 M_{⊙} model represented in Fig. 8. Left panel: zoom on the 0–0.12 Gyr period. Right panel: zoom on the 0.12–0.8 Gyr period. The dashed vertical lines are located respectively at 25.9 Myr, 32.6 Myr, 117.1 Myr, 254.3 Myr, 380.3 Myr, 576.0 Myr, and 726.0 Myr. 
On the main sequence, the competition between the gravitational settling and the radiative acceleration leads to iron accumulations in the outer regions of the model. We note iron enrichment in the surface H/He convective zone and in the Zbump region (at around log (ΔM/M) = −7 (T ≃ 200 000 K). In this region, iron enrichment significantly increases the opacity, hence the radiative gradient, and when this radiative gradient exceeds the adiabatic gradient, the socalled iron convective zone appears. Here it is assumed to be connected and mixed up to the H/He surface convective zone through overshooting. As a result, the bottom of the surface convective zone sinks abruptly to log (ΔM/M) ≃ −7. The dilution in this thick convective zone decreases the iron abundance in the opacity bump: the opacity and the radiative gradient consequently decrease, which leads to the disappearance of the iron convective zone. The surface convective zone then recedes. A new cycle starts as radiative forces proceed in accumulating iron in the opacity bump region. The middle column of Fig. 9 shows a second iron accumulation/deep convection phase.
The alternation between convective and radiative episodes in the Zbump proceeds from ≃ 20 to ≃ 115 Myr. A thick convective zone then stays for the subsequent mainsequence evolution. The right hand column of Fig. 9 illustrates the internal structure variations of our model during the later evolutionary phases.
Figure 10 shows the position of the bottom of the surface convective zone of our model during the main sequence. At the beginning of the mainsequence phase, the convective zone includes the H and He ionization regions. After a few million years and because of the He gravitational settling, the HeII convective zone disappears leading to an abrupt recession of the convective region (whose bottom subsequently lies at log (ΔM/M) ≃ −11.25). After this first receding episode, the convective region undergoes a 100 Myrperiod of rapid variations during which the convective depth oscillates between log (ΔM/M) ≃ −11.2 and log (ΔM/M) ≃ −6.5. After this period, a thick surface convective zone appears that slowly deepens inside the interior during the rest of the mainsequence phase.
The lasting thick convective zone results from two features explained below.

Because of structural variations forced by the star’s evolution, theiron accumulation needed to induce convection in the opacitybump varies along the main sequence. After115 Myr, in particular, smaller ironaccumulations in the Zbump are needed to initiate the ironconvective zone.

The succession of convective sinking and receding episodes, which occurs from 20 to ≃ 115 Myr, leads to rapid variations in the surface iron abundance. As an illustration of this phenomenom, Fig. 11 displays the iron abundance evolution over 30 Myr (including 5 convective sinking/receding episodes). The global effect of these episodes is illustrated in Fig. 12, which presents the timedependent variations of the averaged iron surface abundance (averaged over 20 Myr). The diffusioninduced iron enrichments (in the Zbump but also in the H/He convective zone) combined with the deep convective mixing episodes leads on average to an iron enrichment of the surface convective zone. For high values of this enrichment, the convective dilution in the presence of an iron convective zone may become unable to efficiently decrease the iron abundance in the Zbump and to suppress the iron convective zone. A thick convective zone then persists during the subsequent evolution.
According to Fig. 12, the iron abundance in the surface convective zone continues to increase after 115 Myr (because of radiative accelerations), while the convective zone slowly deepens. After 400 Myr, the bottom of the convective zone reaches ironunderabundant regions, which causes as a slow decrease in the iron surface abundance.
Fig. 11 Time dependent variations in the Fe surface abundance in the 2.1 M_{⊙} model presented in Fig. 8. The iron abundance is shown on a 30 Myr period including 5 convective sinking/receding episodes. This graph is similar to Fig. 16 of Richard et al. (2001) (see text for more details). 
Figure 13 presents the convective zone extension (vs time) in all computed models. For an easy comparison with other masses, the 2.1 M_{⊙} model is shown in this figure. A succession of convective sinking and receding episodes is observed in all the computed models. This period of rapid oscillations of the convective depth is shorter for higher stellar masses. In the 1.5 M_{⊙} model, it persists during most of the mainsequence phase, whereas in the 2.5 M_{⊙} model, a thick convective zone appears in less than 80 Myr.
Fig. 12 Iron surface abundance along mainsequence evolution for the 2.1 M_{⊙} model represented in Fig. 8. During the 0–120 Myr period (i.e. when the external convective zone undergoes rapid and large variations), the iron abundance shown is an averaged value (over 20 Myr). 
The patterns observed during the early mainsequence phase of our models are consistent with the results presented by Richard et al. (2001) for their r303M, 1.5 M_{⊙} model. The alternation between radiative and convective episodes in the opacity bump that occurs in our models is similarly observed in the Montreal model as described in Sect. 3 (paragraph “Model with mixing episodes”) of Richard et al. (2001). The comparison between our Figs. 11 and 16 of Richard et al. (2001) (middle panel) also shows similar behaviour to the surface iron abundances. But while the Montreal code is limited by time consuming issues (their r303M, 1.5 M_{⊙} model is evolved for less than 100 Myr), TGEC is able to manage the rapid structural and compositional variations due to the succession of convective sinking and receding episodes and allows computing complete mainsequence evolutionary tracks for masses up to at least 2.5 M_{⊙}.
5. Conclusion
The implementation of radiative diffusion effects in stellar modeling turns out as a necessary condition for any computations of accurate models of F, A, and B type stars. Different ways of introducing these effects are considered by theoreticians, but few reliable methods are presently known.
The most accurate method consists in computing the timedependent variations of the chemical species present in the stellar mixture while using a precise description of the diffusion process. This includes an accurate determination of the radiative accelerations applied to each species and the resolution of the Burgers equations to obtain the diffusion velocities. Such a method has been successfully introduced in the Montreal code, which has been used for many years for the modeling of mainsequence FAB type stars. This accurate method allows precise determinations of abundance variations due to atomic diffusion, but it leads to highly timeconsuming codes, which makes the computations of complete mainsequence evolutionary tracks difficult.
The development of asteroseismology during the past 30 years and, in particular, the discovery of a wide variety of pulsating F, A, and Btype stars have brought to light the necessity for developing stellar modeling tools able to provide accurate models adapted to seimic analysis for these stars. The seismic modeling of a target star is an iterative process that requires computating numerous models before getting the best one. In this context, using a code like the Montreal one appears particularly unadapted.
Fig. 13 Variation in the surface convective zone extension along the mainsequence phase for the same models as presented in Fig. 8. The vertical axis represents the outer mass fraction, the solid curves represent the bottom of the external convective zone. 
This encouraged us to develop an alternative code that includes a more flexible prescription of the radiative accelerations. In this context we included the radiative diffusion in the TGEC code using the semianalytical prescription proposed by Alecian & LeBlanc (2002) and LeBlanc & Alecian (2004) (called the SVP method). This method, less accurate than the one introduced in the Montreal code, has been proven to give good results with reasonable precision. With the most recent version of TGEC, this method has been successfully introduced for the first time in a stellar evolution code. In TGEC, the diffusion velocities are then computed following the Chapman and Cowling formalism in the testatom approximation. The timedependent variations in the chemical species are computed consistently.
To validate our code, we did a comparison between a model computed with the Montreal code and a similar model computed with TGEC. (The two models include input physics that are as close as possible.) The comparison shows good agreement in the computations of the radiative acceleration for models with similar internal structures.
Nickel is not presently included in the TGEC computations, as its atomic data are not included in TOPBase (Cunto et al. 1993), the atomic data base of the Opacity Project (Seaton et al. 1992) from which the atomic data used in the SVP approximation are taken. This explains its absence in the SVP package. Its introduction in TOPBase^{2}, as well as in the SVP package, is underway.
To demonstrate both the abilities of the TGEC code and its flexibility as compared to the Montreal one, we presented sets of models that include atomic diffusion with small additional mixing. In this context, models including “minimal mixing” were already presented in Théado et al. (2009); in Sect. 4 we chose to present a new set of models (with masses ranging from 1.5 to 2.5 M_{⊙}) including physical ingredients close to the r303M, 1.5 M_{⊙} model of Richard et al. (2001). In these models, only mild mixing is introduced below the surface convective zone, and the iron convection region, when its appears, is assumed to be connected and mixed with the surface convective envelope. Under these assumptions, our models show similar behavior to the Montreal models, at least during the first 100 Myr (where the Montreal computations stop). In particular, the introduced physical hypothesis leads to an alternation of convective and radiative episodes in the opacity bump region in both the Montreal and the TGEC models. We have also shown that the TGEC code allows computing complete evolutionary tracks (up to the subgiant branch) for AF stars with masses up to 2.5 M_{⊙}. It is able to manage rapid variations in the chemical composition, under cover of mild mixing at the transition between convective and radiative regions. The effects of quasipure atomic diffusion can then be evaluated in the context of stellar evolutionary models.
As discussed in the introduction, precise tests of the results obtained by the two codes for the case of iron accumulation inside stars, when thermohaline convection is not taken into account, still have to be performed, since differences in the maximum iron value appear in some cases. Such tests are underway and will be presented elsewhere in the near future.
The improvements brought to the atomic diffusion aspects of the TGEC code presented here and in Théado et al. (2009) makes it an excellent test bed for stellar evolution and asteroseismology. Its application to various types of stars will certainly help us understand them better.
The OPCD_3.3 packages is available on the following website: http://cdsweb.ustrasbg.fr/topbase/op.html
Acknowledgments
We thank Olivier Richard for kindly providing its 5.3D503 1.70 M_{⊙} model. We acknowledge the financial support of the Programme National de Physique Stellaire (PNPS) of CNRS/INSU, France and the Natural Sciences and Engineering Research Council of Canada.
References
 Alecian, G. 1985, A&A, 145, 275 [NASA ADS] [Google Scholar]
 Alecian, G. 1994, A&A, 289, 885 [NASA ADS] [Google Scholar]
 Alecian, G., & Artru, M. 1990, A&A, 234, 323 [NASA ADS] [Google Scholar]
 Alecian, G., & LeBlanc, F. 2002, MNRAS, 332, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Angulo, C. 1999, in AIP Conf. Ser., 495, 365 [Google Scholar]
 Bahcall, J. N., & Pinsonneault, M. H. 1992, Rev. Modern Phys., 64, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., TurckChièze, S., & Morel, P. 1998, ApJ, 506, 913 [NASA ADS] [CrossRef] [Google Scholar]
 Burgers, J. M. 1969, Flow Equations for Composite Gases (New York: Academic Press) [Google Scholar]
 Chapman, S., & Cowling, T. G. 1970, The mathematical theory of nonuniform gases. an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases (Cambridge University Press) [Google Scholar]
 ChristensenDalsgaard, J., & Daeppen, W. 1992, A&ARv, 4, 267 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [CrossRef] [Google Scholar]
 Cunto, W., Mendoza, C., Ochsenbein, F., & Zeippen, C. J. 1993, A&A, 275, L5 [NASA ADS] [Google Scholar]
 Escobar, M., Théado, S., Vauclair, S., et al. 2012, A&A, 543, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse, 15 [Google Scholar]
 HuiBonHoa, A. 2008, Ap&SS, 316, 55 [NASA ADS] [CrossRef] [Google Scholar]
 LeBlanc, F., & Alecian, G. 2004, MNRAS, 352, 1329 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Michaud, G. 1970, ApJ, 160, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Michaud, G., Charland, Y., Vauclair, S., & Vauclair, G. 1976, ApJ, 210, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Michaud, G., Martel, A., & Ratel, A. 1978, ApJ, 226, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Montmerle, T., & Michaud, G. 1976, ApJS, 31, 489 [NASA ADS] [CrossRef] [Google Scholar]
 Paquette, C., Pelletier, C., Fontaine, G., & Michaud, G. 1986, ApJS, 61, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, O., Michaud, G., & Richer, J. 2001, ApJ, 558, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, O., Vauclair, S., Charbonnel, C., & Dziembowski, W. A. 1996, A&A, 312, 1000 [NASA ADS] [Google Scholar]
 Richer, J., Michaud, G., & Turcotte, S. 2000, ApJ, 529, 338 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
 Schatzman, E. 1969, A&A, 3, 331 [NASA ADS] [Google Scholar]
 Schlattl, H. 2002, A&A, 395, 85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seaton, M. J. 1997, MNRAS, 289, 700 [NASA ADS] [CrossRef] [Google Scholar]
 Seaton, M. J. 2005, MNRAS, 362, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Seaton, M. J., Zeippen, C. J., Tully, J. A., et al. 1992, Rev. Mex. Astron. Astrofis., 23, 19 [NASA ADS] [EDP Sciences] [Google Scholar]
 Théado, S., Vauclair, S., Castro, M., Charpinet, S., & Dolez, N. 2005, A&A, 437, 553 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Théado, S., Vauclair, S., Alecian, G., & LeBlanc, F. 2009, ApJ, 704, 1262 [NASA ADS] [CrossRef] [Google Scholar]
 Turcotte, S., Richer, J., Michaud, G., Iglesias, C. A., & Rogers, F. J. 1998, ApJ, 504, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Turcotte, S., Richer, J., Michaud, G., & ChristensenDalsgaard, J. 2000, A&A, 360, 603 [NASA ADS] [Google Scholar]
 Wolszczan, A., & Frail, D. A. 1992, Nature, 355, 145 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Turbulent diffusion coefficient in the 170 Myr, 1.7 M_{⊙}model computed with TGEC. D_{turb} is constant from the surface down to log T = 5.3. Just below this region, D_{turb} is by definition 50 times larger than the He atomic diffusion coefficient (see Eq. (2)). The vertical dashed line locates the region where log T = 5.3. 

In the text 
Fig. 2 Evolutionary tracks of two 1.7 M_{⊙} sequences computed with TGEC or with the Montreal code. The dashed curve represents the 5.3D503 model presented in Richard et al. (2001) (and computed with the Montreal code). The solid curve represents a model computed with TGEC and including similar input physic. The black dots and crosses respectively represent the 30 and 400 Myrmodels. 

In the text 
Fig. 3 Comparison of the internal structures (pressure, temperature and density) of the two 1.7 M_{⊙}models represented in Fig. 2. The comparison is shown at 30 Myr. MTR stands for the model computed with the Montreal stellar evolution code while TGEC stands for the model computed with the ToulouseGeneva stellar evolution code. 

In the text 
Fig. 4 Radiative accelerations (in cgs unit) at 30 Myr in the TGEC and Montrealmodels represented in Fig. 2. The solid curves represents the accelerations in the TGEC model, while the dashed curve represent those in the Montreal model. The longdashed curve represents the local gravity in the TGEC model. (The local gravity in the Montreal model is not represented as it is indistinguishable from that of the TGEC model.) 

In the text 
Fig. 5 Diffusion velocities (in cgs units) below the surface convective zone in the Montreal (dashed curves) and TGEC (solid curves) 30 Myrmodels represented in Fig. 4. The vertical lines locate the bottom of the completely mixed region. The peaks observed for Ca and Fe correspond to sign changes in the diffusion velocities: the Ca diffusion velocity is positive for log (ΔM/M) ≲ −6 and −3.5 ≲ log (ΔM/M) ≲ −2.5. The Fe diffusion velocity is positive for −7 ≲ log (ΔM/M) ≲ −4.5. 

In the text 
Fig. 6 Abundance profiles in the TGEC (solid curves) and Montreal (dashed curves) models shown in Fig. 2 as a function of the outer mass fraction. log (X/X_{0}) represents the ratio between the current element mass fraction and its initial values (at t = 0). The blue and orange curves present the abundances at 30 Myr and 400 Myr after the ZAMS, respectively. 

In the text 
Fig. 7 Turbulent diffusion coefficient (according to Eq. (5)) in a 1.7 M_{⊙}model computed with TGEC at two evolutionary stages. The solid line represents the diffusion coefficient below the surface convective zone (due to the ionization of H and He) at 170 Myr, the dashed line represents the diffusion coefficient below the iron convective zone (assumed connected with the H and He ones) at 588 Myr, the vertical dotted line represents the base of the iron convective zone. 

In the text 
Fig. 8 Evolutionary tracks of 1.5 to 2.5 M_{⊙} models computed with atomic diffusion (including radiative acceleration effects) and a mild mixing below the convective enveloppe. 

In the text 
Fig. 9 Internal structure of the 2.1 M_{⊙}model represented in Fig. 8. The profiles are represented as a function of the outer mass fraction and at various evolutionary steps. Upper panel: iron abundance profiles, ^{56}Fe/^{56}Fe_{0} represents the ratio between the current Fe abundance and its initial value (cf. Table 3). The various curves are defined in the upper panels and are valid for all the plots of the respective column. Middle panel: Rosseland opacity profiles. Lower panel: difference between the radiative and adiabatic gradients. 

In the text 
Fig. 10 Variation with time of the bottom of the surface convective zone in the 2.1 M_{⊙} model represented in Fig. 8. Left panel: zoom on the 0–0.12 Gyr period. Right panel: zoom on the 0.12–0.8 Gyr period. The dashed vertical lines are located respectively at 25.9 Myr, 32.6 Myr, 117.1 Myr, 254.3 Myr, 380.3 Myr, 576.0 Myr, and 726.0 Myr. 

In the text 
Fig. 11 Time dependent variations in the Fe surface abundance in the 2.1 M_{⊙} model presented in Fig. 8. The iron abundance is shown on a 30 Myr period including 5 convective sinking/receding episodes. This graph is similar to Fig. 16 of Richard et al. (2001) (see text for more details). 

In the text 
Fig. 12 Iron surface abundance along mainsequence evolution for the 2.1 M_{⊙} model represented in Fig. 8. During the 0–120 Myr period (i.e. when the external convective zone undergoes rapid and large variations), the iron abundance shown is an averaged value (over 20 Myr). 

In the text 
Fig. 13 Variation in the surface convective zone extension along the mainsequence phase for the same models as presented in Fig. 8. The vertical axis represents the outer mass fraction, the solid curves represent the bottom of the external convective zone. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.