A&A 416, 1023-1036 (2004)
DOI: 10.1051/0004-6361:20031735

Diffusion in stellar interiors: Critical tests of three numerical methods

G. Meynet¹ - A. Maeder¹ - N. Mowlavi^1,2

1 - Geneva Observatory, 1290 Sauverny, Switzerland
2 - Integral Science Data Center, Ch. d'Ecogia, 1290 Versoix, Switzerland

Received 5 May 2003 / Accepted 18 November 2003

Abstract
We describe and discuss the properties of three numerical methods for solving the diffusion equation for the transport of the chemical species and of the angular momentum in stellar interiors. We study through numerical experiments both their accuracy and their ability to provide physical solutions. On the basis of new tests and analyses applied to the stellar astrophysical context, we show that the most robust method to follow the secular evolution is the implicit finite differences method. The importance of correctly estimating the diffusion coefficient between mesh points is emphasized and a procedure for estimating the average diffusion coefficient between a convective and a radiative zone is described.

Key words: diffusion - methods: numerical - stars: interiors

1 Introduction

For many physical processes in stellar interiors, such as heat transfer, transport of angular momentum, mixing of chemical elements through time dependent convection, semiconvection or turbulence, a diffusion equation needs to be solved. Diffusion is thus involved in numerous interesting astrophysical problems. For instance, Vauclair & Charbonnel (1998) discuss the importance of selective diffusion in low-metallicity stars and the consequences for lithium primordial abundance; Charbonnel (1995) has shown how rotational diffusion in low mass stars may lead to the destruction of ³He and thus to new insights on the evolution of the abundance of this cosmological element; most of the results on the effects of rotation in stellar interiors are based on the resolution of a diffusion equation for both the transport of the chemical species and of the angular momentum (Chaboyer & Zahn 1992; Zahn 1992; Talon & Zahn 1997; Denissenkov et al. 1999; Heger et al. 2000; Maeder & Meynet 2001; Meynet & Maeder 2003). Diffusion plays also a key role in our understanding of the structure and the cooling of white dwarfs (cf. Kawaler 1995).

For solving numerically the diffusion equation, different methods are available. Quite generally, they can be classed into two main categories: the finite differences methods and the finite elements methods. Depending on how the time discretisation is performed, one distinguishes three subclasses in each of these two main categories, namely the subclasses of the explicit, implicit and of the "Crank-Nicholson'' type methods. Thus, at least six different methods are available. Here we are interested in modelling the long-term evolution of stars (secular evolution). This immediately prevents us from using explicit methods which requires the use of much too short time steps (see below). Among the four remaining methods, the "Crank-Nicholson'' finite elements method was not tested in view of the results obtained by the "Crank-Nicholson'' finite differences method (see Sect. 5.4). Thus we shall focus our attention on three of them namely

• the implicit finite elements method;

• the "Crank-Nicholson'' finite differences method;

• the implicit finite differences method.

Descriptions of these methods as well as a general discussion of their respective advantages and weaknesses can be found for instance in Press et al. (1992) for the finite differences methods and in Zienkiewicz & Taylor (2000) for the finite elements methods. If these considerations can help in choosing the best method in general, we think that very valuable complementary information can be gained by applying the methods to the resolution of realistic astrophysical cases, where the different timescales involved are very specific. This is what we propose here and this is the main aim of this work. As said in Press et al. (1992) "... differencing partial differential equations is an art as much as a science'', and a numerical scheme, which from a theoretical point of view does appear as very robust, can show some weaknesses when applied to a particular real situation. This is precisely why it is necessary to test realistic astrophysical situations. This allows us to make a new and appropriate analysis to support the choice of the best method to handle diffusion problems.

The interested reader will also find in this paper a detailed description of the discretized equations for resolving the diffusion of chemical elements and of the angular momentum in stars. We also propose a new way to estimate the diffusion coefficient in the interface between a convective and a radiative zone.

In Sect. 2 we briefly recall the physical problem to be resolved. The three numerical methods are described in details in Sect. 3. The estimate of the diffusion coefficient at the interface between a convective and a radiative zone is presented in Sect. 4. Section 5 is devoted to various tests and comparisons. Section 6 summarizes the main results.

2 The physical problem

Diffusion is a process by which components of a mixture move from one part of a system to another as a result of random motion. Let us briefly recall how the diffusion velocity is defined (see e.g. Battaner 1996). In a multicomponent fluid, one defines the mean velocity $\vec{v}_0$ of the mixture as

$$\vec{v}_0={\sum_{i} m_i n_i < \vec{v}_i>\over \sum_{i} m_i n_i},$$

where

$$<v_i> ~={1 \over n_i} \int_p f_i \vec{v}_i {\rm d}\tau_p,$$

is the mean velocity of the ith component, n_iis the partial number density of particles i, m_i their mass, v_itheir velocity and $f_i(\vec{x}, \vec{p}, t)$ the probability distribution function, which gives the probability to find a particle i at the position $\vec{x}$, with the momentum $\vec{p}$, at the time t; $\tau_p$ is the phase space volume for the particle's momentum. One can wonder, why in the expression for $\vec{v}_0$, the weighting factor is the partial mass density and not the partial number density. The question could be formulated in another way: when is the mean velocity of the mixture zero? When the net flux of the number of particles is zero ( $\sum_{i} n_i <\vec{v}_i>~=0$), or when the net flux of the mass is zero ( $\sum_{i} m_i n_i <\vec{v}_i>~=0$)? Obviously when the net flux of the mass is zero, otherwise the centre of gravity of the system is moving and thus the mean velocity of the mixture is not zero.

The peculiar velocity of a particle is defined as

$$\vec{V}_i= \vec{v}_i- \vec{v}_0.$$

The mean value of $\vec{V}_i$ is

$$< \vec{V}_i>~={1 \over n_i}\int_p f_i \vec{V}_i {\rm d}\tau_p=~< \vec{v}_i>~-~ \vec{v}_0.$$ (1)

This quantity is called the diffusion velocity of the ith component. It can be shown very easily, using the expression for $\vec{v}_0$ that the different diffusion velocities must compensate one another, i.e.

$$\sum_i m_i n_i < \vec{V}_i>~=\vec{0}.$$ (2)

There is thus no net mass flux associated with diffusion. In that respect, in a star, the changes of chemical composition due to diffusion can be treated very similarly as those due to the nuclear reactions. Like the nuclear reactions, diffusion modifies the mean molecular weight (and thus the hydrostatic structure of the star), but does not induce any net mass flux.

