{'fraction_non_alphanumeric': 0.02926843089476965, 'fraction_numerical': 0.013940046466821556, 'mean_word_length': 4.627391581632653, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ASTEC -- the Aarhus STellar Evolution Code

The Aarhus code is the result of a long development, starting in 1974, and still ongoing. A novel feature is the integration of the computation of adiabatic oscillations for specified models as part of the code. It offers substantial flexibility in terms of microphysics and has been carefully tested for the computation of solar models. However, considerable development is still required in the treatment of nuclear reactions, diffusion and convective mixing.


Introduction
What has become ASTEC started its development in Cambridge around 1974. The initial goal was to provide an improved equilibrium model for investigations of solar stability, following earlier work by Christensen-Dalsgaard et al. (1974). However, with the initial evidence for solar oscillations and the prospects for helioseismology (Christensen-Dalsgaard and Gough, 1976) the goals were soon extended to provide models for comparison with the observed frequencies. Given the expected accuracy of these frequencies, and the need to use them to uncover subtle features of the solar interior, more emphasis was placed on numerical accuracy than was perhaps common at the time. The code drew some inspiration from the Eggleton stellar evolution code (Eggleton, 1971), which had been used previously, but the development was fully independent of that code. An early description of the code was given by Christensen-Dalsgaard (1982), with further extensive details provided by Christensen-Dalsgaard  (1978); many aspects of this still stand and will only be summarized here. With the increasing quality and extent of the helioseismic data the code was further developed, to allow for more realistic physics; a major improvement was the inclusion of diffusion and settling (Christensen-Dalsgaard et al., 1993). This led to the so-called Model S of the Sun  which has found extensive use as reference for helioseismic investigations and which at the time provided reasonably up-to-date representations of the physics of the solar interior. In parallel, extensions have been made to the code to consider the evolution of stars other than the Sun; these include the treatment of convective cores and core overshoot, attempts to model redgiant evolution and the inclusion of helium burning. This development is still very much under way.
For use in asteroseismic fitting a version of the code has also been developed in the form of a subroutine with a reasonably simple calling structure which also includes the computation of adiabatic oscillation frequencies as part of the computation. The combined package is available as a single tar file, making the installation relatively straightforward, and the code has been successfully implemented on a variety of platforms. Nevertheless, it is sufficiently complex that a general release is probably not advisable.

Formulation of the equations
The equations of stellar structure and evolution are written with x = log 10 q as independent variable, where log 10 is logarithm to base 10; here q = m/M is the mass fraction where m is the mass interior to the point considered and M is the photospheric mass of the star, defining the photosphere as the point where the temperature is equal to the effective temperature. The equations are written on the form Here r is distance to the centre, ρ is density, p is pressure, G is the gravitational constant, T is temperature, L is the luminosity at r, ε is the energy generation rate per unit mass, H is enthalpy per unit mass, X k is the mass fraction of element k, R k is the rate of change of X k due to nuclear reactions, and D k and V k are the diffusion and settling coefficients for element k. Also, the value of ∇ = d log T /d log p depends on whether the region is convectively stable or not. In the case of a convectively stable region, where a is the radiation density constant,c is the speed of light, and κ is the opacity. The factor ψ is discussed following Eq. 11 below. The calculation of the temperature gradient in convective regions is discussed in Section 4. In Eq. 5 the derivatives with respect to m should obviously be expressed in terms of derivatives with respect to x. As discussed below, it is convenient to introduce the diffusion velocity on the form then Eqs 5 can be written as Here M = log 10, log being the natural logarithm. (In the codeỸ k = M −1 Y k is used as variable.) For elements where diffusion and settling are ignored the original form, Eq. 5, with D k and V k set to zero is obviously used.

Boundary conditions
The centre, r = 0, is obviously a singular point of Eqs 1 -4. As discussed in considerable detail by Christensen-Dalsgaard (1982), this is treated by expanding the variables to second significant order in r to obtain the required conditions at the innermost point in the numerical solution. These are written in the form of expressions of the m and L at the innermost mesh point in terms of the other variables. The treatment of diffusion and settling has yet to be included consistently in this expansion; currently, conditions are imposed that essentially set theỸ k to zero at the innermost meshpoint. The outer boundary of the solution of the full evolution equations is taken to be at the photosphere, defined as the point where T = T eff , the effective temperature. Consequently an obvious boundary condition is where σ is the Stefan-Boltzmann constant. To obtain a second condition, a expression is assumed for the dependence of T on optical depth τ, conveniently written as Based on this the equation of hydrostatic equilibrium, on the form dp can be integrated assuming, currently, a constant gravitational acceleration g. The surface boundary condition based on an approximate treatment of the outer parts of the atmosphere, is applied at the top of the atmosphere, τ = τ min , typically taken to be 10 −4 . The boundary condition is obtained by equating the pressure resulting from integrating Eq. 12 to the pressure in the interior solution.
The T (τ) relation and Eq. 12 obviously define the temperature gradient ∇ in the atmosphere. To ensure a continuous transition to the interior I follow Henyey et al. (1965) and include the factor ψ = 1 + q ′ (τ) in Eq. 6.

Numerical scheme
As indicated in Eqs 1 -5 and 8, 9, the dependent variables on the left-hand sides of the equations are expressed in terms of the set {y i } = {log 10 r, log 10 p, , log 10 T, log 10 L, X k ,Ỹ k ′ } , i = 1, . . . , I .
Here the index k runs over all the abundances considered, whereas the index k ′ runs over those elements for which diffusion is included. In the solution of the equations, the computation of the right-hand sides is the heaviest part of the calculation and hence needs to be optimized in terms of efficiency. As discussed below, this may involve expressing the thermodynamic state in terms of a different set of variables {z j }, related to the {y i } by a non-singular transformation. Using also the reformulation, Eqs 8 and 9, of the diffusion equation, the combined set of stellar-evolution equations consists of equations of one of the following three types: The equations of Type III are obviously relevant to the evolution of abundances of elements for which diffusion is not taken into account. To this must be added the transformation In practice, z i = y i for i = 2; the choice of z 2 depends on the specific equation of state considered. The equations are solved on the interval [x 1 , x 2 ], with boundary conditions that can be expressed as The equations are solved by means of the Newton-Raphson-Kantorovich scheme (see Richtmyer, 1957;Henrici, 1962), in the stellar-evolution context known as the Henyey scheme (e.g., Henyey et al., 1959Henyey et al., , 1964. We introduce a mesh x 1 = x 1 < x 2 < . . . x N = x 2 in x and consider two timesteps t s and t s+1 , where the solution is assumed to be known at timestep t s . Also, we introduce with a similar notation for Λ n,s i , as well as for quantities at timestep t s+1 . Then Eqs 15 -18 are replaced by the following difference equations: Here ∆ x n = x n+1 − x n and ∆t s = t s+1 − t s . Also, the parameters θ i allow setting the centralization of the difference scheme in time, θ i = 1/2 corresponding to time-centred differences and θ i = 1 to a fully implicit scheme. The former clearly provides higher precision but potentially less stability than the latter; thus time-centred differences are typically used for processes occurring on a slow timescale, such as the change in the hydrogen abundance, whereas the implicit scheme is used, e.g., for the time derivatives in the energy equation where the characteristic timescale is much shorter than the evolution time and hence short compared with the typical step ∆t in time. The algebraic equations 24 -26, together with the boundary conditions, Eqs 21 and 22, are solved using Newton-Raphson iteration, to determine the solution {z n,s+1 i } at the new time step (see also Christensen-Dalsgaard, 1982). Given a trial solution {z n,s+1 i } the equations are linearized in terms of corrections δ z n,s+1 i and the resulting linear equations are solved using forwards elimination and backsubstitution (e.g., Baker et al., 1971). This process is repeated with the thus corrected solution as trial, until convergence.
The initial ZAMS model is assumed to be static and with a prescribed chemical composition. Thus in this case only Eqs 24 are solved. Time evolution is started with a very small ∆t for the initial non-zero timestep.
The number N of meshpoints is kept fixed during the evolution, but the distribution of points is varied according to the change in the structure of the model. The mesh algorithm is based on the first-derivative stretching introduced by Gough et al. (1975) (see also Eggleton, 1971), taking into account the variation of several relevant variables. An early version of the implementation was described by Christensen-Dalsgaard (1982), but this has subsequently been substantially extended and is still under development. In particular, a dense distribution of points is set up near the boundaries of convection zones, although so far points are not adjusted so as to be exactly at the edges of the zones. After completing the solution at each timestep the new mesh is determined and the model transferred to this mesh using, in general, third-order interpolation; linear interpolation is used near boundaries of convective regions and in other regions where the variation of the variables is not smooth.
Having computed model at timestep t s+1 the next timestep is determined from the change in the model between timesteps t s and t s+1 ; this involves a fairly complex algorithm limiting the change in, e.g., log 10 p, log 10 T and X at fixed m. To compensate for the fairly crude treatment of the composition in a possible convective core, the change of X in such a core is given higher weight than the general change in X. The algorithm correctly ensures that very short steps in time are taken in rapid evolutionary phases. In typical simple calculations, assuming 3 He to be in nuclear equilibrium, roughly 35 (100) timesteps are needed to reach the end of central hydrogen burning in models without (with) a convective core, and 100 (200) steps to reach the base of the red-giant branch. Calculations with more complex physics or requiring higher numerical precision obviously require a substantially higher number of timesteps. Evolution up the red giant branch typically requires a large number of timesteps owing to the rapid changes at fixed m in the hydrogen-burning shell 1 although the timestep algorithm has options to reduce the weight given to this region.

Microphysics
As the code has developed over the years a number of options have been included for the microphysics, although not all of these have been kept up to date or properly verified. As a general principle, the code has been written in a modular way, so that replacing, for example, routines for equation of state or opacity has been relatively straightforward. It should be noted that the use of a different set {z i } of dependent variables on the right-hand side of the equations yields additional flexibility in the replacement of aspects of the physics. Däppen and Guzik (2000) provided a review of the treatment of the equation of state and opacity in stellar modelling. A review of nuclear reactions in the solar interior, of relevance in the more general stellar case, was provided by Parker and Rolfs (1991).

Equation of state
The original version of the code used the Eggleton et al. (1973) equation of state and that remains an often used option. In this case z 2 = log 10 f , where f , introduced by Eggleton et al. and related to the electron-degeneracy parameter, is used as one of the thermodynamic variables; this allows explicit calculation of partial ionization and hence a very efficient evaluation of the required thermodynamic quantities. The formulation also includes a crude but thermodynamically consistent implementation of 'pressure ionization' (which actually results mainly as a result of the high density in deep stellar interiors) which provides apparently reasonable results. The fact that the whole calculation is done explicitly makes it entirely feasible, if somewhat cumbersome, to evaluate analytical derivatives. Up to second derivatives of pressure, density and enthalpy are provided in a fully consistent manner, whereas third derivatives, required for the central boundary condition, assume full ionization. Partial electron degeneracy is included in the form of expansions that cover all levels of degeneracy and relativistic effects.
As an upgrade to the basic EFF formulation Christensen-Dalsgaard and Däppen (1992) included Coulomb effects, in the Debye-Hückel formulation, following Mihalas et al. (1988); unlike some earlier treatments of the Coulomb effect this ensured that thermodynamic consistency was maintained; as a result, a substantial effect was found also on the degrees of ionization of hydrogen and helium, resulting from the change in the chemical potential. The computation of the Coulomb effects consequently requires some iteration, even if the EFF variables are used, and hence some increase in computing time.
However, the resulting equation of state captures major aspects of the more complex forms discussed below.
More realistic descriptions of the equation of state require computations that are currently too complex to be included directly in stellar evolution calculations. Thus interpolation in pre-computed tables is required. The first such set to be included was the so-called MHD equation of state (Mihalas et al., 1990), based on the chemical picture where the thermodynamic state is obtained through minimization of an expression for the free energy including a number of contributions. Application of this formulation to a comparison with observed solar oscillation frequencies showed a very substantial improvement in the agreement between the Sun and the model (Christensen-Dalsgaard et al., 1988). Further updates to the MHD equation of state have been made but they have so far not been implemented in ASTEC. An alternative description is provided by the so-called OPAL equation of state (Rogers et al., 1996), based on the physical picture where the thermodynamic state is obtained from an activity expansion. This has been the preferred equation of state for solar modelling with ASTEC, used, e.g., in Model S.
The early versions of both the MHD and OPAL equations of state suffered from a neglect of relativistic effects in the treatment of the electrons (Elliott and Kosovichev, 1998). This has since been corrected (Gong et al., 2001;Rogers and Nayfonov, 2002). Also, the OPAL formulation suffered from inconsistencies between some of the variables provided (e.g. Boothroyd and Sackmann, 2003;Scuflaire et al., 2007); effects of the inconsistency on the computation of adiabatic oscillations are discussed by Moya et al. (2007). This has been improved in the latest version of the OPAL equation of state, 2 which has also been implemented in ASTEC.
Interpolation in the OPAL tables uses quadratic interpolation in ρ, T and X. Typically, a single representative value of Z is used, even in cases with diffusion and settling of heavy elements, although the code has the option of using linear interpolation between two sets of tables with different Z.

Opacity
The early treatment of the opacities used tables from Cox and Stewart (1970) and Cox and Tabor (1976) with an interpolation scheme developed by Cline (1974). A major improvement was the implementation of the OPAL opacity tables Iglesias et al., 1992). With various updates of the tables this has since been the basis for the opacity calculation in the code. The most recent tables, including a variety of mixtures of the heavy elements, are based on the computations by Rogers and Iglesias (1995). Atmospheric opacities must be supplied separately; here tables by Kurucz (1991), Alexander and Ferguson (1994) or Ferguson etal (2005) have been used, with a smooth matching to the interior tables at log 10 T = 4. Electron conduction may be included based on the tables of Itoh et al. (1983).
Interpolation in the OPAL tables is carried out with schemes developed by G. Houdek (see Rogl, 1993, 1996). The tables are provided as functions of (R, T ), where R = ρ/T 3 . For interpolation in (log R, log T ) two schemes are available. One uses a minimum-norm algorithm with interpolating function defined in a piecewise fashion over triangles in the log R − log T plane (Nielson, 1980(Nielson, , 1983. The second scheme uses birational splines (Späth, 1991). In practice the latter scheme has been used for most of the model calculations with ASTEC. Interpolation in X and log Z is carried out using the univariate scheme of Akima (1970). Note that the relative composition of Z is assumed to be fixed; thus differential settling or changes in heavyelement composition resulting from nuclear burning are not taken into account.
The code includes flexible ways of modifying the opacity, to allow tests of the sensitivity of the model to such modifications. An extensive survey of the response of solar models to localized opacity changes, specified as functions of log 10 T , was made by Tripathy and Christensen-Dalsgaard (1998).

Nuclear reactions
The calculation of the nuclear reaction rates is based on the usual approximations to the reaction integrals (e.g., Clayton, 1968) using a variety of coefficients (Bahcall and Pinsonneault, 1995;Adelberger et al., 1998;Angulo et al., 1999). Electron screening is computed in the Salpeter (1954) approximation. Electron capture by 7 Be is treated according to Bahcall and Moeller (1969).
The nuclear network is relatively limited and is one of the points where the code needs improvement. This is to some extent a heritage of its origin as a solar-modelling code, as well as a consequence of the fact that pre-mainsequence evolution is not computed. In the pp chains 2 D, 7 Li and 7 Be are assumed to be in nuclear equilibrium. On the other hand, the code has the option of following the evolution of 3 He, although in many calculations it is sufficient to assume 3 He to be in nuclear equilibrium. To simulate the evolution of the 3 He abundance X( 3 He) during the pre-main-sequence phase the initial zero-age main-sequence mode assumes the X( 3 He) that would have resulted from evolution over a specified period t3 He at constant conditions, as described by Christensen-Dalsgaard et al. (1974) but generalized to allow a non-zero initial abundance. A typical value is t3 He = 5 × 10 7 yr.
In the CNO cycle the CN part is assumed to be in nuclear equilibrium and controlled by the rate of the 14 N( 1 H, γ)  Helium burning has been included in the code using the expressions of Angulo et al. (1999), and including also the reaction 12 C( 4 He, γ) 16 O. However, the code is unable to deal with helium ignition in a helium flash. Thus models with helium burning are restricted to masses in excess of 2.3M ⊙ where ignition takes place in a relatively quiet manner. Also, as in cases of hydrogen burning (cf. Sect. 4), the treatment of semiconvection that may be required in this phase has not been implemented.

Diffusion and settling
Diffusion and settling are treated in the approximations proposed by Michaud and Proffitt (1993), with revisions kindly provided by Proffitt. For completeness I give the complete expressions in Appendix A. If included, diffusion of heavy elements assumes that all elements behave as fully ionized 16 O; this is a reasonable approximation in the solar case where the outer convection zone is relatively deep, but becomes increasingly questionable in more massive mainsequence stars. Here, also, effects of selective radiative levitation should be taken into account (e.g., Richer et al., 1998;Turcotte et al., 1998); there are no current plans to do so in the code. Various approximations to turbulent diffusion can be included, partly inspired by Proffitt and Michaud (1991). In addition, the code has the option to suppress settling in the outer parts of the star, to allow modelling of diffusion and settling in the cores of relatively massive stars where otherwise settling beneath the thin outer convection zone would result in a complete depletion of the surface layers of elements heavier than hydrogen.
At present, diffusion and settling is coupled to nuclear evolution in the consistent manner of Eq. 5 only for helium. For the remaining elements taking part in the nuclear network diffusion is neglected. Correcting this deficiency is an obvious priority.

Treatment of convection
The temperature gradient in convection zones is usually computed using the Vitense (1953) and Böhm-Vitense (1958) version of mixing-length theory; for completeness, the expressions used are provided in Appendix B. The mixing length is a constant multiple, α ML H p , of the pressure scale height H p . In addition, emulations of the Canuto and Mazzitelli (1991) formulation, established by Monteiro et al. (1996), can be used.
Convective regions can obviously, at least in stars that are not extremely evolved, be assumed to have uniform composition. This can in principle be achieved by including a very high diffusion coefficient in such regions. In ASTEC this is used in convective envelopes, ensuring that they are chemically uniform. The treatment of convective cores remains a concern and an area of active development, however. Given the lack of a proper treatment of the diffusion of all elements an explicit calculation of the chemical evolution is required. This is characterized by the (assumed homogeneous) abundances X k,c of the elements. The rate of change of these abundances can be written here q cc is the mass fraction in the convective core, x cc = log 10 (q cc ), and (28) is the reaction rate averaged over the convective core. In the second term in Eq. 27 X k (x cc ) is evaluated just outside the core; this term only has an effect if there is a composition discontinuity at the edge of the core, i.e., if the core is growing and diffusion is neglected. Finally, the term inỸ k (x cc ) is obviously only included in cases where diffusion is taken into account.
In models with a growing convective core, and neglecting diffusion, a discontinuity is set up in the hydrogen abundance at the edge of the core (see also Fig. 1). This situation arises in intermediate-mass stars (with masses up to around 1.7M ⊙ ) where the gradual conversion during evolution of 16 O into 14 N causes an increase in the importance of the CNO cycle in the energy generation (for more massive stars the CNO cycle dominates even with the initial 14 N abundance). Since pressure and temperature are obviously continuous, there is also a discontinuity in density and opacity κ, leading also to a jump in the radiative temperature gradient ∇ rad (cf. Eq. 6). Since κ increases with X and ρ, while ρ decreases with increasing X, it is not a priori clear how κ and hence ∇ rad react at the discontinuity; in practice the common behaviour is that ∇ rad increases going from the value of X in the convective core to the higher value in the radiative region just outside it. This raises the question of the definition of convective instability: if the border of the convective core is defined using the composition of the core, the region immediately outside the core is therefore convectively unstable. As a consequence, ASTEC defines the border of the convective core by the abundance in the radiative region, leaving a small convectively stable region just below the border, which nevertheless is assumed to be fully mixed in the standard ASTEC implementation. This may be regarded as an example of semiconvection, of somewhat uncertain physical consequences (e.g., Merryfield, 1995). A different scheme, now under implementation, is discussed in Section 6.
Various options exist for the treatment of convective overshoot. Simple formulations involve compositionally fully mixed overshoot regions from the convective core or convective envelope, over a distance α ov H p below a convective envelope, or α ov min(r cc , H p ) above a convective core, where r cc is the radius of the core. The overshoot region may be taken to be either adiabatically or radiatively stratified. A more elaborate treatment of overshoot from a convective envelope has been implemented, following Monteiro et al. (1994), where various forms of the temperature gradient can be specified, still assuming the overshoot region to be fully mixed. This is being extended to emulate the overshoot simulations by Rempel (2004).

Implementation details
When computing models of the present Sun it is important to adjust the input parameters such as to obtain a model that precisely matches the observed solar radius, luminosity and ratio Z s /X s between the abundances of heavy elements and hydrogen, at the present age of the Sun. This is achieved by adjusting the initial hydrogen and heavy-element abundances X 0 and Z 0 as well as the mixing-length parameter α ML (or another parameter characterizing the treatment of convection). In ASTEC the iteration to determine these parameters is carried out automatically, making the computation of solar models, and the exploration of the consequences of modifications to the input physics, rather convenient.
The ASTEC code has grown over three decades, with a substantial number of different uses along the way. This is clearly reflected in the structure of the code as well as in the large number of input parameters and flags that control its different options. These are provided in an input file, in many cases using simply the defaults provided in the source of the code. Also, several different executables can be produced, reflecting partly the evolution of the code and partly different versions of the physics, in particular the equation of state and the opacity, as well as the option to include diffusion and settling. To make the code somewhat more user-friendly, scripts have been made which allow simply to change a few key parameters, such as mass, heavy-element abundance and number of timesteps, by editing templates of the input files. Consequently, the code has been used with success by several users, including students and postdocs at the University of Aarhus.
In addition to summary output files listing global properties of the models in the evolution sequence, output of detailed models, on the full mesh of the calculation, can be made in three different forms: the so-called emdl files, including just the variables {z i } actually used in the calculation, as well as a complete listing of the input parameters, to provide full documentation of the calculation; the amdl files which provide the variables needed for the Aarhus adiabatic pulsation code (see Christensen-Dalsgaard, 2007a); and the gong files, giving an extensive set of variables for use in further investigations of the models or plotting. A convenient way to use the code, without overloading storage systems with the large gong files, is to store the full emdl file and subsequently read in models at selected timesteps, to output the corresponding amdl or gong files.
For asteroseismology it is evidently crucial to compute oscillation frequencies of the computed models. The full calculation of frequencies, for given input parameters to the evolution calculation, often needs to be carried out as part of a larger computation, e.g., when fitting observed frequencies to determine the properties of the model. To facilitate this a version of the code has been made where the evolution calculation and all parts of the adiabatic oscillation calculation is carried out by a single subroutine call, with internal passing of the intermediate products of the calculation. This subroutine can then be called by, for example, a fitting code. An example of such use is the application of the code in genetic-algorithm fitting (e.g., Metcalfe and Charbonneau, 2003), under development by T. Metcalfe, High Altitude Observatory.
The code has been implemented on a variety of platforms and appears to be relatively robust. To simplify the installation, a complete tar package including all the required files, with a setup script and a full makefile, has been established. However, the complexity of the code and the lack of adequate documentation makes it unrealistic to release it for general use.

Further developments
From the preceding presentation it is obvious that there are significant deficiencies in ASTEC, and work is ongoing to correct them. A fairly trivial issue is the restricted nuclear network and the failure to include all elements undergoing nuclear reactions in the full diffusive treatment. Rather more serious problems concern the treatment of the borders of convective regions; even though this obviously also involves open issues of a basic physical nature, the code should at least aim at treating these regions in a numerically consistent, even if perhaps not physically adequate, manner. A serious problem is the failure of the code for models with convective cores, when diffusion and settling of both helium and heavy elements are included (cf. Christensen-Dalsgaard, 2007b); on the other hand, the case of just helium diffusion can be handled. This problem may be related to issues of semiconvection where convective stability is closely related to the details of the composition profile .
In ASTEC I have considered two cases of semiconvection, although the implementation is still under development. One case concerns growing convective cores in nondiffusing models (see also Sect. 4), as illustrated in Fig. 1. The dashed curves illustrate the properties of a model computed with the normal ASTEC implementation, where the hydrogen abundance X is discontinuous at the boundary of the core (see panels a and b). As discussed above, the extent of the convective core is determined by the behaviour of the radiative temperature gradient ∇ rad in the radiative region Fig. 1 Hydrogen-abundance profile and ∇ rad − ∇ ad against fractional radius in two 1.5M ⊙ models of age 1.36 Gyr and with Z = 0.02. The dashed curves illustrate the usual treatment in ASTEC, where the boundary of the convective core is determined by the composition in the radiative region just outside it, and the hydrogen abundance X is discontinuous. In the model illustrated by the solid curves a gradient in X is set up such that ∇ rad − ∇ ad = 0 in the outer parts of the convective core.
(see panel c), 3 leading to a convectively stable region just beneath the boundary which nonetheless is assumed to be fully mixed. A possibly more reasonable treatment, illustrated by the continuous curves, assumes that a hydrogen-abundance profile is established such that ∇ rad = ∇ ad in the outermost parts of the convective core. Since this affects a very small part of the core of the model the effects on its global properties are modest; however, it may have some influence on its pulsational properties, particularly for g modes or interface modes predominantly trapped near the boundary of the core, possibly offering the potential for asteroseismic diagnostics of this behaviour. Significant sensitivity to the treatment of this region was indeed found by Moya et al. (2007); a similar model, but computed without such treatment of semicon-vection (and with somewhat inadequate mesh resolution of the critical region) displayed far larger frequency differences than did the model considered in Fig. 1 between oscillation codes using different numerical techniques, particularly for modes with g-mode like behaviour and partly trapped near the edge of the core.
The second aspect of semiconvection concerns the base of the convective envelope in models with diffusion and settling of helium and heavy elements. As noted by Bahcall et al. (2001) and discussed in more detail by Christensen-Dalsgaard and Di Mauro (2007) the gradient in the heavy-element abundance Z established by settling beneath the convection zone leads to convective instability, e.g., in 1M ⊙ models of age somewhat higher than the present Sun. Complete mixing of the unstable region removes the gradient and hence the cause of the instability, leading to an uncertain situation, characteristic of semiconvection. As above, a perhaps reasonable assumption might be that partial mixing takes place to establish composition gradients that ensure neutral convective stability in the critical region, with ∇ rad = ∇ ad . Since the opacity depends both on X and Z the mixing must consistently affect both abundances. I have attempted to implement this by including a turbulent diffusivity, obviously common to all elements, determined iteratively as a function of depth beneath the convection zone such that the resulting profiles of X and Z lead to neutral stability.

Concluding remarks
No stellar evolution code is probably ever finished or fully tested. ASTEC has certainly proved useful in a number of applications, and the results for the Sun, as applied to helioseismology, are perhaps reasonably reliable, at least within the framework of 'standard solar modelling'. It is obvious, however, that application to the increasingly accurate and detailed asteroseismic data will require further development.
For a trace element k, with charge Z k and atomic mass A k , D k = (4πρr 2 ) 2 2BT 5/2 5 1/2 ρZ 2 where Here C k H is the collision integral between element k and hydrogen, which Michaud and Proffitt (1993) approximate, fitting numerically determined values, as where Λ k H = Λ HHe − log Z k + log 2 ; similarly the collision integral C k H between element k and helium is evaluated as C k He = 1 1.2 log[exp(1.2 logΛ k He ) + 1] , where Λ k He = Λ HHe − log Z k . Also, is the reduced mass in the element k, hydrogen system, A H being the atomic mass of hydrogen; the reduced mass A k,He involving helium is defined similarly. As presented by Michaud and Proffitt (1993) the diffusion velocity of trace elements depends on the gradient in the hydrogen abundance. To incorporate this in the formalism described by Eqs 29 and 30, I express ∂ X/∂ r in terms ofỸ 1 and write the coefficient V k as