Table Of ContentAstrophysics and Space Science (CoRoT/ESTA Volume) manuscript No.
(will be inserted by the editor)
P. Morel · Y. Lebreton
CESAM: a free code for stellar evolution calculations
8
0
0
2
n
a
Received date / accepted date
J
4
1 Abstract The cesam code is a consistent set of pro- 1 Introduction to CESAM.
grams and routines which perform calculations of 1D
]
h quasi-hydrostaticstellar evolutionincluding microscopic Within the limitations due to electronic degeneracy,ce-
p diffusionofchemicalspeciesanddiffusionofangularmo- samallowsthecomputationofthequasi-staticevolution
o- mentum. The solution of the quasi-static equilibrium is ofstellarmodelsaslongastheassumptionofquasi-static
r performed by a collocation method based on piecewise equilibrium remains valid, that is to say, until the ex-
st polynomials approximationsprojected ona B-spline ba- haustion of oxygen in the core. The modular structure
a sis; that allows stable and robust calculations, and the of cesam facilitates the choices among several physical
[ exactrestitutionof the solution,notonly at gridpoints, formalisms,forequationofstate(hereafterEOS),convec-
1 even for the discontinuous variables. Other advantages tion,opacities,diffusion coefficients,etc... Many nuclear
v are the monitoring by only one parameter of the accu- networksandinitialmixturesareavailablewhichallowto
9 racy and its improvement by super-convergence.An au- optimise the physical description according to the kind
1 tomaticmeshrefinementhasbeendesignedforadjusting of model and evolutionary phase of interest. Mass loss
0 the localisations of grid points according to the changes and infall of planetoids are also implemented.
2 of unknowns. For standard models, the evolution of the
.
1 chemical composition is solved by stiffly stable schemes
0 of orders up to four; in the convectionzones mixing and
The available packages. Earlier versions, cesam2-3-4,
8 evolution of chemical are simultaneous. The solution of
were programmed in F77 and cesam5 in F90. Though
0
the diffusion equation employs the Galerkin finite ele-
: obsolete, cesam4 and cesam5 are available at:
v mentsscheme;themixingofchemicalsisthenperformed
http://www.obs-nice.fr/morel/CESAM
i by a strong turbulent diffusion. A precise restoration of
X In the early 2000’s, cesam was re-programmed in F95
the atmosphere is allowed for.
and named cesam2k. Three versions are now available:
r
a
– The“fixed version”andthe“fixed COROT version”,
Keywords Methods: numerical sun: evolution sun: limited to 3α burning, available respectively at:
interior stars: evolution stars: i·nterior · http://www.obs-nice.fr/cesam/
· · http://perso.obspm.fr/∼lebreton/Modeles/CESAM.html
– The“β version”morecomplete,notfixed,stillinde-
velopment(not free of bugs)is hereafter signalledby
PACS 97.10.Cv 97.10.Sj 95.75.Pq
· · the flag(β). Itincludes evolutionuptooxygenburn-
ing, diffusion of the angular momentum and other
developments ofminor importance.It is availableat:
P. Morel http://www.obs-nice.fr/morel/CESAM
Cassiop´ee,URACNRS1362,ObservatoiredelaCˆoted’Azur,
Nice, France Each package contains five directories and two files:
Tel.: +033-492003039
– “SOURCE”,contains the fortran sources.
Fax: +033-492003121
E-mail: [email protected] – “EXPLOIT”,contains programs to exploit the models
and examples of input files.
Y.Lebreton
Observatoire de Paris, GEPI, CNRS UMR 8111, 5 Place – “SUN STAR DATA”, contains physical data and pro-
Janssen, 92195 Meudon, France grams for their implementation.
E-mail: [email protected] – “TESTS”,containsprogramsperformingvariouschecks.
2 P.Morel, Y.Lebreton
– “SCRIPTS”, contains scripts to be used for the im- Theotherinputfileshaveonlyspecificfunctions,among
plantation and operating, and a MAKEFILE. them the most useful are:
– “aide mem2k.ps”,ashortguideofdirectionsforuse,
– “mixture” allows the use of an initial mixture not
the “aide-m´emoire” (hereafter Paper 2).
implemented in cesam.
– “cesam2k.ps”,a complete description of the numer-
– “modif mix”allowstomodifytheabundancesofsome
ical aspects and physics implemented, the “notice”
species in a mixture already implemented.
(hereafter Paper 3).
– “reglages” allows to personalise the kind of preci-
The source is structured in 14 modules: numerical rou-
sion to be used (see below).
tines, opacities, convection, etc. All the F95 routines of
– “planet”containsthecharacteristicsofinfallofplan-
the source are compiled once. The requirements for the
etoids.
calculations are read in external files to be supplied by
– “rap iso” allows to modify the isotopic ratios.
the user. The run is interactive. Messages displayed in
– “vent”defines the chemical compositionof the wind
French, or in English, allow the control of calculations.
when it differs from the atmosphere composition.
Theextentofconvectionzones,aH–Rdiagram,andthe
– “zoom” allows to fit the resolution of the display for
profiles of temperature, pressure, luminosity and abun-
the plot on line.
dancesaredisplayed1 online.cesamhasbeenespecially
– “langue”allows to have the messages in English.
designedtofacilitatetheimplementationofvariousphys-
ical constants, opacities, EOS, atmosphere, nuclear net- The meaning of all items are explained in Paper 2. Ex-
works etc. So, its overall structure is separated in two amples of files are given in the directory EXPLOIT.
spaces:
1. A “physical space” where the coefficients of the dif- Output files. At the end of each time step, a “return
ferentialequationsarewritteninaformclosetotheir binaryfile”(hereafterRBF)iscreated.Itcontainsallthe
physical formalism. dataneededtoinitialiseorpursueacomputation.There
2. A “numerical space” wherethe differentialequations is the possibility, either to save all RBF, or to keep only
are formally solved. the last RBF created. On request, “output data files”
Thereforecesamallowstoimplementphysicalprocesses, (hereafter ODF) are created. Three ODF are designed
and physical data, without any knowledge of numerical for adiabatic, non adiabatic and inversion asteroseismic
methods involved for the solution of the equations. For investigations. These ODF also serve as input for some
thephysics,theuseofgenericroutinesmakesthereading programs of the directory EXPLOIT. An ODF concerns
of the algorithms easier. the diffusion of the angularmomentum (β). All items of
output files are detailed in Paper2. It is also possible to
Units and values of physical constants. cesam uses cgs create personalised ODF.
units except for the mass, radius and luminosity ex-
pressedinsolarunits.Twosetsoffundamentalconstants
Thekindsofprecision. Tooptimisethecalculations,sets
areimplementedwhichcorrespondtowidelyusedvalues
of parameters, named “r´eglages”, are fixed according to
(Clayton,1968;Christensen-Dalsgaard,1988;Lide et al.,
the kind of models to be calculated and to their subse-
1994; Cox, 2000). For each calculation, one of these sets
quent use. The most useful r´eglages are:
ischosenastheunique sourceoffundamentalconstants.
Other constants are initialised locally, e.g. the mass ex- – “realistic precision”for standard evolutions.
cesses,intheroutineperformingthecalculationsofther- – “super precision”usedwhenahighlevelaccuracy
monuclear reaction rates. is needed.
– “solar accuracy”, close to super precision, but
Input files. Very often only the “input data file” (here- especially designed for seismological investigations;
afterIDF)isneeded.Itisreadattheonsetoftherunand the number of shells of the last model is increased
collects allthe requirements needed for the calculations: up to its maximal value.
– physicalparameters:mass,chemicalcomposition,mi- – “corot”, close to super precision, but especially
xing-length parameter, etc. designed for investigations connected to CoRot.
– numerical parameters: maximum number of shells, – “advanced”, for the computation of early type star
kind of precision, etc. models evolved up to the oxygen burning.
– criteriaforhaltingthecomputations:agetobereached, – “normal precision”,for exploratory work.
value of the hydrogen abundance at centre, etc. – “low mass”, for the computation of late type star
– names and locations of the external data files con- models.
taining the data of the tabulated EOS and opacity, – “reglages”, in that case, the parameters, designed
names of the physical routines to be used, name of bytheuser,arereadonaninputfilenamed“reglages”.
the model, of the set of units to be used, etc.
The largerthe expected accuracy,the largerthe compu-
1 Use of thePGPLOT package. tational expense.
CESAM: a free code for stellar evolution calculations 3
Operating data and programs. The directory EXPLOIT ofthe relativechangesofheliumonthe wholestar.Oth-
contains: erwise, during the pre-main sequence, the limitation of
the change of mechanical energy controls the time step.
– examplesofinputfilesasIDF,“reglages”,“planet”,
In case of divergence of any iterative algorithmthe time
“modif mix”, etc.
step is halved.
– ASCII files of preliminary models for the initialisa-
tions of PMS and ZAMS models.
Stop criteria. According to flags, read in the IDF, the
– miscellaneous programs to make plots of the chem-
computations may be stopped:
ical composition profile, or an extension of the grid
– when the expected age is reached.
on a given set of radius, or create the IDF for solar
– when the central temperature reaches a given value.
calibration, etc.
– assoonastheabundanceofhydrogenatcentrereaches
The flow chart of cesam is described in Sect. 2. As the a given value.
numericalfeaturesaredetailedin the appendix ofMorel – at the exhaustion of hydrogen at centre.
(1997)2 (hereafter Paper I) they are only succinctly re- – when the helium core reaches a given extent.
called in Sect. 3, except for the automatic allocation of – whenthe effective temperature crossesagivenvalue.
meshpointsdescribedinSect. 3.1.Therestitutionofthe – at the ignition of the 3α cycle.
atmosphere is outlined in Sect. 4. The algorithms per- – at the ignition of the carbon cycle.
forming the temporalevolutionare describedinSect. 5. – at the ignition of the oxygen cycle.
The nuclear network is detailed in Sect. 6 and the im-
Formostkindsofprecision,thelasttimestepisadjusted
plementationofthe rotationisdescribedinSect. 7.The
in order to fulfil the required stop condition.
variousformalismsofconvectionimplantedincesamare
described in Sect. 9. Mass loss formalisms and infall of
planetoids are described in Sect. 8. EOS and opacities
3 Numerical methods.
data available are listed in Sect. 10.
Choice of variables. Forthenumericalintegrationofthe
stellar structure equations, the Lagrangian form is the
2 The flow chart. most convenient as the discretisation on mass is readily
expressed. However, it presents a singularity at centre
Initialisations. At the onset of the run, the IDF is read. and the core needs to be integrated apart.The Eulerian
Chemical composition is initialised according to the ini- formof the equationsdoes notsuffer fromsuch inconve-
tial mixture and to the isotopes used by the chosen nu- nience,butsincethestellarradiusvarieswithtime,there
clear network. Then, using fit-formulas (see Sect. 6.1), isafreeboundary(seeStoer & Bulirsch,1979,par.7.3).
thethermonuclearreactionsratesaretabulatedonarel- With the Lagrangianvariables:µ M2/3, R2,L2/3, the
evant interval of temperatures3. Then, evolutionbegins: centralsingularitydisappears(see≡Paper I,par.B1)and
there is no need of a special treatment for the core (M
– either it starts from zero age on PMS or ZAMS: an
is the mass, L the luminosity and R the radius). To be
initial modelhavingthe requiredspecificationsis de-
consistent, the chemical species are taken as functions
duced from a model taken from a RBF or from a
of M2/3. The pressure, P, and the temperature, T, are
model in ASCII chosen in the directory EXPLOIT.
expressed in logarithms, (ξ = lnP, η = lnT), on the
– or it pursues a previous calculation, then the input
ground that they change by more than six magnitudes
is one RBF of the evolution going on.
from the centre to the atmosphere.
The evolution. The number of shells is updated as ex- Solvingthedifferential equations. Theunknownsareap-
plained in Sect. 3.1. Then, taking possible overshooting proached by piecewise polynomials of order defined ac-
intoaccount,thelimitsbetweenradiativezonesandcon- cording to the required accuracy;the mostly used are of
vective mixed zones (hereafter LMR) are localised. The order1,i.e.linearpiecewise,andoforder2,i.e.parabolic
angular velocity (β) and the chemical composition are piecewise. For the stellar modelling, such a flexible rep-
thenupdated.In fine,the equationsofquasi-staticequi- resentationiswelladaptedtothepresenceofdiscontinu-
librium for the interior and the atmosphere are solved. ities resulting from the mixing of the convection zones.
The process is repeated until convergence. For the calculations, the piecewise polynomials are pro-
When the nuclear engine is at work, the time step jected on a local linear basis of normalised B-splines
control is first based on the local accuracy achieved for (De Boor, 1978; Schumaker, 1981). That allows to find
the numerical integration of species of relevant interest back exactly the solution at any location. Moreover, B-
(see Sect. 6.2andSect. 6.3)and,secondona limitation splines basis are also used for solving:
2 Only available on electronic form. – thetwopointsboundaryinitialvalueproblemsofthe
3 The errors introduced by these interpolations remain stellarstructureandoftheatmospherebycollocation
within theerror bars of thedata. (De Boor, 1978, ch. XV).
4 P.Morel, Y.Lebreton
– the diffusion equations of chemicals and angularmo- Choice of Q. Q should be a strictly monotonous, two
mentumbyfiniteelements(Quarteroni & Valli,1994). times differentiable function as simple as possible. By
experiments,ithasbeenfoundthatthemostconvenient
The linear systemsinvolvedby the resolutionofimplicit compromise is:
equations are band-diagonal.
However, due to the non-trivial and unfamiliar algebra Q=∆ξξ+∆ηη+∆µµ, ∆ξ =∆η = 1, ∆µ =15.
−
of B-splines, the algorithms are much more elaborated
where ξ and η have been defined above.
than with the finite differences. Furthermore, efficient
and stable algorithms have been constructed for inte-
Mesh refinement. The initial value of C(t) is fixed ac-
gration, differentiation, integration of differential equa-
cordingtotheexpectedlevelofaccuracy.Alongtheevo-
tions by collocation and for interpolation. In cesam the
lution, its value is kept within 2%, of its initial value
routines, especially constructed to manage the calcula- ±
by increasing (or decreasing) the total number of shells.
tions with B-splines, are derived from the algorithms of
Schumaker(1981, chap. 4). Details are given in Paper I.
Setting a grid point on a LMR. At the limit between a
mixed zone and a radiative one, there may be a singu-
larity of chemicals. Therefore each LMR has to coincide
3.1 Moving grid, mesh refinement and discontinuities precisely with a grid point. To do that cesam uses a
tracking. weighted repartition function. On each side of a LMR,
theweightsarecomputedinsuchawaythattheyadjust
An automatic mesh refinement is implemented. At time locallythe values ofC(t) tothe amountsjust neededfor
t, the locations of the n(t) mesh points are set by ful- iteratively “pushing” the closest grid point on the limit.
filling the condition that, from a grid point to the next, In most cases, including mixing, the distances from the
thejumpofastrictlymonotonous“repartitionfunction”, LMR locations to the closest grid points are lower than
Q(µ,t), is equal to a “repartition constant” C(t) (see a few per cents of the characteristic local grid size. The
Eggleton, 1971; Press et al., 1986, sect. 16.5). The loca- more well defined the location of the LMR, the most
tions of grid points, µ ,i = 1,...,n, known at the issue efficient the algorithm.
i
of the computations, satisfy:
The grids. cesam uses several grids for the B-splines
Q(µi+1,t) Q(µi,t) C(t), i=1,...,n 1. (1) representationsofquasi-staticvariables,atmosphereand
− ≡ −
rotationvariables.Asitcoincideswiththeadjustablegrid
ThechoiceofQ(µ,t)isbasedonanaprioriknowledgeof
of quasi-static variables, the Lagrangian grid used for
the behaviour of the solution. For each t, one defines an
the abundances of chemicals, is not fixed with respect
“index”function q(µ,t) mapping [0,µ ] on[1,n];the in-
b to time. There is the possibility to use a “fixed grid”, to
dex1(resp. n)correspondstothecentreandtheindexn
avoidthe numericaldiffusion resulting fromthe variable
2
tothesurface,i.e. µ =M3 .Therefore,theintegration Lagrangiangrid.
b ext
ismadeonanequidistant grid.Intermsofthederivative
of Q with respect to q, Eq. (1) reads:
4 Atmosphere.
∂ψ ∂Q
=0, with ψ(t) .
∂q ≡ ∂q The atmosphere connects the convective optically thick
(cid:18) (cid:19)t (cid:18) (cid:19)t
outer part of the envelope to the optically thin inter-
The change of variables µ q(µ,t):
→ stellarmedium. Inthe interior,the diffusionapproxima-
tion (Kippenhahn & Weigert, 1991) is used to simplify
∂µ ∂Q
ψ(t)=θ , with θ(µ,t) the calculation of the radiative flux. In the outermost
∂q ≡ ∂µ
(cid:18) (cid:19)t (cid:18) (cid:19)t parts, this approximation is no longer valid, as soon as
theRosselandopticaldepthislowerthan 10:aspecial
iscalculatedfromtheanalyticformofQ(µ,t).Thereare ≃
treatment is needed to restore the atmosphere.
two more unknowns: ψ(t) and µ(q,t); they fulfil a sys-
cesam restores the atmosphere from a T(τ,T ,g) law
temofdifferentialequationsoffirstorderwithboundary eff
here, τ is the Rosseland optical depth, T the effective
conditions: eff
temperature andg the gravity.The stellar radius,R ,is
⋆
∂µ ψ ∂ψ q =1, µ=0 definedasthe bolometricone,i.e. the radiusatthe level
= , =0, with
∂q θ ∂q q =n, µ=µb. wherethelocaltemperatureisequaltoTeff (Morel et al.,
(cid:18) (cid:19)t (cid:18) (cid:19)t (cid:26) 1994). For genuine radiativeT(τ) laws,whatever is T ,
eff
The differential equations of internal structure and of τ is located at a fixed optical depth, e.g. τ = 2/3 for
⋆ ⋆
atmosphere, written with respect to q, are detailed in the Eddington’s law. For more precise laws, the value of
Paper I.Theequationsarethensolvedonanequidistant τ isnotfixed.Therefore,thelocationwhereτ isdefined
⋆ ⋆
grid, that allows numerical optimisations. is a free boundary.
CESAM: a free code for stellar evolution calculations 5
Asoscillationmodesarereflectedintheoutermostparts Hayashitrack,thestarisfullyconvective,thereforeisen-
of stars, the pressure, the temperature and their gradi- tropic and chemically homogeneous. The energy equa-
entsshouldbecontinuousatthelimitwiththeenvelope. tion is reduced as in Iben (1965):
The continuity of the pressure gradient is trivially in-
∂L ∂S
suredbytheequationofquasi-staticequilibriumverified =ǫ = T =c(t)T (3)
G
on both sides of the limit. The continuity of the tem- ∂M − ∂t
perature gradient is more intricate to fulfil as the con- where c(t) is the “contraction constant”, and S the en-
nection occurs in zones of convective instability. For the tropy.Alongtheintervaloftime∆t,theenergyradiated
restorationof the atmosphere,the temperature gradient equals the change of gravitational energy, therefore, at
isderivedfromtheT(τ)lawitself.Fixingthegravitythe first order:
T(τ) law is expressed as: T4 = 3T4 f(τ).
4 eff L(t)+L(t+∆t) GM2 GM2
∆t
3 PκL 2 ∼ R(t+∆t) − R(t) ⇒
Using dτ = κρdr and = , (cid:18) (cid:19)
− ∇rad 16πacGMT4 2GM2[R(t) R(t+∆t)]
∆t − .
after some calculations: ∼ (L(t)+L(t+∆t))R(t)R(t+∆t)
dlnT PκdT R2 here R(t) is the stellar radius. An estimate of the initial
= = (2)
time step is deduced from two models computed with
∇≡ dlnP T dτ GM
2 2 Eq. (3), with close values of the contraction constant,
3 PκL R df L R df
⋆ = c(t)andc(t+∆t).cesamuses:c(t+∆t)=1.1 c(t).With
16πacG MT4 (cid:18)R⋆(cid:19) dτ ∇radL⋆ (cid:18)R⋆(cid:19) dτ c(t)=0.02cgsthetemperatureatcentreisa×bout105K;
itistentimeslargerwithc(t)=0.00008cgs.Changesof
HereκistheRosselandmeanopacity,ρthedensity,athe
the contraction constant allow to choose between differ-
radiation constant, G the gravitationalconstant, L the
⋆
entinitialPMSmodels.MostofthePMSmodels canbe
luminosity of the star, c the speed of light and the
rad
∇ initialised with a preliminary model in ASCII available
radiativegradient.Inprinciple,thecontinuityisensured
in the directory EXPLOIT. cesam assumes that a PMS
when the same convectiontheory prevails onboth sides.
model becomes a ZAMS model as soon as the release of
As in most cases the values differ, the continuity of the
gravitationalenergybalances the thermonuclearnuclear
temperaturegradientisinsuredbyaweightedmeanwith
one. Such a model is chemically inhomogeneous.
respect to lnτ.
Eq. (2) is no longer valid for a genuine radiative T(τ)
law, as it ignores convection. In such a case, following Initial homogeneous ZAMS model. A model of ZAMS
the prescription of Henyey et al. (1965) (M. Gabriel, J. with homogeneous chemical composition is not a physi-
Christensen-Dalsgaard, priv. comm.), the temperature calrealityasthenuclearenginedoesnotworkatequilib-
gradientin the convective atmosphere is computed with rium.However,itisaveryconvenientshortwayas,after
a modified radiative gradient: a few time steps, the model is very close to the model
at the end of the PMS. Several ZAMS initial models in
df
∗ = . ASCII are available in the directory EXPLOIT.
∇rad ∇raddτ
For a genuine radiative T(τ) law:
6 Evolution of chemicals.
df
lim =1,
τ≫1dτ 6.1 The nuclear networks.
therefore, at the limit between the atmosphere and the
envelope,theradiativegradientiscontinuousandconse- At the onset of the computations, the abundances of
quently the convective temperature gradient. chemicalsareinitialisedaccordingtotheinitialhydrogen
The numerical integration of the differential equations and helium mass ratios and mixture. The initial abun-
fulfilled in the atmosphere is made by collocation, see dance of each chemical species is split between its iso-
Paper I for details. The number of shells in the atmo- topes, according to the isotopic ratios of nuclides. Sev-
sphere is fixed from 50 to 100 according to the required eral mixtures are implemented, the most useful are: the
level of accuracy. Anders & Grevesse (1989) meteoritic mixture and the
solarmixturesofGrevesse & Noels(1993)andGrevesse & Sauval
(1998). If necessary an IDF allows to modify the initial
abundances of specific species or to use a mixture not
5 Evolution of the internal structure. yet implemented.
Upto16nuclearnetworksarepresentlyavailable.Hence,
Initial PMS model. TheenergysourceinaninitialPMS one can follow the evolution using only the chemical
model is only of gravitationalorigin.At the onset of the species and the thermonuclear reactions of interest.
6 P.Morel, Y.Lebreton
The nuclear reaction rates are tabulated on relevant in- here, is the gradient operator and v the advection
i
∇
tervals of temperatures. The rates are computed using velocity.Thesymbol“ ”meansthevectorinnerproduct
•
theformulasofCaughlan & Fowler(1988)orofNACRE
T
∂x ∂x
compilations(Angulo et al.,1999).Forsolarmodels,the and = 1,..., nX .
M
improved rates of Adelberger et al. (1998) are available. ∇ X ∂M ∂M
(cid:18) (cid:19)
TheweakscreeningofSalpeter(1961)andtheweakand Theturbulentdiffusioncoefficients,areaddedtothei-th
intermediate screenings of Mitler (1997) are available. component of the vector, D , of diffusion coefficients of
i
the species i, see Eq. (16).
For the integration, the abundances are approached by
piecewisepolynomialsexpressedonaB-splinebasiswith
6.2 Evolution without diffusion.
discontinuousderivativesateachLMR.Afinite-elements
method (see e.g. Quarteroni & Valli, 1994) is used to
The time scales involved in the temporal evolution of
solve the diffusion equation. That allows an integration
chemicalsdifferbyalargenumberofmagnitudes.Froma
bypartswhichreducestounitytheorderofthediffusion
mathematicalpointofviewitisastiffproblem(Hairer & Wanner,
equation.The scheme is fully implicit. The nuclearterm
1991) and algorithms have been especially designed for
is evaluated as for the implicit Euler’s formula.
it.Withoutmicroscopicdiffusion,L-stable(Hairer & Wanner,
The mixing is made by turbulent diffusion with coef-
1991, par. 4.3) implicit Runge-Kutta schemes, are avail-
ficient D >> 1. At each LMR, the abundances x
able for the chemical evolution. MZ i
and fluxes F are continuous functions with discontinu-
In the radiative zones, the equations to be solved are i
ous first derivatives, owing to the jumps of the diffusion
formally written:
coefficients. Therefore Eq. (6) holds everywhere. Two
∂x
i =Ψ (P,T, ,t), 1 i n . (4) formalisms are available for the calculation of the diffu-
i X
∂t X ≤ ≤ sion vector D :
i
Here,x istheabundancepermoleofthechemicalspecies
i
i, Ψ the rate of change of x , the vector of chemical – ThecoefficientsarecalculatedaccordingtoMichaud & Proffitt
i i
abundances,t thetime andn Xthenumberofchemicals; (1993).Themetalsare“testelements”,theirdiffusion
X
x is given by: onlyresultsfromcollisionsagainstprotons.Basedon
i
the presence of protons, this formalism is only valid
X
i
xi = , for the main sequence.
ν
i – The diffusion coefficients are computed according to
X is the abundance mass fraction and ν the atomic
i i Burgers (1969), this formalism is outlined beneath.
mass. In mixed zones (hereafter MZ), convective eddies
homogenise the chemical composition. There, the mix-
Boundaryconditions. Attheoutermostlimit,M =M ,
ing and updating of chemicals are done simultaneously, ext
it is assumed that there is neither input nor output of
therefore the changes of mean abundances read:
matter, then F (M ) = 0 for any particle i. At centre
i ext
dx¯
i = Ψi(P,T, ¯,t)dM dM. (5) M =0, because of the spherical symmetry, Fi(0)=0.
dt X
ZMZ (cid:30)ZMZ
As a grid point is defined on each LMR (see Sect. 3.1), Control of the gravitational settling. For stellar models
the discontinuities of the abundances are explicitly cal- with mass larger than 1.4M , the use of microscopic
⊙
≈
culated. diffusion alone produces an important depletion of he-
lium and metals at the surface and a concomitant en-
Control of the accuracy. A good estimate of the nu- hancement of hydrogen.Different waysto overcomethis
merical accuracy of an integration is obtained with the problemhavebeenused.Eggenberger et al.(2005)intro-
Fehlberg method (Stoer & Bulirsch, 1979, par. 7.5.2). duce some turbulence due to rotation. Di Mauro (2004)
It needs to triple the calculations. As it is prohibitive, suppress diffusion in the outer layers. Chaboyer et al.
falling anything better, the time step is simply adjusted (1999) include a wind mass loss which reduces the dif-
insuchawaythat,overatimestep,therelativechanges fusion in outer layers,Turcotte et al. (1998) introduce a
oftheabundancesremainwithinfixedlimits.Thelargest turbulent mixing. cesam allows to control gravitational
the expected accuracy, the narrowerthe limits.. settlingwitharadiativeturbulenceofcoefficientd pro-
ν
portionaltotheradiativekineticviscosity;itresultsfrom
the energy exchangesbetween thermal collisions leading
6.3 Evolution with diffusion.
to excitationandionisationofatoms and ions (Thomas,
1930; Mihalas & Weibel-Mihalas, 1984, p. 461-472).
With microscopic diffusion, the equations of the evolu-
tion of chemicals have the form: 4 aT4
d =Re , (7)
∂xi ∂Fi ∂xi ν ν15cκρ2
= + (6)
∂t ∂M ∂t
(cid:18) (cid:19)nucl. c is the speed of light in vacuum. The phenomenolog-
F =4πR2ρ 4πR2ρD +v x ical parameter Re has been found close to unity by
i i M i i ν
•∇ X
(cid:0) (cid:1)
CESAM: a free code for stellar evolution calculations 7
Morel & Th´evenin (2002).Thephysicalmeaningofthis, K are the so-called resistance coefficients, they repre-
ij
asefficientassimplesourceofturbulentmixing,hasbeen sent the effects of collisions between particles i and j
questioned by Al´ecian & Michaud (2005). (Michaud & Proffitt, 1993):
16
K = n n m Ω(11), (12)
Burgers’sflowequations. Withrespecttotheabscissax, ij 3 i j ij ij
the density number n and w , the diffusion velocity of
i i with m = (m m )/(m +m ) as the reduced mass of
ij i j i j
the particle i, are related by Iben & MacDonald (1985):
particles i and j. The heat flux terms involve additional
collision integrals:
∂n 1 ∂ ∂n
i = (x2n w )+ i ,
∂t x2∂x i i (cid:18) ∂t (cid:19)nucl. z =1 2Ωi(j12),
With the formalism of Burgers (1969), the diffusion ve- ij − 5Ω(11)
ij
locity is expressed as:
25Ω(12) Ω(13) Ω(22)
∂xj z′ =2.5 ij − ij , z′′ = ij ,
wi = bij ∂x +vi. ij − 5 Ω(11) ij Ω(11)
ij ij
j
X
where the diffusion coefficient of the particle i with re- Paquette et al.(1986)showedthatΩ(kl) canbe written:
ij
spect to the particle j, b , and the advection velocity,
ij
vi, come from the solution of a linear system. The dif- Ωi(jkl) ≡Fi(jkl)ǫij, (kl)=(11), (12), (13), (22),
fusion velocities, for ions and electrons, (Burgers, 1969;
Cox et al., 1989; Thoul et al., 1994), satisfy4: e4 π z¯i2z¯j2 1
ǫ .
dP ρ dP ij ≡ 4 r2k3 √mij √T3
i i
n z¯eE = K (w w )+
i i ij j i
dR − ρ dR − − For attractive and repulsive screened Coulomb poten-
Xj tials, the quantities ln(F(kl)) have been tabulated by
m r m r ij
+ K z j i− i j (8) Paquette et al. (1986). The equations of dynamical con-
ij ij
m +m
i j servation of the mass and charges respectively are:
j
X
5n kdT = 5 K z mj(wj −wi) 2K z′′r nX nX+1
2 i dR − 2 ij ij mi+mj − 5 ii ii i xiνiwi =0, z¯ixiwi =0, (13)
j
X i=1 i=1
3m2+m2z′ + 4m m z′′ X X
−ri Kij i (mjii+j mj5)2 i j ij (w1h9e8r6e)naXn+d1otishethrewionrdkesx,oe.fge.leIcbternon&s.MFoallcoDwoinnaglPda(q1u9e8t5t)e,et al.
j6=i
X Cox et al. (1989), Thoul et al. (1994), the diffusion ve-
m m (3+z′ 4z′′)
+ K i j ij − 5 ij r . (9) locitieswi comefromthesolutionofthesystemof2nX+
j6=i ij (mi+mj)2 j 3linearequations,formedbythe2nX Eq. (8)andEq. (9)
X for the ions, the Eq. (9) for the electrons and the two
e is the electron charge, E the electric field, Pi the par- Eq. (13).The2nX+3unknownsarewi, i=1,...,nX+
tial pressure, mi = νimu is the mass of particle i, k the 1, ri, i=1,...,nX +1 and E.
Boltzmann constant, and ri the magnitude of the resid- Forwantofsomethingbetter,withtheidealgaslaw,the
ual heat flow vector; ρi the partial density is given by: pressure and the partial pressures respectively are:
ρ =ρX =ρx ν =n ν m , (10) ρ T
i i i i i i u P = R , P =n kT =ρx N kT =Pµx , (14)
i i i 0 i
µ
with N as the Avogadro number, the inverse of the
0
here istheperfectgasconstant,andµthemeanmolec-
atomic mass unit m :
u R
ularweight.Insphericalsymmetry,the pressureandthe
Xi temperature gradients are given by:
n =ρN =ρN x . (11)
i 0 0 i
ν
i dP dT T ∂lnT d T
= ρg, = P = ρg. (15)
The force equation (8) represents the pressure and the dR − dR ∂ln dR −P∇
P P
concentrationdependenceofthediffusionvelocity,while
WorkingwithEq. (10)toEq. (15),withrespecttoµand
Eq. (9) prevails for its thermal dependence. The charge
x , the Burgers’s equations (8), (9), may be rewritten:
of any isotope is taken as the averaged5 charge, z¯, over i
i
all its ionisation states. An unique mean charge is used Aω =γ+GD , D =(∂x1,...,∂xnX,0,...,0)T.
for all the isotopes of a given chemical. The quantities x x ∂R ∂R
The solutions are:
4 Without conduction current and magnetic field.
5 Weighted bythe ionisation rates. ω =V +BDx, V =A−1γ, B =A−1G.
8 P.Morel, Y.Lebreton
For abridgment, neither the derivation of above equa- Table 1 Comparison between the mean charges from the
tions, nor the complicated forms of vector γ and matrix approximate solution and the precise equation of state of
Gabriel (1997) (label “G”) beneath the convection zone and
AandG,arereproduced,alldetailsaregiveninPaper3.
at centre,for a calibrated solar model, .
The diffusion velocities of the ions are expressed as:
nX ∂x T=2.2MK centre
j
w = b +v , i=1,...,n ,
i ij i X
j=1 ∂R elements Z¯G Z¯ Z¯G Z¯
X
and, owing to Eq. (11), the diffusion vector reads:
C 5.89 5.92 6.00 6.00
D =(x b ,...,x b )T, (16) N 6.75 6.82 7.00 7.00
i i i1 i inX
O 7.47 7.56 8.00 8.00
for i,j =1,...,n , b and v , are respectively the coef- Fe 17.0 16.7 24.2 24.2
X ij i
ficients of matrix B and vector V.
Calculation of mean charges. TheBurgersformalismin- the g are the statistical weights of fundamentals. For
j,i
volves the charges of the isotopes. To simplify, cesam the smoothing function f(x), cesam uses the piecewise
considers a unique mean charge for all the isotopes of cubic polynomial with zero derivatives at x = 0 and
eachchemical.Forthecalculationoftheionisationrates, x=p (p=4):
theSaha-Boltzmannequation(Cox & Giuli,1968,eq.15-
0 if x 0,
30) has been adapted in the following way. Let n be
j,i 2 ≤
the number density of atoms in ionisation state j of the f(x)= x 2 x +3 if x [0,p],
chemicalspeciesi.Theratioofthetotalnumberofatoms (cid:16)1p(cid:17) h− (cid:16)p(cid:17) iif x∈ p,
in successive stages of ionisation can be written: ≥
nj−1,i = Uj−1,i exp η+ χj,i , (17) gToheesn,toaszesoroonanasdχtih,je≤levχeCi,lj,jthe1qbuaecnotmityesgjf−ul1l,yi/giojn,ifis(exd).
n U kT −
j,i j,i ThesetofSaha-Boltzmannequationsforallionsiswrit-
(cid:16) (cid:17)
here χ is the ionisation potential, U the partition ten as in Mihalas (1978, eq. 5-17). They are solved by
j,i j,i
function and η the electron degeneracy, related to the iterationsusingasecondorderNewton-Raphsonscheme.
numberdensityoffreeelectronsandtemperature,through For the mean charges, quantities of main interest in the
thehalfintegerFermi-Diracfunction(Clayton,1968,eq. present investigation, Table 1 reveals agreements better
2-57). Eggleton et al. (1973) have introduced a conve- than 3% with the data of Table 1 of Gabriel (1997).
±
nient approximate treatment of pressure ionisation: a
numerical correction is postulated to ensure that the
plasma remains completely ionised at sufficiently high 7 Rotation.
density. In the inner solar radiative zone, this approxi-
mate treatment leads to too large mean charges for the Rotationisconsideredincesamundertheassumptionof
ions and, at the centre, iron is fully ionised, while it is sphericalsymmetry.With non zeroangular velocity,the
only 85%accordingto Table 1 inGabriel(1997). A sim- meancentrifugalaccelerationaffectsthelocalgravity.In
ilar behaviour is observed with the modified parameters the initial model, rigid rotation is assumed. The initial
recommended by Proffit & Michaud (1991, eq. 4). angular velocity can be read from the IDF in different
Incesam, the partitionfunctions arelimited to the sta- units:
tisticalweightsoffundamentallevelsand,assoonasthe
– radian/s.
mean distance between the ions becomes of the order
– km/s,it correspondsto the rotationalvelocityof the
of the size of the ion cloud, R , the ratio of statistical
D outer part of the star. As the outer radius depends
weights is smoothly reduced to zero. The Debye-Huckel
on the outer gravity,the initial model has to be iter-
radius writes (Clayton, 1968, eq. 2-235) *:
atively adjusted to have the required outer velocity.
– days, that corresponds to the rotation period.
kT
R = , ζ = z¯(z¯ +1)x ,
D s4πe2ρN0ζ i i i
i6=e
X 7.1 Rotation without diffusion of angular momentum.
The Saha-Boltzmann Eq. (17) is then written:
Several options are available:
n g χ
j−1,i j−1,i j,i
= f(x)exp η+ , (18)
n g kT – no rotation: the initial angular velocity must also be
j,i j,i
(cid:16) (cid:17) zero.
x= χi,j 1, χ′ =max(0,χ χC ), χC = je2, – solid-bodyrotationwheretheangularvelocityiskept
χ′ − j,i i,j − i,j i,j R to its initial value.
j,i D
CESAM: a free code for stellar evolution calculations 9
– rotation is solid and angular momentum is globally border of a convective core, the discontinuities of chem-
conserved.Theangularvelocitychangeswithrespect icals are assumed to be eroded only by diffusion. With-
to time, according to structural changes. out diffusion (β), the abundances profiles are spatially
– angular velocity changes according to the conserva- linearlyinterpolatedbetweentheirvaluesonthe convec-
tion of the local angular momentum. The rotation is tive core and at the former location of the LMR. That
not solid, but in the mixed zones. avoidsanoisy behaviourofchemicalgradientsand,con-
sequently,ofthe profileof the Brunt-Va¨issa¨l¨afrequency.
7.2 Rotation with diffusion of angular momentum.
10 Equation of state and opacity.
ThetwoformalismsofTalon et al.(1997)andMathis & Zahn
(2004)ofthediffusionofangularmomentum(β)areim- Four analytical EOS are implemented. The most useful
plemented in cesam. The diffusion coefficients of the areeff(Eggleton et al.,1973)andceff(Christensen–Dalsgaard& Dappen
angular momentum are computed according either to 1992). Numerical EOS are available. The mhd tables
Palacios et al. (2003) or to Mathis et al. (2004). (Mihalas et al.,1988)andtheopal1993andopal2001
tables are used with the OPAL interpolation scheme for
tables with Z =0.01 and 0.02.
8 Mass loss and infall of planetoids. OPAL1996opacitytablesareimplemented,withKurucz
low-temperature values. The ratios between the abun-
Several formalisms of mass loss are implemented. The dancesofheavy-elementsarefixedattheirinitialvalues,
mass loss rate is either negative or positive, i.e. increase regardless of changes due to nuclear reactions and dif-
ofmass.With diffusion,the chemicalcompositionofthe fusion. All OPAL tables of type 2, available today, are
input (resp. output) can differ from those of the outer implemented altogether with z14xcotrin21, the pack-
convective zone, from where the output (resp. input) is age of A.I. Boothroyd.
assumedtocomefrom(resp.vanish).Themasslossrate, The neutrinos emitted in nuclear reactions are taken
in unit of M .yr−1 is read in the IDF. The following into account, it is assumed that they can freely escape
⊙
options are implemented: the star. Other physical processes related to neutrino
such as Urca process, plasma neutrino, pair neutrino
– standard mass loss.
(Kippenhahn & Weigert,1991,par.18.6)haveonlybeen
– solarmassloss,the masslossishaltedassoonasthe
implemented in private codes derived from cesam.
mass of the model reaches the solar value.
For temperature values above T 7107K 7Kev, the
– mass changes due to nuclear energy generation. ≥ ≃
plasma is fully ionised, so the Rosseland mean opacity
– infall of planetoids (β). The characteristicsof the in-
is reduced to the Compton scattering by free electrons
fall, duration,amount of terrestrialmasses,chemical
(Cox & Giuli, 1968, par. 16.6).
composition of the planetoids, are read in a IDF.
Angular momentum losses are implemented but in vali-
dation (β). 11 Calibration of the Solar model.
We have performed calibrations of the solar model with
9 Convection. cesam2k with various sets of input physics and ini-
tialparameters.Thisconsistsinadjustingthe initialpa-
Two formalismsfor the computationof the temperature rameters of the model (helium abundance and mixing-
gradient in the convection zones are available: the stan- length parameter) in order to satisfy the observational
dardB¨ohm-Vitense(1958)mixing-length(MLT)formal- constraints on the solar global parameters at solar age.
ism is considered with the optical thickness of the con- We have adopted the values of the astronomicaland
vective bubble, and the Canuto & Mazzitelli (1991) for- physicalconstantsspecifiedforthecalculationofthestel-
malism.Themixing-lengthparameterisreadintheIDF. lar models compared in the different ESTA tasks (see
Overshootingbeneathand/orabovetheconvectionzones Lebreton et al., 2007). For the solar global parameters,
canbeaccountedfor.Theovershootingparameters,scaled we therefore took R = 6.9599 1010 cm (radius), L =
⊙ ⊙
by the local pressure scale height, are read in the IDF. 3.8461033erg.s−1(luminosity)andM =1.989191033g
⊙
Inovershootregions,the temperature gradientisset, ei- (mass). The value R refers to the radius of the model
⊙
ther to the adiabatic or to the radiative gradient. The layerwhereT =T =5777K.Forsolarage,weadopted
eff
convective zones and their extents by overshooting are the value t =4.57 Gyr (Bahcall & et al., 1995).
⊙
homogenised, see Sect. 6.2. Table 2 presentsthe sets ofinput physicsusedin the
Uptonowspecialisedtreatmentofthesemi-convectionis different calibrations of solar models. All models have
not implemented in cesam. With microscopic diffusion, beencalculatedwiththeOPAL2001EOS(Rogers & Nayfonov,
in areas swept across by the backward movement of the 2002)andthe1995OPALinterioropacities(Iglesias & Rogers,
10 P.Morel, Y.Lebreton
Table 2 The sets of parameters used for the calibration of solar models with CESAM2k. We used (i) two solar heavy
elementsmixtures,theGN93(Grevesse & Noels,1993)ortheAGS05(Asplundet al.,2005)mixture;(ii)twoformalismsfor
convection treatment, the BV (B¨ohm-Vitense, 1958) and CGM (Canuto et al., 1996) formalism, (iii) two T(τ) laws for the
atmosphere,theEddingtongreylawandT(τ)lawsderivedfromATLAS9modelatmospheres(Kurucz,1992)calculatedwith
the same convection formalism as the interior models (see Samadi et al., 2006) and (iv) two formalisms for the calculation
of the microscopic diffusion of the elements, theMP93 (Michaud & Proffitt, 1993) and B69 (Burgers, 1969) formalism.
Model Mixture Convection T(τ)law Diffusion
A GN93 BV ATLAS9 MP93
B GN93 BV EDDINGTON MP93
C GN93 BV ATLAS9 B69
D GN93 CGM ATLAS9 MP93
E AGS05 BV EDDINGTON MP93
Table 3 Properties of the five CESAM2k calibrated solar models. Y0, Z0 and (Z/X)0 are the initial values of the helium
andheavyelementsmass fractions and(Z/X) ratio. Ye,Ze arethepresentvaluesintheconvectiveenvelopewhile Re isthe
radiusat thebasis oftheconvectiveenvelopeinunitsofthesolar radius. ρc,Tc and Xc arethecentralvaluesofthedensity
in g.cm−3, temperature in 106 K and hydrogen abundance.αconv,int is the valueof theconvection parameter in theinterior
(α=l/Hp where l is the mixinglength and Hp is thepressure scale-height)
Model Y0 Z0 (Z/X)0 Ye Ze Re ρc Tc Xc αconv,int
A 0.2735 0.0196 0.0278 0.2447 0.0181 0.7143 153.1 15.72 0.3387 2.42
B 0.2735 0.0196 0.0278 0.2447 0.0181 0.7145 153.1 15.72 0.3387 1.76
C 0.2741 0.0198 0.0280 0.2460 0.0180 0.7152 153.0 15.72 0.3388 2.40
D 0.2735 0.0196 0.0278 0.2447 0.0181 0.7143 153.1 15.72 0.3387 0.77
E 0.2637 0.0141 0.0196 0.2339 0.0129 0.7297 151.0 15.52 0.3578 1.72
erations)butweconsideredeithertheMichaud & Proffitt
0.015
(1993, hereafter MP93) or the Burgers (1969, hereafter
A
B B69) formalism. Convection is treated either according
C totheclassicalMLT(B¨ohm-Vitense,1958,hereafterBV)
D or to the Canuto et al. (1996, hereafter CGM) formal-
0.01 E
ism. In the CGM formalism, like in the CM formalism
(Canuto & Mazzitelli, 1991), the contribution of eddies
with different sizes is taken into account in the calcu-
c lation of the convective flux and velocity. In addition
δc/ 0.005 Canuto et al. (1996) take into account the feedback of
theturbulenceontheenergyinputfromthesourcewhich
generates turbulent convection. For the atmosphere cal-
culation we considered either the classical Eddington
0 grey T(τ)-law or a law derived from Kurucz’s ATLAS 9
1-D model atmospheres (Kurucz, 1992). We have taken
thesameT(τ)-lawsasusedintheworkbySamadi et al.
(2006). These laws are based onmodel atmospheres cal-
-0.005 culated with either the BV or the CGM convection for-
0 0.2 0.4 0.6 0.8 1 mulation.In both cases,the atmosphere calculationwas
r/R performedadoptingavalue ofthe mixing-lengthparam-
SUN
eter α = 0.5 which allows to fit at best the ob-
Figure 1 Relative differences between the seismic sound conv,atm
servedprofilesofthesolarBalmerlines(seevan’t Veer-Menneret and Megessier
speedderivedbyBasu et al.(2000)andthesolarmodelspre-
sented in Table 2. 1996). Therefore αconv,atm is different from the value of
α inthe interior,this latterbeing adjustedto cal-
conv,int
ibratethesolarmodel.Finally,weadoptedtheGN93so-
1996). Two sets of low temperature opacities have been larmixtureofheavyelements(Grevesse & Noels,1993)
used: the Alexander & Ferguson (1994) tables given for inallmodelsbutonewhereweusedthenewAGS05mix-
theGrevesse & Noels(1993)solarmixtureandtheFergusontuertea(Al.splund et al.,2005)whichisderivedfromatime-
(2005) tables givenfor the new solarmixture derivedby dependent,3-Dhydrodynamicalmodelofthesolaratmo-
Asplund et al. (2005) (see below). All models take into sphere. The abundances of C, N, O of the AGS05 mix-
account the diffusion of chemicals due to pressure, tem- turearesmallerthaninthe GN93one whichleadstoan
peratureandconcentrationgradients(noradiativeaccel-