In the context of the Boltzmann microscopic theory, it is possible to deduce from the equations of motion for the different components, expressions for the  $<\vec{V}_i>$'s (see Battaner 1996). We shall not repeat these developments here. Instead we consider that diffusion is equivalent to a displacement of particles whose velocities may be related to the diffusive coefficient, D, by the expression (see also Sect. 3.2)

$$< \vec{V}_i>~=-{D \over X_i}\vec \nabla X_i,$$

where X_i is the mass fraction of element i. The minus sign results from the fact that the velocity is directed in the opposite direction to the mass fraction gradient. This expression automatically satisfies the condition seen above that $\sum_i m_i n_i < \vec{V}_i>~= \vec{0}$, indeed

$$\sum_i m_i n_i < \vec{V}_i>~=-\rho D \vec \nabla \left( \sum_i X_i\right)= \vec{0}$$

since $\sum_i X_i=1$ and the relation between n_iand X_i is

$$n_i={\rho X_i \over A_i m_{\rm H}}$$

where $\rho$ is the density, A_i the atomic weight expressed in units of proton mass $m_{\rm H}$.

We are interested in resolving the diffusion equation in spherically symmetric systems. Let us consider a one-dimensional problem, where particles i may diffuse along the radial direction r. Thus

$$<V_i> ~=-{D \over X_i} {\partial X_i \over \partial r}\cdot$$ (3)

The variation of the number of particles in the element of volume $\upsilon$and of surface S, due only to diffusion, is

$${\partial \over \partial t} \int_{\upsilon} n_i {\rm d}\upsilon =-\int_{S} n_i <\vec{V}_i>\cdot ~{\rm d}{\vec{S}}.$$ (4)

From the divergence theorem and the expression of <V_i> just seen above one obtains

$${\partial n_i \over \partial t}= -{1 \over r^2} {\partial \ov... ...al \left({n_i\over \rho}\right) \over \partial r}\right )\cdot$$ (5)

Replacing n_i by its expression as a function of X_iin Eq. (5) gives

$${\partial \over \partial t}(\rho X_i )={1 \over r^2} {\partia... ...r}\left(\rho r^2 D {\partial X_i \over \partial r}\right)\cdot$$ (6)

Thus Eq. (6) results from the expressions for the <V_i> (Eq. (3)) and from the continuity equation for the n_i's (Eq. (4)).

Sometimes, Eq. (6) is written as below

$${\partial \over \partial t}(\rho c_i )={1 \over r^2} {\partia... ...ial r}\left(\rho r^2 D {\partial c_i \over \partial r}\right),$$ (7)

with partial concentration c_i instead of mass fraction, without indication on the precise meaning of what concentration means. Partial concentration in mass is equivalent to the mass fraction, while partial concentration in number is equal to n_i/n where $n=\sum_i n_i$. Replacing c_i in Eq. (7), by n_i/n and expressing n_i as a function of X_i, one obtains

$${\partial \over \partial t}(\rho \mu X_i )={1 \over r^2} {\pa... ...left (\rho r^2 D {\partial \over \partial r}(\mu X_i)\right ),$$ (8)

where $\mu=\rho/(n m_{\rm H})$ is the mean molecular weight of the ions. This equation is identical to Eq. (6) only when $\mu$remains constant as a function time and is constant as a function of r. Thus Eq. (8) is equivalent to Eq. (6) only when applied to minor constituents, the abundances of which do not affect $\mu$ and when $\mu$ has no gradient. In more general cases, Eq. (8) is not equivalent to Eqs. (5) and (6).

This can also be seen in a slightly different way. When, in Eq. (7), the c_i's are identified with the number fractions, one obtains

$${\partial \over \partial t}\left( \rho {n_i \over n}\right ) ... ...al \over \partial r} \left(r^2 \rho {n_i \over n} w_i \right),$$ (9)

where we have introduced a velocity $w_i=-{D \over n_i/n} {\partial \over \partial r} \left ( {n_i \over n}\right )$. This velocity is not a diffusive velocity in the sense that $\sum_i m_i n_i w_i=-nD m_{\rm H} \partial \mu/\partial r$is not equal to zero except in zones where $\mu$ is constant. Equation (9) expresses the change of $\rho n_i/n=\mu n_i m_{\rm H}$ in an element of volume resulting from the transport of the quantity $\rho n_i/n$. But physically, this is not $\rho n_i/n$ which diffuses but the particles themselves. As above, there are two conditions for Eq. (9) to be equivalent to Eqs. (5) and (6): 1) the particles which diffuse must have a very small abundance, so small that their diffusion does not affect $\mu=\rho/(n m_{\rm H})$; 2) and $\mu$ does not vary with the radius. We conclude thus that when the diffusive transport is described by an equation of the form given in Eq. (7), the concentrations are equivalent to mass fractions. The identification with number fractions results in unphysical description of the process except in very particular situations.

Starting from Eq. (6) and summing over all the chemical species i, one obtains

$${\partial \rho \over \partial t}={1 \over r^2} {\partial \ove... ...al r}\left(\rho r^2 D {\partial \over \partial r}(1) \right)=0.$$

This is quite consistent with the fact that the diffusive velocities must compensate each other. Thus the density can be put outside from the time derivative in Eq. (6):

$$\rho{\partial X_i \over \partial t}\bigg\vert _{m_r} = {1 \ov... ... r} \left (\rho r^2 D {\partial X_i \over \partial r}\right ),$$ (10)

where m_r is the Langrangian mass coordinate. This is the equation we shall numerically resolve in this paper. The conditions at the centre and the surface are:

$${\partial X_i \over \partial r }\bigg\vert _{m_r=0}={\partial X_i \over \partial r }\bigg\vert _{m_r=M}=0,$$ (11)

where M is the total mass of the star.

The equation expressing the diffusion of the angular momentum is (see e.g. Endal & Sofia 1978)

$$\rho {\partial (r^2 \Omega) \over \partial t}\bigg\vert _{m_r... ...eft (\rho r^4 {D'} {\partial \Omega \over \partial r}\right ),$$ (12)

where $\Omega$ is the angular velocity and D' the diffusion coefficient for the angular momentum. In our numerical experiments, the radii do not change with time. Equation (12) then becomes

$$\rho r^2{\partial \Omega \over \partial t}\bigg\vert _{m_r} =... ...t )={\rm div}\left(\rho r^2 {D'} \vec \nabla (\Omega)\right ).$$ (13)

The conditions at the centre and at the surface are

$${\partial \Omega \over \partial r }\bigg\vert _{m_r=0}={\partial \Omega \over \partial r }\bigg\vert _{m_r=M}=0.$$ (14)

In general the transport equation for the angular momentum may contain other terms expressing the advection of angular momentum by meridional circulation and the effects of magnetic braking. As these last two terms bring no specific problems with respect to diffusion, we do not consider them here. However, the advection of angular momentum needs a particularly careful treatment as well (cf. Meynet & Maeder 2000).

3 The three tested methods

3.1 The implicit finite elements method

The detailed description of the implicit finite elements method used to resolve the diffusion equation for the chemical species is presented in Schatzman et al. (1981) (see the appendix by Glowinsky & Angrand). We present here the procedure for the case of the angular momentum diffusion. The basic idea of this method is to decompose the unknown function, here  $\Omega(m_r,~t)$, as a linear combination of well chosen independent functions. This decomposition can be performed in many different ways. We adopt here the same decomposition as in Schatzman et al. (1981). Of course, the conclusions concerning the ability of this method to provide physical solutions will only refer to this particular choice.

We decompose the star in K shells, with the langrangian mass coordinate of the ith mesh point being m_i (m_i is the mass inside the sphere of radius r_i). In the following all the quantities with an indice i are evaluated at the mesh point i. The shells are numbered from 1 at the surface to K at the centre. Let us introduce K functions b_i(m_r) defined by

$$\begin{array}{clll} & &{m_r-m_{i-1} \over m_i-m_{i-1}} & {\r... ...& 0 & {\rm if}\ \ m_r \notin [m_{i+1},m_{i-1}]. \\ \end{array}$$

The function b_i is equal to one at m_r=m_i, is equal to zero at m_i+1 and m_i-1 and varies linearly as a function of m_r inbetween. By multiplying Eq. (13) by each of the functions b_i(m_r), one obtains K equations. The integration over the whole volume of the star, V, of each of these Kequations gives

$$\int_V \rho r^2 {\partial \Omega \over \partial t} b_i {\rm d... ...div}\left(\rho r^2 {D'}\vec{grad}(\Omega) \right)b_i {\rm d}V.$$ (15)

Using the general relations:

$${\rm div}(a\vec{v})=a {\rm div}(\vec{v})+\vec{grad}(a)\cdot\vec{v},$$ (16)

$$\int_V {\rm div}(\vec{v}){\rm d}V=\int_S \vec{v}\cdot{\rm d}\vec{S},$$ (17)

where a is a scalar, $\vec{v}$ a vector and S the surface of the volume V, one obtains that the right hand term of Eq. (15) becomes

$$\int_S \rho r^2 {D'} \vec{grad}(\Omega) b_i\ {\rm d}\vec{S}- ... ...ho r^2 {D'} \vec{grad}(\Omega)\cdot \vec{grad}(b_i)\ {\rm d}V.$$

The integral on the surface is null, since $\vec{grad}(\Omega)={0}$on the surface. Thus one has K equations of the type

$$\int_V \rho r^2 {\partial \Omega \over \partial t} b_i {\rm d... ...\over \partial r} {\partial b_i \over \partial r} {\rm d}V=0.$$ (18)

Using ${\rm d}m_r=\rho {\rm d}V=4\pi r^2 \rho {\rm d}r$, one obtains

$$\int_0^M \left[ r^2 {\partial \Omega \over \partial t} b_i +r... ... m_r} {\partial b_i \over \partial m_r}\right ] {\rm d}m_r=0.$$ (19)

The functions b_i constitute a set of independent functions, therefore the unknown function $\Omega(m_r, t)$ can be expressed as a linear combination of b_i's. One can write

$$\Omega(m_r, ~t)=\sum_{j=1}^{K}\Omega_j b_j(m_r),$$ (20)

where $\Omega_j=\Omega(m_j, t)$. Equation (20) simply says that to obtain $\Omega(m_r, t)$, with $m_r \in [m_{i+1},m_i]$, one simply interpolates linearly as a function of mass between  $\Omega_{i+1}$ and $\Omega_i$.

Expressing $\Omega(m_r, t)$ as indicated in Eq. (20), Eq. (19) becomes

$$\sum_j A_{ij}{\partial \Omega_j \over \partial t} +\sum_j B_{ij} \Omega_j =0,$$ (21)

where

$$A_{ij}=\int_0^M r^2 b_i b_j\ {\rm d}m_r,$$ (22)

$$B_{ij}=\int_0^M r^2 {D'} \left(4\pi r^2 \rho\right)^2 {\parti... ... \partial m_r} {\partial b_j \over \partial m_r}\ {\rm d}m_r.$$ (23)

We have used the facts that the $\Omega_j$s do not depend on m_rand the b_js do not depend on time. When |j-i| > 1, there is no overlap between the functions b_i and b_j, thus

$$A_{ij}=0,\ \ \ {\rm for}\ \ \ \vert j-i\vert > 1.$$ (24)

For a given value of i only, A_i, i-1, A_i, i and A_i+1, i are different from zero. The same is true for the B_ij. Let us estimate  A_i, i-1:

$$A_{i,~i-1}=\int_{m_{i}}^{m_{i-1}} r^2 {m_r-m_i \over m_{i-1}-m_i}{m_r-m_{i-1} \over m_{i}-m_{i-1}}\ {\rm d}m_r,$$ (25)

one approximates r² by $0.5\left(r^2(m_{i-1})+r^2(m_i)\right)$, then, by trivial integration, one obtains

$$A_{i,i-1}=-0.5\left ( r^2(m_{i-1})+r^2(m_i)\right ){m_i-m_{i-1}\over 6}\cdot$$ (26)

The other matrix elements are obtained in the same way:

$$A_{i+1,~i}=-0.5\left ( r^2(m_{i+1})+r^2(m_i)\right ){m_{i+1}-m_i\over 6},$$ (27)

A_i, i=2(A_i, i-1+A_i+1, i). (28)

At the centre and at the surface, one has

$$A_{K,~K} =0.5 r^2(m_{K-1}){m_{K-1}\over 3},$$ A_K, K-1=0.5 A_K, K, A_2, 1 =0.5 A_1,1,

$$A_{1,~1} =0.5\left ( r^2(1)+r^2(2)\right ){m_1-m_2 \over 3}\cdot$$

Let us set $\overline D=r^2 (4 \pi r^2 \rho)^2D'$ and $\overline D_{i,~i-1}$ an appropriate mean value of $\overline D$ between the shells i and i-1(see Sect. 4), then the B_i, jbecomes

$$B_{i,~i-1}={\overline D_{i,~i-1}\over m_i-m_{i-1}},$$ (29)

$$B_{i+1,~i}={\overline D_{i+1,~i}\over m_{i+1}-m_i},$$ (30)

B_i, i=-B_i, i-1-B_i+1, i. (31)

At the centre and at the surface, one obtains

$$\begin{array}{ll} B_{K,~K} &={\overline D_{K,~K-1}\over m_{K... ... B_{1,~1} &={\overline D_{2,1}\over m_1-m_2}\cdot \end{array}$$

Adopting an implicit discretisation in time, we have for Eq. (21)

$$\sum_j A_{i,~j}\left({\Omega_j^{n+1}-\Omega_j^n \over \Delta t} \right) +\sum_j B_{i,~j}\Omega_j^{n+1}=0,$$ (32)

or

$$\sum_j \left({A_{i,j}\over \Delta t}+B_{i,~j} \right)\Omega_j^{n+1}=\sum_j {A_{i,~j}\over \Delta t}\Omega_j^n,$$ (33)

where $\Omega_j^n$ is equal to $\Omega(m_j,t^n)$. This linear system of equations is then solved by standard procedure.

3.2 The Crank-Nicholson finite differences method

Let us consider in a first step the case of the diffusion of the chemical elements. Multiplying Eq. (10) by $4\pi r^2$ and integrating with respect to the spatial coordinate r from $r-\Delta r/2$ to $r+\Delta r/2$, one obtains

$$\Delta m {\partial X_i \over \partial t}\bigg\vert _{r}=\left... ...\partial r}\right ) \bigg\vert _{r-\Delta r/2}^{r+\Delta r/2},$$ (34)

where $\Delta m=\int_{r-\Delta r/2}^{r+\Delta r/2} 4 \pi r^2 \rho dr$. The quantity $\Delta r$ has been chosen sufficiently small for ${\partial X_i \over \partial t}$ to remain constant over the spatial range $\Delta r$ around r. This integration reduces the equation to a first order differential equation, which simply expresses the fact that the time variation of the mass of element i in a spherical shell is equal to the difference between the mass of element i which enters in the shell and the mass of this same element which goes out from the shell. As indicated in Sect. 2, one can interpret the quantity,

$$-{D \over X_i} {\partial X_i \over \partial r}=-D {\partial \ln X_i \over \partial r}=~<V_i>,$$ (35)

as a diffusion velocity, <V_i>, for the element i (let us recall that the minus sign appears because, when the spatial abundance gradient is positive, the velocity is directed inward). In that case, Eq. (34) writes

$$\Delta m {\partial X_i \over \partial t}\bigg\vert _r=-4\pi \rho r^2 X_i <V_i>\vert _{r-\Delta r/2}^{r+\Delta r/2},$$

which makes clearly appear the difference of the mass fluxes between the two borders of the shell. This is a more elementary description of diffusion.

One can also derive Eq. (35) from the continuity equation. Indeed supposing that the bulk gas velocity is zero, that $\partial \rho / \partial t =0$(stationary situation), and that there is no sink/source terms of matter, the continuity equation becomes

$$\rho{\partial X_i \over \partial t}\bigg\vert _{r} =- {1 \ove... ... {\partial \over \partial r} \left( \rho r^2 <V_i> X_i\right).$$ (36)

By equating the right hand sides of Eqs. (10) and (36), one obtains for <V_i> Eq. (35) above.

	Figure 1: Discretized distribution of the abundance of element i in a star model as a function of radius. Only a few mesh points are represented. K is the total number of mesh points. Numbering increases from outside to inside.
Open with DEXTER

Let us now consider the discretized abundance profile of element i sketched in Fig. 1. The mass of element i removed from point k+1 and added to point k per unit time is (see Eq. (34)):

$$-4 \pi r^2_{k+1/2} \rho_{k+1/2} D_{k+1/2} {\partial X_i \over \partial r}\bigg \vert _{k+1/2},$$ (37)

where the physical quantities with a lower subscript k+1/2 are estimated at the midpoint between the mesh points k+1 and k. For instance, one takes r_k+1/2 =0.5(r_k+r_k+1). For the diffusion coefficient, we use the procedure exposed in Sect. 4. A similar expression can be written for the mass of element iremoved from point k and added to point k-1. Let us set X_i, k, the mass fraction of element i evaluated at the mesh point k. Thus one has,

 

$$\Delta m_k {\partial (X_{i,~k}) \over \partial t} - 4 \pi r^2_{k+1/2} \rho_{k+1/2} D_{k+1/2} {\partial X_i \over \partial r}\bigg \vert _{k+1/2} +4 \pi r^2_{k-1/2} \rho_{k-1/2} D_{k-1/2} {\partial X_i \over \partial r}\bigg \vert _{k-1/2},$$ (38)

where $\Delta m_k=\int_{r_{k+1/2}}^{r_{k-1/2}} 4 \pi r^2 \rho {\rm d}r.$Let us now discretize this equation as a function of time. The left hand term of Eq. (38) can be written

$$\Delta m_k {\partial (X_{i,~k}) \over \partial t}\rightarrow \Delta m_k {X_{i,~k}^{n+1}-X_{i,~k}^{n} \over \Delta t},$$ (39)

where the superscript n indicates that the quantity is estimated at the time tⁿ. The Crank-Nicholson method consists in evaluating the right hand term at the same time as the left hand term, i.e. at time t^n+1/2. The system is centered in time and thus second order accurate in time. One obtains

 

$${X_{i,~k}^{n+1}-X_{i,~k}^{n} \over \Delta t} -{\sigma_{k+1/2} \over \Delta m_{k}} D_{k+1/2} {{X_{i,~k+1}^{n+1}... ...{i,~k+1}^{n}\over 2}-{ X_{i,~k}^{n+1}+ X_{i,~k}^{n}\over 2}\over r_{k+1}-r_{k}}$$ =

$$+{\sigma_{k-1/2} \over \Delta m_{k}} D_{k-1/2} {{X_{i,~k}^{n+1}+ ... ...,~k}^{n}\over 2}-{ X_{i,k-1}^{n+1}+ X_{i,~k-1}^{n}\over 2}\over r_{k}-r_{k-1}},$$ (40)

where one has introduced the quantity $\sigma_{k+1/2}$ defined by

$$\sigma_{k+1/2}=4\pi r^2_{k+1/2}\rho_{k+1/2}.$$ (41)

Let us define the quantity [k+1,k] by

$$[k+1,k]={1 \over 2} \sigma_{k+1/2} {D_{k+1/2} \Delta t \over r_{k+1}-r_k},$$ (42)

and let us separate the abundances evaluated at time tⁿ from those estimated at time tⁿ⁺¹. Equation (40) then becomes

 

$${[k,k-1] \over \Delta m_k} X^{n+1}_{i,~k-1}+ \left (1- {[k,~k-1] ... ...\Delta m_k}\right )X_{i,~k}^{n+1} +{[k+1,~k] \over \Delta m_k}X_{i,~k+1}^{n+1}=$$

$$-{[k,~k-1] \over \Delta m_k} X^{n}_{i,~k-1}+ \left (1+ {[k,~k-1] ... ...er \Delta m_k}\right )X_{i,~k}^{n} -{[k+1,~k] \over \Delta m_k}X_{i,~k+1}^{n} .$$ (43)

Following a similar line of reasoning, one obtains for the equations at the center and at the surface (see Fig. 1),

 

$${[K,K-1] \over \Delta m_K }X_{i,K-1}^{n+1}+ \left (1-{[K,K-1] \over \Delta m_K }\right )X_{i,K}^{n+1}=$$

$$\quad\quad\quad\quad-{[K,K-1] \over \Delta m_K }X_{i,K-1}^{n}+ \left (1+{[K,K-1] \over \Delta m_K }\right )X_{i,K}^{n} ,$$ (44)

 

$${[2,1] \over \Delta m_1 }X_{i,2}^{n+1}+ \left (1-{[2,1] \over \Delta m_1 }\right )X_{i,1}^{n+1}=$$

$$\quad\quad\quad\quad\quad\quad-{[2,1] \over \Delta m_1 }X_{i,2}^{n}+ \left (1+{[2,1] \over \Delta m_1 }\right )X_{i,1}^{n} ,$$ (45)

where $\Delta m_K=\int_{0}^{r_{K-1/2}} 4 \pi r^2 \rho {\rm d}r$ and $\Delta m_1=\int_{r_{1-1/2}}^{r_1} 4 \pi r^2 \rho {\rm d}r$. Equations (43-45) constitute a system of linear equations whose unknowns are the  X_i, kⁿ⁺¹. It is solved by using classical methods of tridiagonal matrix inversion.

The discretized equations describing the diffusion of the angular momentum are obtained in a similar way. The final result are equations similar to Eqs. (43)-(45) with X_i replaced by $\Omega$, $\Delta m_k$ replaced by $\Delta m_k r^2_k$, $\sigma_{k+1/2}$ replaced by $\sigma_{k+1/2} r^2_{k+1/2}$ and  $\sigma_{k-1/2}$ replaced by $\sigma_{k-1/2} r^2_{k-1/2}$. Of course the diffusion coefficient must be the one describing the diffusion of the angular momentum.

This method is not fully implicit and thus, one can expect that some sort of Courant's condition will limit its domain of validity (see Sect. 5).

3.3 The implicit finite differences method

In the implicit finite differences method, the right hand term of Eq. (38) is estimated at time tⁿ⁺¹ (in the explicit method the right hand term would be estimated at time tⁿ). In this case, one obtains

 

$${X_{i,~k}^{n+1}-X_{i,~k}^{n} \over \Delta t} -{\sigma_{k+1/2} \over \Delta m_{k}} D_{k+1/2} {X_{i,~k+1}^{n+1}- X_{i,~k}^{n+1}\over r_{k+1}-r_{k}} +{\sigma_{k-1/2} \over \Delta m_{k}} D_{k-1/2} {X_{i,~k}^{n+1}-X_{i,~k-1}^{n+1}\over r_{k}-r_{k-1}}\cdot$$ (46)

Let us define the bracket terms [k+1,k] by

$$[k+1,k]=\sigma_{k+1/2} {D_{k+1/2} \Delta t \over r_{k+1}-r_k},$$

doing so, we obtain, separating the abundances at time tⁿ from those estimated at time tⁿ⁺¹,

$${{[k,~k-1] \over \Delta m_k} X^{n+1}_{i,~k-1}+ \left (1- {[k,~k-1] \over \Delta m_k}-{[k+1,~k] \over \Delta m_k}\right )X_{i,~k}^{n+1}}$$

$$\quad\quad+{[k+1,~k] \over \Delta m_k}X_{i,~k+1}^{n+1}= X_{i,~k}^{n} .$$ (47)

Following a similar line of reasoning, one obtains for the equations at the center and at the surface,

$${[K,K-1] \over \Delta m_K }X_{i,~K-1}^{n+1}+ \left (1-{[K,~K-1] \over \Delta m_K }\right )X_{i,~K}^{n+1}= X_{i,~K}^{n} ,$$ (48)

$${[2,1] \over \Delta m_1 }X_{i,~2}^{n+1}+ \left (1-{[2,~1] \over \Delta m_1 }\right )X_{i,~1}^{n+1}=X_{i,~1}^{n}.$$ (49)

One can easily check that the system of discretized equations (both in the cases of the Crank-Nicholson scheme and in the implicit finite differences method) conserves the integrated quantities of the elements over the mass of the star. The system also keeps equal to one the sum of the X_iand does not induce any diffusion when the chemical gradients are flat.

For obtaining the equations for the angular momentum, we apply the same recipe as indicated at the end of the previous section. The system of equations describing the transport of the angular momentum also conserves the total angular momentum.

	Figure 2: Schematic representation of the mass shell between the radii rk and rk-1. The diffusive coefficient Dk operates in the zone between rk and $r_{\rm c}$. The diffusive coefficient Dk-1 operates in the zone between $r_{\rm c}$ and rk-1. The quantity $\Delta r$ is equal to rk-1-rk.
Open with DEXTER

4 Estimate of ${{\vec D}}_{{{\vec k}}-{\vec {1/2}}}$

Let us consider two mesh points as in Fig. 2. Let us suppose that the diffusion coefficient D_k-1 is the coefficient for the whole region between r_k-1 and $r_{\rm c}$. Similarly the diffusion coefficient D_k is the coefficient for the whole region between $r_{\rm c}$ and r_k. When both mesh points are in a radiative or convective region, the radius $r_{\rm c}$ is equal to (r_k-1+r_k)/2, otherwise $r_{\rm c}$ is taken as the radius of the limit between the radiative and the convective zone.

Over the path $f\Delta r$ (see Fig. 2), the particles have an average diffusive velocity V_k(see Eq. (35)) and over the path $(1-f)\Delta r$ they have a diffusive velocity <V_k-1>. The total time for going from k to k-1is equal to

$$\Delta t = f\Delta r/<V_{k}> + ~(1-f)\Delta r/<V_{k-1}>.$$ (50)

The mean velocity <V_k-1/2> over the whole interval is given by $\Delta r/\Delta t$, thus one has

$$<V_{k-1/2}>~={1 \over {f \over <V_{k}>} + {(1-f) \over <V_{k-1}>}} ={<V_{k}> <V_{k-1}>\over f <V_{k-1}>+(1-f)<V_{k}> }\cdot$$

Now, the diffusive coefficient is proportional to the diffusive velocity (see Eq. (35)). This suggests the way to compute the diffusive coefficient between two shells:

$$D_{k-1/2}={D_{k-1} D_{k} \over f D_{k-1} + (1-f) D_{k}}\cdot$$ (51)

This expression is more physical than

D_k-1/2=(D_k-1 + D_k)/2,

that simply averages the diffusion coefficients. Equation (51) implies that if $D_{k} \gg D_{k-1}$ and (1-f) is of the order of 0.5 then $D_{k-1/2} \sim D_{k-1}$. As expected, the smallest diffusion coefficient governs the diffusion between the two mesh points. The simple algebraic mean would give $D_{k-1/2} \sim D_{k}/2$ a much greater diffusion coefficient, which is not physically justified. One sees also that when D_k-1 = D_k = D then D_k-1/2=D whatever the value of f. In Sect. 5 below we shall illustrate by a numerical example the importance of correctly evaluating the diffusion coefficient at the interface between a convective and a radiative zone.

When the finite elements method is used, we need to evaluate the mean value between two mass shells of the diffusion coefficient, D, multiplied by a factor $\alpha$, which for the transport of the chemical elements is $(4\pi r^2 \rho)^2$ and for that of the angular momentum is $r^2(4\pi r^2 \rho)^2$ (see Sect. 3.1). In that case, one adopts the following procedure:

$$\alpha_{k-1/2} D_{k-1/2}=0.5\left (\alpha_{k}+\alpha_{k+1} \right){D_{k-1} D_{k} \over f D_{k-1} + (1-f) D_{k}}\cdot$$

Figure 3: Comparison of the abundances obtained with the implicit finite differences method using two different ways of evaluating the diffusion coefficient Dk-1/2between two mesh points (see Sect. 4). The dotted line shows the initial situation ($\tau =0$). A time step $\Delta t$ of 10 years was adopted and the computation was stopped at an age of $\tau = 1000$ years. The dashed line show the result obtained when one takes for Dk-1/2 a simple algebraic mean, the continuous line presents the result obtained using Eq. (51).

Open with DEXTER

5 Tests and comparisons

5.1 Initial configuration and Courant's condition

In this section, we study how the different methods described above behave in a very simple configuration. Let us consider a uniform density sphere of 1 $M_\odot$, with a radius of 1 $R_\odot$, composed of two chemical elements X and Y. The initial distribution of these two elements in the sphere is a quasi step function. The variation as a function of the radius of the abundance of one of these two elements, let us call it X, at the beginning of the computation is shown in Fig. 3 by the dotted line (cf. also the dotted lines in Figs. 5 to 8). At the center, its initial value is 0.99, in the outer regions X=0.01. The abundance Y is simply 1-X. The diffusion coefficient is set equal to 10¹⁴ cm² s^-1in the interior region (i.e. for a radius r inferior to 0.7963 $R_\odot$ or a mass inferior to 0.5049 $M_\odot$) and to $1.5 \times 10^6$ cm² s^-1in the outer zone. Such values for the diffusion coefficient imply that the interior has a very small mixing timescale (of the order of $R^2/D \sim 1$ yr) and therefore will be always homogenized by mixing, while the outer one will have a much longer mixing timescale (of the order of 4.26 Myr). Such a situation is quite analog to a convective core surrounded by a radiative envelope in a star. In the following we shall define the interior region as "the core'', and the outer one as "the envelope''. In our numerical experiments we keep the diffusion coefficient constant as a function of time. The sphere is decomposed in 100 shells of equal mass.

From these data one expects that the whole star will be completely mixed in less than about 10⁷ yr and that the final abundance of element Xin the homogeneous star will be 0.5048. From this initial structure, on can also estimate the "Courant condition''. Let us recall that the "Courant condition'' imposes a superior limit $\Delta t_{\rm c}$to the time step for an explicit method to be stable (see e.g. Press et al. 1992). To estimate this limit, one has to compute for each element and for each mass shell, the time required for the element to diffuse through the width of the shell i.e.

$$\Delta t_X = {\Delta r \over v_X} = {\Delta r \over {D \over X}\left\vert {\partial X \over \partial r}\right\vert}\cdot$$

Then one has to take the smallest value of all these times. In a discretized form, $\Delta t_{\rm c}$ may be written

$$\Delta t_{\rm c} = {\rm Min}\left ({(r_i-r_{i+1})^2 \over {D_... ...X_{i+1} \over {(X_i+X_{i+1}) \over 2}}\right\vert}\right )\cdot$$

For the simple initial structure considered here, the diffusion time through a shell takes a non infinite value only at the interface between the core and the envelope. In this case, one has that r_i-r_i+1 is equal to 0.000199 $R_\odot$, D_i+1/2, estimated from Eq. (51) with f taken equal to 0.5, is equal to $3 \times 10^6$ cm² s, X_i+1-X_i = 0.98 and ${X_i+X_{i+1} \over 2}$ is equal to 0.5. This gives

$$\Delta t_{\rm c} \sim 1\ {\rm yr.}$$

Let us stress that implicit methods are generally stable for any time step, and thus are not limited by the Courant condition exposed above (cf. Press et al. 1992). In the following we shall test this point.

5.2 Comparisons of the methods

Let us begin by illustrating the importance of correctly estimating the diffusion coefficient at the interface between a convective and a radiative zone. In Fig. 3, the resulting distributions of the element X inside the star after 1000 yr is shown. The time step used is $\Delta t = 10$ yr The continuous line shows the results obtained using Eq. (51) for D_k-1/2. The dashed line represent the solution obtained using a simple algebraic mean. We see that the results are significantly different, in particular the algebraic mean tends to slightly increase the convective core. This may have important consequences when such increases are repeatedly applied over the whole evolution of a star. The use of Eq. (51) which results from physical considerations is thus recommended, and this is the one we used in the numerical experiments we now describe.

We have computed the diffusion of the chemical elements (and of the angular momentum, see Sect. 5.7) for different durations $\tau$ of the period during which the diffusion operates and for different time steps $\Delta t$. In Fig. 4 the set of values ( $\log \Delta t$, $\log \tau$) explored are indicated. For each couple of values, we performed computations with the three methods described above, namely the implicit finite elements method, the Crank-Nicholson finite differences method and the implicit finite differences method. Most of the results are displayed in Figs. 5 to 8. On these figures, the results obtained for increasing durations, using the same time step, are ordered horizontally (from left to right), while the results obtained for the same duration, with increasing time steps, are arranged vertically (from top to bottom). Since the Courant time step is equal to one year, the time step $\Delta t$ expressed in years gives directly the time step in units of the Courant time step.

Figure 4: Computations have been performed with the three methods for each couple of ( $\log \Delta t$, $\log \tau$) values corresponding to the position of a square in the above figure; $\Delta t$ is the time step, $\tau$ is the duration over which the computation was performed, (no computation, of course, have been performed in the grey zone where $\log \Delta t > \log \tau$). The zone I (square with a cross inside) shows the "non-physical'' region for the implicit finite elements method. By "non-physical'', we mean here that negative abundances are obtained. The zone II (square with a filled triangle inside) corresponds to the zone of "non-physical'' solutions for the Crank-Nicholson finite differences method. The implicit finite differences method gives physical solution in all the cases considered here (indicated by squares filled or empty).

Open with DEXTER

5.3 The implicit finite elements method

The solutions obtained with the implicit finite elements method are shown in Figs. 5 and 6 (the dashed-dotted lines). This method gives reasonable results in all the cases explored in this work except for small values of $\tau$. We note that for $\tau < 1000$ yr, whatever is the time step, the solution leads to negative values of X just above the "convective core''. If the integration is performed over a sufficiently long time, the system eventually reaches a reasonable solution, even when instabilities appear at earlier time (see for instance in Figs. 5 and 6 the evolution when $\tau$ increases with $\Delta t = 0.5$ yr). One notes that the star is completely mixed for $\tau > 10^7$, i.e. for durations more than twice the mixing timescale of the envelope (about 4.3 Myr, see Sect. 5.1), which is not very satisfactory. The final abundance in the homogeneous star is on the other hand equal to that expected ($\sim$0.5).

As expected from an implicit scheme, reasonable solutions are obtained even if the time steps are much greater than the time step given by the Courant condition. Inspection of the results for $\tau \geq 1000$ yr show in general a great stability of the solution with respect to the choice of the time step. The solutions obtained with the implicit finite elements method are in general less mixed that those obtained with the implicit finite differences method (see Figs. 7 and 8).

5.4 The Crank-Nicholson finite differences method

This method is a kind of compromise between a fully implicit and an explicit scheme. In order to better understand what happens here, let us briefly recall a few generalities about implicit and explicit schemes (see Press et al. 1992, pp. 838-842).

Explicit finite differences schemes are stable only if the time steps satisfy the Courant condition. Such methods are not suited for the computation of the secular evolution of stars. Indeed, we are interested in modeling the evolution of features with spatial scales of the order of the radius of the star R, much greater than the distance between two mesh points $\Delta r$. If we are limited to time steps satisfying the Courant condition, we will need of the order of $R^2/(\Delta r)^2$ steps before anything interesting begins to happen. We would like instead use timesteps of the order of R²/D or, maybe, for purpose of accuracy, somewhat smaller.

With such great timesteps, it is no longer possible to describe accurately what happens at small spatial scales. However at small scales, the differencing scheme must do "something stable, innocuous, and perhaps not too physically unreasonable'' write Press et al. (1992). One possibility is to use a Crank-Nicholson differencing scheme. One of its main property is to let small-spatial-scale features maintain their initial amplitudes. In that case, according to Press et al. (1992), the evolution of the larger-scale features, in which we are interested in, take place superposed with a kind of "frozen in'' (though fluctuating) background of small-scale features. This is what happens in our numerical experiments. Indeed looking at the continuous lines in Figs. 5 and 6, one sees that the amplitude of the initial step in chemical composition is more or less maintained, although with great fluctuations.

The method gives reasonable solutions for not too big time steps. More precisely, physical solutions are obtained when $\Delta t$ in years is inferior to $10\sqrt{\tau}$, or when $\tau$in years is inferior to 100 yr (see Fig. 4). This restriction on the time step is reminiscent of a kind of Courant's condition. This is not so surprising given that the Crank-Nicholson scheme is a mixture of both an implicit (not submitted to Courant condition) and an explicit method (submitted to Courant condition).

Let us finally note that in its domain of validity in the log $\Delta t$ versus log $\tau$ plane this method gives the same solution as the implicit finite differences method, and since, it is centered in time, this scheme is second-order accurate in time.

5.5 The implicit finite differences method

Another possibility for imposing an "inoccuous'' behaviour to the small-scale features is to adopt an implicit method. Such schemes drive small-scale features to their equilibrium form, i.e. imposes that at small scales

$${\partial \over \partial r} \left (\rho r^2 D {\partial X_i \over \partial r}\right )\rightarrow 0.$$

This can be seen from Eq. (46) with $\Delta t \rightarrow \infty$.

The solutions obtained with the implicit finite differences method are shown in Figs. 7 and 8 (continuous lines). The results show in general a great stability of the solution with respect to the choice of the time step. Only when the time step is of the same order of magnitude as $\tau$ can we notice some differences (compare for instance the continuous line for $\tau = 10^7$ yr and $\Delta t = 10^6$ yr with the case $\tau = 10^7$ yr and $\Delta t = 10^7$ yr in Fig. 8).

This method leads to a completely mixed star after a time less than 10⁷ in agreement with the estimate made in Sect. 5.1. The final abundance of the element X in the homogeneous star is also equal to that expected.

Although, in that case, the differencing scheme is only first-order accurate in time, it has the advantage over the Crank-Nicholson finite differences method and the implicit finite elements method to propose a physical solution (i.e. without negative abundance values) for all the cases investigated here (see Fig. 4). Moreover, as seen just above, it predicts a timescale for the star to be completely homogenized in agreement with the usual analytical estimate (see Sect. 5.1). In that respect this method does appear as the best one.

Table 1: Maximum value over the star of the quantity $\Delta \chi =1-X(k)-Y(k)$ and of $\Delta B={B_{\rm final}-B_{\rm initial} \over B_{\rm initial}},$for different values of $\tau$, $\Delta t$ (in years) and for different numerical schemes. X and Y are the mass fraction of the two elements composing the "star'', B is the total angular momentum of the star. The labels CN and FI are for "Crank-Nicholson'' and "Fully Implicit'' respectively (see text).

Figure 5: Abundance of the element X as a function of the distance to the center. The dotted line show the situation at the beginning of the computation. The continuous line refers to the solution obtained after a time $\tau$ by the Crank-Nicholson finite differences method. The dashed-dotted line shows the result obtained after a time $\tau$ by the implicit finite elements method. The time step $\Delta t$used in each case is indicated. In all the cases, the Courant condition time is about 1 yr. This means also that the value of $\Delta t$ is equal to the ratio $\Delta t/\Delta t_{\rm c}$.

Open with DEXTER

	Figure 6: Same as Fig. 5 for other values of $\tau$ and $\Delta t$.
Open with DEXTER

	Figure 7: Same as Fig. 5 except that the continuous line corresponds to the result obtained after a time $\tau$ by the implicit finite differences method.
Open with DEXTER

	Figure 8: Same as Fig. 6 except that the continuous line corresponds to the result obtained after a time $\tau$ by the implicit finite differences method.
Open with DEXTER

5.6 Conservation of the sum of the abundances

Let us now have a look on the ability of the different schemes to conserve the sum of the abundances (in mass fraction). In each shell and at any time step, one should have that X+Y=1. In Table 1, we indicate for all the cases where the three methods give physical solutions, the maximum over the star of the quantity $\Delta \chi =1-X(k)-Y(k)$ at the end of the computation. As expected we see that in general  $\Delta \chi$ becomes greater when $\tau$ increases. The Crank-Nicholson finite differences method appears to give, in most of the cases, the smallest values for  $\Delta \chi$ (inferior to 10^-4 for the cases considered). This is likely related to the fact that this method is second order accurate in time. The implicit finite elements method gives in general the highest values (in the worst case  $\Delta \chi$ is of the order of 10^-3), while the implicit finite differences method gives in general results inbetween. Thus the implicit finite differences method enables to avoid unphysical solutions and keeps reasonably well the sum of the abundances equal to one.

5.7 Conservation of the angular momentum

We have performed similar tests and comparisons for the diffusion of the angular momentum. We started from an initial configuration where the "convective core'' defined in Sect. 5.1 has an angular velocity equal to $1\times 10^{-5}$ s^-1 and the radiative envelope an angular velocity equal to $1\times 10^{-6}$ s^-1. The other variables were taken as described in Sect. 5.1. In such situation the Courant time step, which writes

$$\Delta t_{\rm c} = {\rm Min}\left( {(r_i-r_{i+1})^2 \over {D_... ...} \over {(\Omega_i+\Omega_{i+1}) \over 2}}\right\vert}\right ),$$

is equal to 1.24 year, not very different from the one we obtained for the diffusion of the chemical species. It is thus not surprising that the results we obtained are qualitatively similar to those presented in Figs. 5 to 8. The zones where unphysical solution are encountered are the same as those presented in Fig. 4. We show in Table 1 the value of the quantity $\Delta B$ which is equal to

$$\Delta B={B_{\rm final}-B_{\rm initial} \over B_{\rm initial}},$$

where $B_{\rm initial}$, $B_{\rm final}$ are the total angular momentum of the star at the beginning, respectively at the end of the computation. The relative error on the total angular momentum is much greater than the error on the sum of the abundances. This arises because the total angular momentum is an integrated value over the whole system, implying that the errors can accumulate. The sum of the abundances, instead, is evaluated locally in a given shell.

We see that the finite differences method appear to better conserve the total angular momentum. The Crank-Nicholson finite differences method gives the best results, followed by the implicit finite differences method which reaches the same level of accuracy in most cases. Thus the implicit finite differences method that we recommended on the base of the results obtained for the diffusion of the chemical species gives also very satisfactory results for the diffusion of the angular momentum.

6 Conclusion

From the above numerical experiments, it appears that the implicit finite differences method seems to be the most robust one, giving physical solutions in all the cases studied here spanning more than nine orders of magnitude in $\tau$ and $\Delta t$. It has moreover the following characteristics: 1) it reduces the problem to a first order differential equation, 2) it enables an easy and clear interpretation of what happens physically in the system, 3) it is quite easy to implement in a code, 4) it conserves reasonably well the sum of the abundances at each mesh point as well as the total angular momentum.

On the basis of the new tests and analysis made here, we can recommend this method to resolve the diffusion equation in stellar interiors. Of course many more tests could be performed by changing the initial conditions. We restrained our discussion here on a case with a very sharp gradient in order to test the different numerical methods in some extreme conditions. Adopting initially shallower gradients would facilitate the finding of a solution for all three methods for all three methods and they would give identical results.

The diffusion coefficient between two mesh points has to be evaluated correctly to obtain reliable results. The mean diffusion coefficient is equal neither to the algebraic nor to the geometric mean. Its expression is given in Eq. (51).

References

• Battaner, E. 1996, Astrophysical Fluid Dynamics (Cambridge University Press) (In the text)
• Chaboyer, B., & Zahn, J.-P. 1992, A&A, 253, 173 [NASA ADS] (In the text)
• Charbonnel, C. 1995, ApJ, 453, L41 [NASA ADS] (In the text)
• Denissenkov, P. A., Ivanova, N. S., & Weiss, A. 1999, A&A, 341, 181 [NASA ADS] (In the text)
• Endal, A. S., & Sofia, S. 1978, ApJ, 220, 279 [NASA ADS] (In the text)
• Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] (In the text)
• Kawaler, S. D. 1995, in Stellar Remnant, Saas Fee Advanced Course 25, ed. G. Meynet, & D. Schaerer (Springer), 1 (In the text)
• Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [NASA ADS] (In the text)
• Meynet, G., & Maeder, A. 2000, A&A, 361, 101 [NASA ADS] (In the text)
• Meynet, G., & Maeder, A. 2003, A&A, 404, 975 [NASA ADS] (In the text)
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in Fortran, The Art of Scientific Computing (Cambridge: University Press) (In the text)
• Schatzman, E., Maeder, A., Angrand, F., & Glowinski, R. 1981, A&A, 96, 1 [NASA ADS] (In the text)
• Talon, S., & Zahn, J. P. 1997, A&A, 317, 749 [NASA ADS] (In the text)
• Vauclair, S., & Charbonnel, C. 1998, ApJ, 502, 372 [NASA ADS] (In the text)
• Zahn, J.-P. 1992, A&A, 265, 115 [NASA ADS] (In the text)
• Zienkiewicz, O. C., & Taylor, R. L. 2000, The Finite Element Method, Butterworth-Heinemann (Oxford) (In the text)