A&A 366, 166-173 (2001)
DOI: 10.1051/0004-6361:20000219

Nonradial nonadiabatic stellar pulsations: A numerical method and its application to a $\beta$ Cephei model

M. A. Dupret

Institut d'Astrophysique et de Géophysique de l'Université de Liège, avenue de Cointe 5, 4000 Liège, Belgium

Received 25 February 2000 / Accepted 7 November 2000

Abstract
A new general method for the computation of nonradial nonadiabatic oscillations of a given stellar model is presented for a linear approximation. A simple and useful modelling of the atmosphere is included, allowing to obtain credible values for the eigenfunctions in the atmosphere. Some of the results obtained for a $10~M_{\odot}$ model are shown as an illustration. Our study opens the way to different applications. Better theoretical line-profile variations could be obtained from our method, allowing a more detailed comparison with observations. More generally, our study is relevant for asteroseismology, giving a way for a better knowledge of stellar interiors.

Key words: stars: interiors - stars: oscillations

1 Introduction

The problem of finding precise theoretical values for the nonradial nonadiabatic eigenfunctions of a given stellar model has become very interesting and useful for astronomers. Great progress in the study of pulsating stars has been made from an observational point of view. On the one hand, precise periods for multiperiodically pulsating stars have been recently detected from photometric data (e.g. Breger et al. 1999). On the other hand, the observation and analysis of line-profile variations for variable stars is in great development. In those observations, a trace of nonradial oscillations is now clearly admitted, and methods of mode identification have been developed (Aerts 1996; Telting & Schrijvers 1997). So, a confrontation with theory, followed by the improvement of stellar interiors' knowledge (asteroseismology) becomes possible. First steps in this way have already been made (Dziembowsky & Jerzykiewicz 1996; Shibahashi & Aerts 2000).

Concerning the study of line profiles, methods have been developed in order to obtain theoretical line-profile variations comparable with observations. As data for these methods, theoretical values of the nonradial eigenfunctions in the atmosphere are necessary. Until now, only adiabatic eigenfunctions have been used. It is well known, however, that in the atmosphere of a star, the adiabatic approximation gives inaccurate results for some eigenfunctions ( $\delta T/T$, $\delta \rho/\rho$, $\delta \vec{F}$...). The development of a code computing the nonradial eigenfunctions of a stellar model in the linear nonadiabatic case and giving credible values in the atmosphere implies an important step forward. This is the goal of our study.

Dziembowski (1977) was one of the first to study the computation of nonradial nonadiabatic modes. Since then, other authors have worked on the problem (Saio & Cox 1980; Pesnell 1990). In comparison with those previous works, the particularity of our study is, on the one hand, that it puts the emphasis on a new, clear and complete presentation of the problem (Sects. 2 and 4) and of the method of solution (Sect. 5), using a pure Lagrangian formalism; and on the other hand, that it includes a useful treatment of the perturbed equations in the atmosphere (Sect. 3). The results in Sect. 6 are given as an illustration of the quality and the possibilities of our method. We have chosen a $10~M_{\odot}$ model near the turn-off for this application, which is a good candidate for the modelling of a $\beta$ Cephei variable star. We do not propose a modelling of the perturbed convection in this study, so it is not applicable for stars having an important superficial convection zone as the Sun. The presence of a central convection zone (as for $\beta$ Cephei stars) does not pose a problem, for it is located in the quasi-adiabatic region.

   
2 Basic equations

We examine the equations governing the behaviour of a given nonradial mode for which the angular dependence can be expressed by the spherical harmonic $Y_{l m}(\theta,\phi)$and the time dependence by the factor ${\rm e}^{{\rm i} \omega t}$. We use the notation $\delta X$ for the Lagrangian variation of X. First we have the conservation of momentum equations:

   

$$\omega^2 \delta r \frac{\partial\delta\psi}{\partial r} + \frac{\partial(\delta P/P... ... + \left( \frac{\delta\rho}{\rho} - \frac{\delta P}{P} \right) \frac{G m}{r^2},$$ = (1)

$$\omega^2 \delta \theta \frac{1}{r^2} \frac{\partial}{\partial\theta} \left( \delta\psi + \frac{\delta P}{\rho} \right),$$ = (2)

$$\omega^2 \delta \phi \frac{1}{r^2 \sin^{2}\theta} \frac{\partial}{\partial\phi} \left( \delta\psi + \frac{\delta P}{\rho} \right)$$ = (3)

($\psi$, P, $\rho$ and m are respectively the gravitational potential, the pressure, the density and the mass of the sphere of radius r). Next we have the conservation of mass equation:

$$\delta\rho + \rho \left[ \frac{1}{r^2}\frac{\partial}{\partia... ...ight) + \frac{\partial\delta\phi}{\partial\phi} \right] = 0 ,$$ (4)

which gives, using the values of $\delta\theta$ and $\delta\phi$ given by Eqs. (2) and (3):

$$\frac{\delta\rho}{\rho} + \frac{1}{r^2}\frac{\partial}{\part... ...2 r^2} \left( \delta\psi + \frac{\delta P}{\rho} \right) = 0 .$$ (5)

Then we have the perturbed Poisson equation: 2pt

 

$$4 \pi G \ \delta\rho \frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \frac{\partial\delta\psi}{\partial r} \right) - \frac{l(l+1)}{r^2} \delta\psi$$ =

$$\quad +\, \frac{G m}{r^2} \left[ 2 \frac{\partial}{\partial r} \left(\f... ...rac{\partial^2 \delta r }{{\partial r}^2} + \frac{l(l+1)}{r^2} \delta r \right]$$

$$\quad -\, 8 \pi G \rho \frac{\partial\delta r}{\partial r}\cdot$$ (6)

And finally, perturbing the conservation of energy equation $T\,{\rm d} S / {\rm d} t = \epsilon-\nabla\cdot\vec{F}/\rho$with $\vec{F}=-\nu\,\nabla T$, we find: 2pt

 

$${\rm i}\,\omega\,T c_{\rm v}\left(\frac{\delta S}{c_{\rm v}}\right) \epsilon\left(\frac{\delta\epsilon}{\epsilon} +\frac{\delta\rho}{\rho}+2\frac{\partial\delta r}{\partial r} \right)$$ =

$$\quad +\frac{1}{\rho r^2}\frac{\partial}{\partial r} \left[r^2\left(\de... ...c{\partial T}{\partial r} +\nu\frac{\partial\delta T}{\partial r}\right)\right]$$

$$\quad +\frac{L}{4\pi\rho r^3}\left[r^2\frac{{\partial}^2}{\partial r^2} \left(\frac{\delta r}{r}\right)-l(l+1)\frac{\delta r}{r} \right]$$

$$\quad - l(l+1)\frac{\nu}{\rho r^2}\ \delta T$$ (7)

($\epsilon$ is the rate of gain of heat per unit mass due to the thermonuclear reactions, T is the temperature, L the luminosity at radius r and, in the diffusion approximation, $\nu=4 a c T^3 / (3 \kappa \rho)$).

Equations (6) and (7) could seem complicated in comparison with their Eulerian formulation. But we think, our goal being a precise numerical solution of the problem, that a Lagrangian formulation of the equations is better. The reasons are that the surface boundary conditions and the perturbed equations of state appear naturally in a Lagrangian description. Let us consider for example the relation $\delta P=P'+({\rm d}\ln P / {\rm d}\ln r)\ \delta r$. The factor ${\rm d}\ln P / {\rm d}\ln r$ is very large in some parts of a star, going up to values of the order of 10³. So a small value of $\delta P$ could lead to a large value of P' and a first order approximation would no longer be good in an Eulerian description.

In what follows, the star will be subdivided in two parts: the interior from the center to a chosen optical depth $\tau\approx 1$, and the atmosphere from this optical depth to $\tau\approx 0.01$. The frontier between interior and atmosphere will be named the connecting layer and the last layer of the model (at $\tau\approx 0.01$) the surface.

   
3 Treatment of the atmosphere

In the interior of a massive near main sequence star, except for the central convection zone, we can use the diffusion approximation

$$\vec{F}=-\frac{4 a c T^3}{3 \kappa\rho}\nabla T ,$$ (8)

which gives $\nu$ and $\delta\nu$ used in Eq. (7). The problem is that this approximation is not valid in the atmosphere. So we have to adopt another approximation there. A very complex way would be to perturb the radiative transfer equations for each light frequency, and by a judicious integration with respect to the light frequency to find integro-differential equations governing the behaviour of $\delta T$ and $\delta \vec{F}$. The a posteriori determination of the perturbed flux at a given light frequency, knowing the other perturbed variables ( $\delta\vec{r}$, $\delta T$, ...) is already a complex problem (see for example Toutain et al. 1999). So it is easy to see that the general problem of determining simultaneously the different perturbed variables and the perturbed flux is very complicated and goes beyond current studies.

3.1 Our approximation

We propose to treat the atmosphere in the following way. The starting approximation is to suppose that the atmosphere is always in thermal equilibrium during the pulsation. This approximation can be justified by comparing the fundamental pulsation period of a star, $\tau_{\rm p}$, to the thermal relaxation time of its atmosphere,

$$\tau_{\rm th a}=\int_{m_{\rm a}}^{M}\, T\,c_{\rm v}\,{\rm d}m\ /\ L\ ,$$ (9)

where $m_{\rm a}$ is the mass between the center and the basis of the atmosphere and M is the total mass of the star. We find that for a typical massive star near the main sequence,

$$\tau_{\rm th a}\simeq 10^{-4}\ \tau_{\rm p}\ .$$

In an equilibrium model we can write

$$T^4(\tau)=T_{\rm eff}^4\ f(\tau)$$ (10)

(this is not an approximation but a definition of $f(\tau)$).
Now, the approximation we take is to suppose that Eq. (10) is also valid in the perturbed model, with the same $f(\tau)$as for the equilibrium model. As we have said this corresponds in a way to the assumption that the atmosphere is always in thermal equilibrium during the pulsation. In this approximation, we can perturb Eq. (10) and find:

$$4\frac{\delta T}{T}=4\frac{\delta T_{\rm eff}}{T_{\rm eff}} +\frac{{\rm d}\ln f}{{\rm d}\tau}\ \delta\tau .$$ (11)

In Eq. (11), $\delta T_{\rm eff} / T_{\rm eff}$ is proportional to $Y_{l m}(\theta,\phi)$. Making a sufficiently precise atmosphere model at equilibrium, we find sufficiently precise values for the $T(\tau)$ law and its first and second derivatives, as they will be required in the next equations.

Next, we examine the definition of $\tau$ in the perturbed model. More precisely, we have ${\rm d}\tau = - \kappa\rho\ {\rm d}s$, but what is the direction of ${\rm d}s$ we take? We could choose that ${\rm d}s$ has the same direction as the flux vector or that it is perpendicular to the constant T surface or that it is purely radial, etc. But in all these definitions, the angle between ${\rm d}s$ and ${\rm d}r$ is of the first order and so (using the "$\, '\,$'' for Eulerian variations)

$$\frac{{\rm d}\tau '}{{\rm d}s}-\frac{\partial\tau '}{\partial... ...}\delta\tau}{{\rm d}s} -\frac{\partial\delta\tau}{\partial r}$$

are of the second order and can be neglected in our first order approximation. We obtain then after some algebra,

$$\frac{\partial\delta\tau}{\partial r}=-\kappa\rho \left(\fr... ...lta\rho}{\rho} +\frac{\partial\delta r}{\partial r} \right) .$$ (12)

Next, we isolate $\delta\tau$ in Eq. (11), take the derivative with respect to r, and eliminate $\partial\delta\tau / \partial r$ using Eq. (12). We find then 2pt

 

$$\frac{\partial (\delta T / T)}{\partial r} \frac{{\rm d}\ln T}{{\rm d}r} \left[\left(\frac{\delta T_{\rm eff... ...{{\rm d}}{{\rm d}\tau}\left(\frac{1}{{\rm d}\ln T / {\rm d}\tau}\right) \right.$$ =

$$\left. \ \qquad\qquad\qquad\qquad +\frac{\delta\kappa}{\kappa} +\frac{\delta\rho}{\rho}+\frac{\partial\delta r}{\partial r} \right].$$ (13)

This is the equation we will use instead of Eq. (7) in the atmosphere.

3.2 Surface boundary conditions

Let us consider first the thermal surface boundary condition. An accurate and rigorous formulation of it in the general case of nonradial oscillations has been proposed by Gabriel (1989). The principle he adopted is to impose that in the very outer layers of a star, the progressive waves of the radiation field are outgoing ones, with no incoming component. A problem is, that the very outer layer, where the condition has to be imposed, is far away from the photosphere. Thus, in order to use Gabriel's formulation, we have to know precisely the perturbed radiation field between the photosphere and this outer layer. As we have said, to derive equations independent of light frequency governing the field is very difficult and goes beyond the approach we have adopted here.

In our approximation, the flux can be obtained by the diffusion approximation in the interior of the star, but not in the atmosphere. So the only place where a boundary condition can be imposed on the flux is at the connecting layer between atmosphere and interior. From the diffusion approximation, we have on the one hand:

$$\vert\vec{F}\vert^2=\nu^2 \vert\nabla T\vert^2 ,$$ (14)

and from Stephan's law on the other hand:

$$\vert\vec{F}\vert=\sigma T_{\rm eff}^4 .$$ (15)

Perturbing these two equations, we find: 2pt

 

$$\frac{\delta \vert\vec{F}\vert}{\vert\vec{F}\vert} \frac{\delta\nu}{\nu} +\frac{\partial\delta T/\partial r}{{\rm d}T/{\rm d}r} -\frac{\partial\delta r}{\partial r}$$ =

$$4\frac{\delta T_{\rm eff}}{T_{\rm eff}}\cdot$$ = (16)

At the surface, we have also the very natural condition $\delta\tau =0$, which gives, using Eq. (11):

$$\frac{\delta T}{T}=\frac{\delta T_{\rm eff}}{T_{\rm eff}}\cdot$$ (17)

Equations (16) and (17) are the thermal boundary conditions we adopt. The first is imposed at the connecting layer, the second at the surface.

Let us consider now the dynamic boundary conditions. We define R and M as the values of r and m at the surface. For the Eqs. (1) and (5), we proceed as suggested by Cox (1980, Sect. 17.6b). We obtain the boundary conditions assuming that

$$\frac{P}{\rho}\frac{R}{G M} \ \to\ 0$$

and that $\delta P / P$ has a finite value at the surface, which gives for the conservation of radial momentum:

$$\omega^2 \delta r = \frac{\partial\delta\psi}{\partial r} +... ...elta\rho}{\rho} - \frac{\delta P}{P} \right) \frac{G M}{R^2}$$ (18)

and for the conservation of mass:

$$\frac{\delta\rho}{\rho} + \frac{1}{r^2}\frac{\partial}{\part... ...ta r \right) - \frac{l (l+1)}{\omega^2 r^2} \ \delta\psi = 0 .$$ (19)

Finally, for the Poisson equation, we use the classical condition (in Eulerian description):

$$\frac{\partial\psi '}{\partial r}+\frac{l+1}{R}\ \psi ' =-4\pi G\rho\ \delta r .$$ (20)

It is obtained by imposing a first order continuous match (continuity of $\psi$ and its first derivatives) between the interior solution of the Poisson equation and the decreasing solution of the Laplace equation (Ledoux & Walraven 1958, Sect. 75). Equation (20) in a Lagrangian description becomes:

$$\frac{\partial\delta\psi}{\partial r}+\frac{l+1}{R}\ \delta\... ...a r + \frac{G M}{R^2} \frac{\partial\delta r}{\partial r}\cdot$$ (21)

   
4 Dimensionless formulation and central boundary conditions

The central boundary conditions are obtained by imposing the solutions to be regular at the center (finite values for the perturbed variables and their derivatives). After some algebra, it can be shown that the different variables can be rewritten in the following form:

$$\delta r R\ x^{l-1}\ \zeta , \quad \frac{\delta\rho}{\rho}=x^l\ \gamma , \quad \frac{\delta T}{T}=x^l\ \vartheta ,$$ =

$$\frac{\delta S}{c_{\rm v}} x^l\ \eta , \quad \delta\psi=\frac{G M}{R}\ x^l\ \varphi \quad{\rm and}\quad \frac{\delta P}{P}=x^l\ \beta$$ = (22)

where $\zeta$, $\gamma$, $\vartheta$, $\eta$, $\varphi$ and $\beta$are finite at the centre. A justification of those results can be found, for example, in Unno et al. (1989). Defining the following dimensionless symbols: 2pt

$${\lambda}_{\rm p} \frac{R}{G M}\ \frac{P}{\rho} , \quad \Omega^2=\frac{R^3}{G M}\ \omega^2 , \quad q=\frac{m}{M} ,$$ =

$$\frac{r}{R} \quad{\rm and}\quad {\rho}_{\rm m} = \frac{4\pi R^3}{M}\rho ,$$ x =

the Eqs. (1), (5), (6), (7) and (13) can be rewritten in a dimensionless form. We keep only the r dependent part of the variables, so the partial derivatives are transformed into simple derivatives. In what follows, a prime will denote the derivative with respect to x(for example: $\zeta'\equiv{\rm d}\zeta/{\rm d}x$). We assume that all the variables have a derivative equal to zero at the centre, which gives (here for $\zeta$): $\lim_{x\to 0}\zeta'/x=\zeta''$.
Equation (1) gives then:

$$\Omega^2 \zeta = l(\varphi+{\lambda}_{\rm p} \beta ) +x(\varphi'+{\lambda}_{\rm p}\beta') +(\gamma - \beta)q/x ,$$ (23)

which gives the central boundary condition:

$$\Omega^2 \zeta = l(\varphi+{\lambda}_{\rm p} \beta ).$$ (24)

Equation (5) gives:

$$\Omega^2\left(x^2\gamma+x\zeta'+(l+1)\zeta\right) =l(l+1)(\varphi+{\lambda}_{\rm p}\beta).$$ (25)

Substituting the value of $\Omega^2\zeta-l(\varphi+{\lambda}_{\rm p}\beta)$given by Eq. (23) in Eq. (25) and dividing by x², we find the central boundary condition:

$$\Omega^2\left(\gamma+\zeta''\right) +(l\!+\!1)\left(\varphi'... ..._{\rm p}\beta'' +(\gamma\!-\!\beta){\rho}_{\rm m}/3\right)=0.$$ (26)

Equation (6) gives: 1pt

 

$$\varphi''+2(l+1)\frac{\varphi'}{x}\ =\ {\rho}_{\rm m}\gamma \frac{q}{x^3}\zeta''$$ +

$${}+2\left({\rho}_{\rm m}+(l-2)\frac{q}{x^3}\right) \frac{\zeta'}{x} 2(l-1)\frac{{\rho}_{\rm m}-3 q/x^3}{x^2}\zeta .$$ + (27)

It can be shown that

$$\lim_{x\to 0}\frac{{\rho}_{\rm m}-3 q/x^3}{x^2} =\frac{{\rho}_{\rm m}''(0)}{5}\cdot$$ (28)

So we find at the centre:

$$(2 l+3)\varphi''={\rho}_{\rm m}\gamma +\frac{2 l+3}{3}{\rho}_{\rm m}\zeta'' +\frac{2(l-1)}{5}{\rho}_{\rm m}''\zeta .$$ (29)

For the conservation of energy (Eq. (7)), we introduce the following dimensionless symbols and variables: 2pt

$$\sqrt{\frac{R^3}{G M}}\,\frac{L(r)}{M T c_{\rm v} {\rho}_{\rm m} x^3} , \ \epsilon_0=\frac{M x^3}{L(r)}\epsilon ,$$ k =

$${\epsilon}_{1} \frac{\epsilon_0{\rho}_{\rm m}-3}{x^2} , \ T_1 = \frac{x}{{\rm d}\ln T/{\rm d}x}, \ T_2=\frac{1}{x}\frac{{\rm d}\,T_1}{{\rm d}x},$$ =

$$\delta\epsilon_x \frac{1}{x^l}\frac{\delta\epsilon}{\epsilon} \quad {\rm and}\quad \delta\nu_x = \frac{1}{x^l}\frac{\delta\nu}{\nu}\cdot$$ =

Equation (7) gives then: 2pt

 

$$\frac{i\,\Omega\,\eta}{k} \epsilon_0\,{\rho}_{\rm m}\, (\delta\epsilon_x+\gamma) +2(l-1){\epsilon}_{1}\,\zeta$$ =

$$\quad +\, 2(\epsilon_0\,{\rho}_{\rm m}+l-2)\,\zeta'/x +\zeta''$$

$$\quad -\, l\left({\epsilon}_{1}T_1+T_2\right)\vartheta -(\epsilon_0\,{\rho}_{\rm m}+l) (\delta\nu_x+\vartheta)$$

$$\quad -\, (\epsilon_0{\rho}_{\rm m}+2 l-1) T_1\, \vartheta'/x - x(\delta\nu_x'+\vartheta')$$

$$\quad -\, x\, T_2\,\vartheta' - T_1\,\vartheta'' .$$ (30)

Noting that k, ${\epsilon}_{1}$, T₁ and T₂ have finite values when $x\to 0$, we find at the center: 2pt

 

$$\frac{i\,\Omega\,\eta}{k} 3(\delta\epsilon_x+\gamma) +2(l-1){\epsilon}_{1}\,\zeta + (2 l+3)\,\zeta''$$ =

$$\quad -\, l\left({\epsilon}_{1}T_1+T_2\right)\vartheta -(l+3)(\delta\nu_x+\vartheta)$$

$$\quad -\, (2 l+3)T_1\,\vartheta'' .$$ (31)

The appearance of first and second derivatives with respect to x in the coefficient T₁ and T₂ needs to be very careful in order to avoid numerical noise in their estimation. A good way to proceed is to use the analytical values given either by the diffusion approximation in radiative zones or by the adiabatic gradient in convective zones. More precisely, we have in a radiative zone:

$$T_1=-\frac{16 \pi x^3 a c\, T^4 R}{3 \kappa\rho L(r)},$$ (32)

and in a convective zone:

$$T_1=-\frac{\Gamma_1}{\Gamma_3-1}\frac{x^3}{q}\,\lambda_{\rm p}.$$ (33)

For our approximation in the atmosphere, defining $\delta\kappa_x=(\delta\kappa/\kappa)/x^l$ and using Eq. (17), Eq. (13) gives: 2pt

 

$$\vartheta' \frac{l}{x}\,\vartheta - \frac{{\rm d}\ln T}{{\rm d}x} \left(\lef... ...{{\rm d}}{{\rm d}\tau}\left(\frac{1}{{\rm d}\ln T / {\rm d}\tau}\right) \right.$$ +

$$\left. \qquad\qquad\qquad {}+ \delta\kappa_x +\gamma+\frac{\zeta'}{x}+\frac{l-1}{x^2}\zeta \right)=0 .$$ (34)

We use the notation $\vartheta_{\rm s}$ for the value of $\vartheta$at the surface and $\vartheta_{\rm p}$ for the value of $\vartheta$ at the connecting layer. Finally, Eq. (16) gives:

 

$$4\vartheta_{\rm s}/x_{\rm p}^l {(\delta\nu_x)}_{\rm p} +\left(1+\frac{l}{x_{\rm p}\, ({\rm d}\ln T/{\rm d}x)_{\rm p}}\right) \,\vartheta_{\rm p}$$ =

$$\quad +\, \frac{\vartheta'_{\rm p}} {({\rm d}\ln T/{\rm d}x)_{\rm p}} - \frac{\zeta'_{\rm p}}{x_{\rm p}} - \frac{l-1}{x_{\rm p}^2}\zeta_{\rm p} .$$ (35)

To all these differential equations we have to add some algebraic equations. They consist of the perturbed state equations:

   

$$\beta P_{\rm T}\,\eta+\Gamma_1\,\gamma , \quad \vartheta=\eta+(\Gamma_3-1)\,\gamma,$$ = (36)

$$\delta\kappa_x \kappa_{\rm T}\, \eta+\kappa_{\rm\rho S}\ \gamma , \quad \delta\nu_x=3\,\vartheta -\delta\kappa_x-\gamma\quad{\rm and}$$ = (37)

$$\delta\epsilon_x \epsilon_{\rm T}\,\eta +\epsilon_{\rm\rho S}\ \gamma .$$ = (38)

The way to obtain $\epsilon_{\rm T}$ and $\epsilon_{\rm\rho S}$ is explained in Ledoux & Walraven (1958, Sect. 66). As we can see, $\delta\nu_x$ was obtained from the perturbation of the diffusion approximation (Eq. (8)). When a convection zone is present, we use the same equation. Such a simplification does not pose a problem in central convection zones, for they are located in a quasi-adiabatic region where the conservation of energy equation plays a very small role.

   
5 Method of solution

Equations (23), (25), (27) and (30) (or (34) in the atmosphere) form a system of 4 linear ordinary differential equations of first and second order with 4 boundary conditions at the centre (Eqs. (24), (26), (29) and (31)), 3 boundary conditions at the surface (dimensionless form of Eqs. (18), (19) and (21)) and 1 condition at the connecting layer (Eq. (35)). Using Eqs. (36), (37) and (38), we can express these 4 differential equations in terms of the 4 unknown functions $\zeta$, $\varphi$, $\gamma$ and $\eta$. Because of the presence of the factor $\Omega$, this system of equations can be seen as an eigenvalue problem. It can be shown that the number of boundary conditions is in agreement with the number of differential equations and their order (some of the central conditions, Eqs. (26), (29) and (31), are not really boundary conditions in the usual sense, their order being the same or larger than the one of the differential equations).

Let us examine now how we solve numerically this problem. We use a finite differences method in order to express our problem as a finite dimensional eigenvalue problem. So first we have to choose a good grid $\ x_1=0,\,\ldots$, $x_i=r_i/R,\,\ldots$, $x_N=1\$ where the discrete variables will be defined. Then we have to approximate at each point of the grid the derivatives appearing in the differential equations by finite differences between these variables, and our problem will be rewritten as a difference scheme. Finally, introducing 2 additional variables, it will take the canonical form

$$\tens{A} \vec{x}\ =\ \Omega\,\tens{B}\,\vec{x}$$ (39)

where $\tens{A}$ and $\tens{B}$ are complex square matrices. Equation (39) is a generalization of the classical eigenvalue problem $\tens{A} \vec{x}=\lambda \vec{x}$. In order to calculate the eigenvectors ($\vec{x}$) and the eigenvalues ($\Omega$) which interest us, we will use then a generalization of the classical inverse iteration algorithm adapted to the more general problem Eq. (39).

5.1 The discretization

Let us examine first the choice of the number of layers (N) and of the distribution of the different points $x_1,\,\ldots\!,\,x_N$. In order to have sufficiently precise results, the number of layers has to be larger than the one usually used in evolution codes. A number $N\simeq 2000$ turns out to be appropriate and leads to calculation times quite reasonable ($\simeq$10 s for one mode).

Let us examine now the discretization of the differential equations. The most important quality criterion of a difference scheme is its stability, which is even more important than its precision. We insist on this point because in our case the difference scheme obtained is easily subject to instabilities, which can lead to completely wrong results. This is the case, for example, if we use the estimation

$$\frac{{\rm d}f}{{\rm d}x}(x_i) \simeq \frac{f(x_{i+1})-f(x_{i-1})}{x_{i+1}-x_{i-1}}$$ (40)

or its generalization for the case of non-equidistant points.
On the other hand, the estimations

$$\frac{{\rm d}f}{{\rm d}x}(x_i) \simeq \frac{f(x_{i+1})-f(x_{i})}{x_{i+1}-x_{i}}$$ (41)

or

$$\frac{{\rm d}f}{{\rm d}x}(x_i) \simeq \frac{f(x_{i})-f(x_{i-1})}{x_{i}-x_{i-1}} ,$$ (42)

are proved to be more stable and thus much better despite their lower precision order. Note that for the case of radial oscillations, it is possible to combine precision and stability using an interlaced difference scheme (Castor 1971) but it is not so easy for the nonradial case. Using Eq. (41) for all the equations gives a bad linking with the central boundary conditions, and using Eq. (42) for all the equations gives a bad linking with the surface boundary conditions. The intermediate solution we have adopted is to use Eq. (41) for the discretization of Energy and Poisson equations (Eqs. (30), (34) and (27)) because they have to be well connected to the surface boundary conditions Eqs. (21) and (35), and to use Eq. (42) for the others. For the boundary condition on the flux (Eq. (16)), we use Eq. (42), in order to have a good linking between the interior and the atmosphere.

Finally, our problem can be rewritten as a difference scheme consisting of a system of 4 N+2 linear equations with 4 N+2 unknowns. More precisely, let $N_{\rm c}$ be the index of the connecting layer, we have the following variables: $\zeta_i\ {\scriptstyle (1\le i\le N)}$, $\,\varphi_i\ {\scriptstyle (1\le i\le N)}$, $\,\gamma_i\ {\scriptstyle (1\le i\le N)}$, $\, \eta_i\ {\scriptstyle (1\le i\le N_{\rm c}+1)}\$ and $\vartheta_i\ {\scriptstyle (N_{\rm c}\le i\le N)},$ and concerning the equations, we have: N equations for the discretization of the Poisson equation and its boundary conditions (Eqs. (27), (29) and (21)), N equations for the discretization of the radial component of the equation of motion and its boundary conditions (Eqs. (23), (24) and (18)), N equations for the discretization of the continuity equation and its boundary conditions (Eqs. (25), (26) and (19)), $N_{\rm c}$ equations for the discretization of the conservation of Energy equation in the interior and its central boundary condition (Eqs. (30) and (31)), $N-N_{\rm c}-1$ equations for the discretization of the thermal equation we use in the atmosphere (Eq. (34)), 1 equation for the discretization of the flux boundary condition (Eq. (35)) and 2 equations for the linking between $\eta$ and $\vartheta$at the connecting point between interior and atmosphere. An inaccurate way to proceed is to eliminate the N equations corresponding to the Poisson equation and the N variables $\gamma_i$, substituting in the other equations the value of $\gamma_i$ given by the discretization of Eq. (27). The principal reason is that this substitution corresponds to a very bad pivot choice in the superficial layers of a star ( $\rho_{\rm m}\simeq 10^{-8}$ in these layers).

As we have said, we want to express our problem in the canonical form $\tens{A} \vec{x}=\Omega\,\tens{B}\,\vec{x}$. In order to do this, we have to introduce two additional variables at each layer. More precisely, an $\Omega^2$ dependence appears in the discretization of the radial component of the equation of motion and of the continuity equation. For the former equation, we eliminate it introducing the new variables $\xi_i=\Omega\zeta_i$, and for the latter equation, we eliminate it by splitting at each layer i the equation

$$\vec{a}_i^{\rm T} \vec{x}=\Omega^2\,\vec{b}_i^{\rm T} \vec{x}$$ (43)

in the 2 equations

$$\vec{a}_i^{\rm T} \vec{x}=\Omega\,\mu_i \quad {\rm and} \quad \mu_i=\Omega\,\vec{b}_i^{\rm T} \vec{x}.$$ (44)

In this way, our problem finally takes the required form $\tens{A} \vec{x}=\Omega\,\tens{B}\,\vec{x}$, where $\tens{A}$ and $\tens{B}$ are square complex matrices of order 6 N+2 with a bandwidth of 18 (some lines of $\tens{B}$ are equal to zero).

5.2 The numerical algorithm

In order to solve Eq. (39) we use a generalization of the inverse iteration algorithm. Departing from the estimation $\vec{x}_k$ of the eigenvector at the step k and the initial estimation $\Omega_0$ of the eigenvalue, the next estimation is obtained by the formula

$$\vec{x}_{k+1}=(\tens{A} -\Omega_0 \tens{B})^{-1}\tens{B}\,\vec{x}_k.$$ (45)

Under the hypothesis that $\tens{A}^{-1}\tens{B}$ is diagonalizable, it is easy to prove that this algorithm has the same rate of convergence as the classical inverse iteration algorithm. We do not calculate explicitly the inverse matrix but resolve the linear system

$$(\tens{A}-\Omega_0 \tens{B})\,\vec{x}_{k+1}=\tens{B}\,\vec{x}_k.$$ (46)

Because $\tens{A}-\Omega_0 \tens{B}$ is very badly conditioned, we begin by doing a suitable preconditioning of it. Then we do a $\tens{L} \tens{U}$ factorization of it with partial pivoting. And finally we iterate, solving the 2 triangular systems at each step (see for example Wilkinson (1965) for an explanation of this numerical algorithm). The eigenvector and eigenvalue initial estimations can be obtained by solving first the problem in the adiabatic approximation. Note that if the eigenvalue initial estimation is good, the algorithm converges quickly even with a very bad eigenvector initial estimation. As we have said, the calculation time for the computation of one mode is very short ( $t\approx 10$ s). So it is easy to examine the results obtained for a lot of models.

	Figure 1: Real parts of the nonadiabatic eigenfunctions $\delta S/c_{\rm v}\,({\rm d}S)$, $\delta T/T\,({\rm d}T)$ and $\delta T/T\vert _{\rm ad}$ ( ${\rm d} T_0)$ corresponding to the mode p1 with l=1
Open with DEXTER

   
6 Results

We now describe some of the results obtained for a $10~M_{\odot}$ model with Z=0.02. We have chosen a model at the end of the central hydrogen burning, for it corresponds to a region of the HR diagram where $\beta$ Cephei variable stars are found. More precisely, it represents a star with $L=4.35~10^{37}~{\rm erg~s^{-1}}$, $T_{\rm eff}=22\,640~{\rm K}$ and for which the age after the beginning of the hydrogen burning is $1.845\,10^7$ years. As it is well known, many $\beta$ Cephei stars have simultaneous radial and nonradial pulsations.

In Figs. 1 and 2, we give the graphs of the real parts of the nonadiabatic eigenfunctions $\delta S/c_{\rm v}\,({\rm d}S)$ and $\delta T/T\,({\rm d}T)$. We give also for a comparison the graphs of $\delta T/T$ obtained in the adiabatic approximation (refered to by ${\rm d} T_0$). Figure 1 corresponds to the nonradial mode p₁ with l=1and Fig. 2 to the fundamental radial mode. The labels at the bottom correspond to the logarithm of the equilibrium temperature. The labels at the top correspond to the logarithm of the thermal relaxation time of the upper layers ( $\tau_{\rm th}=\int_{m}^{M}\,T\,c_{\rm v}\,{\rm d}m/L$) divided by the dynamical time of the star ( $\tau_{\rm dyn}=\sqrt{R^3/(G M)}$). Those functions are normalized in such a way that $\delta r/r=1$ at the surface. In Fig. 2, we have also given $\vert\delta L/L\vert$(divided by 5). The graphs of $\delta S/c_{\rm v}$ show that from the center to $\log T\simeq 5.5$, the adiabatic approximation is valid. The comparison between the adiabatic and nonadiabatic values of $\delta T/T$ shows that the adiabatic approximation gives wrong results from $\log T\simeq 5.5$to the surface. From a qualitative point of view, the shapes of the eigenfunctions $\delta S/c_{\rm v}$ and $\delta T/T$ (positions of the zeros and sign of the derivatives) are approximately the same for the different modes we have examined (compare for example Figs. 1 and 2).

	Figure 2: Real parts of the nonadiabatic eigenfunctions $\delta S/c_{\rm v}\,({\rm d}S)$, $\delta T/T\,({\rm d}T)$ and $\delta T/T\vert _{\rm ad}$ ( ${\rm d} T_0)$, and $0.2\, \vert \delta L / L \vert\,(0.2 {\rm d}L)$ corresponding to the fundamental radial mode
Open with DEXTER

If the discretization is done carefully, the derivatives of the eigenfunctions are continuous at the connecting layer between interior and atmosphere and the solutions are not dependent on the choice of this layer. In our application, we took it at $\tau=3$ but the results are virtually the same (indistinguishable on a graph) if we take it for example at $\tau=5$ or $\tau=1$. In Fig. 3, we give the graphs of $\vert\delta T/T\vert$ in the atmosphere for the first 3 nonradial p modes of order l=1 (p₁, p₂ and p₃). The abscissa corresponds to the logarithm of the equilibrium optical depth (Rossland mean). In Figs. 1, 2 and 3, the vertical line corresponds to the connecting layer between interior and atmosphere.

	Figure 3: Values of $\vert\delta T/T\vert$ in the atmosphere for the fundamental radial mode (F) and for the 3 first nonradial p modes of order l=1 (p1, p2 and p3)
Open with DEXTER

 

 

Table 1: Comparison between nonadiabatic and adiabatic amplitudes at the surfaces for different modes

mode	$\left\vert\frac{\delta S}{c_{\rm v}}\right\vert$	$\left\vert\frac{\delta T_{\rm eff}}{T_{\rm eff}}\right\vert$	$\left\vert\frac{\delta T}{T}\right\vert _0$	$\phi_T$	f  ($\mu$Hz)
$l \!=0\! \quad F$	2.62	2.39	3.78	2.818	55.89
$l\!=0\! \quad H_1$	4.66	2.92	5.40	3.092	74.27
$l\!=0\! \quad H_2$	7.75	3.52	6.85	-2.909	91.66
$l\!=1\! \quad p_1$	4.56	2.79	5.97	3.131	78.02
$l\!=1\! \quad p_2$	6.39	3.09	7.09	-3.015	90.35
$l\!=1\! \quad p_3$	8.62	3.30	8.21	-2.876	102.70
$l\!=2\! \quad f$	6.35	3.10	6.97	-3.016	90.36
$l\!=2\! \quad p_1$	7.50	3.23	7.69	-2.942	96.93
$l\!=2\! \quad p_2$	9.52	3.35	9.12	-2.825	107.27
$l\!=2\! \quad p_3$	14.86	3.29	7.90	-2.558	130.22
$l\!=3\! \quad f$	7.26	3.23	7.52	-2.956	95.88
$l\!=3\! \quad p_1$	10.87	3.40	9.24	-2.754	113.69
$l\!=3\! \quad p_2$	13.25	3.37	9.04	-2.636	123.84

In Table 1, the results obtained for different modes of the same model are examined. In the second column, we give the value of $\vert\delta S/c_{\rm v}\vert$ at the photosphere ($\tau=2/3$). In the third column, we give the value of $\vert\delta T_{\rm eff}/T_{\rm eff}\vert$. For the comparison, we give in the fourth column the value of $\vert\delta T/T\vert$ obtained in the adiabatic approximation at the photosphere. In the fifth column, we give the phase of $\delta T_{\rm eff} / T_{\rm eff}$in Radians ( $\phi_{T_{\rm eff}}$). In the last column, we give the frequency of the mode (real part) in $\mu$Hz. We see that, systematically, the value of $\vert\delta T_{\rm eff}/T_{\rm eff}\vert$is smaller than the value of $\vert\delta T/T\vert$ obtained in the adiabatic approximation. These differences are significant and could have a great effect on the obtaining of accurate theoretical line-profile variations. We see also that the phase of $\delta T/T$ is not exactly equal to $\pi$ as it was the case in the adiabatic approximation.

Another interesting function is the work integral. Demonstrating the Hermiticity of the linear adiabatic wave equation (see Cox 1980 or Unno et al. 1989), it is shown that the imaginary part of the eigenvalue $\Omega=\Omega_{\rm r}+{i}\,\Omega_{\rm i}$ can be obtained by an integral expression similar to the well known one of the radial case. Using our dimensionless formulation, we have

$$W(x)=\frac{-1}{2\,\Omega_{\rm r}} \frac{\int_{0}^{x}\,\lambda... ...^4}\vert\varphi+\lambda_{\rm p} \beta\vert^2\right){\rm d} x},$$ (47)

with $W(1)=-\Omega_{\rm i}$.
Multiplied by an appropriate constant, the function W(x) corresponds to the work done by the sphere of radius $r=x\,R$ during one pulsation. We give in Fig. 4 the graph of W obtained for the mode p₁ with l=1 and compare it with the one of the fundamental radial mode. The similarity between the two functions shows that the driving mechanisms are essentially the same for radial and nonradial oscillations. We see that the destabilizing region (where W increases) is located between $\log T=5$ and $\log T=5.5$. A first bump of the opacity explained by the tremendous number of transition iron lines is present in this region (see the graph of $\kappa$ in Fig. 4). In this example, the radial mode is unstable and the nonradial one is stable. All this is in agreement with the now admitted fact that the unstability of $\beta$ Cephei stars is due to the $\kappa$-mechanism located there. Such an explanation was first proposed by Cox et al. (1992).

	Figure 4: Work integrals of the fundamental radial mode ($W_{\rm r}$) and the nonradial mode p1 with l=1 ( $W_{\rm nr}$) compared with the opacity ($\kappa$)
Open with DEXTER

The results shown in this section permit to see clearly that a pulsating star can be subdivided in three different parts. From the center to the surface, we have first the adiabatic region where the thermal capacity is so large that luminosity imbalances have not the time to change significantly the entropy of the matter (see the thermal relaxation times and the graphs of $\delta S/c_{\rm v}$ given in Figs. 1 and 2). Then we have the transition region (between $\log T=5.5$ and $\log T=5$ in this example), in which the thermal relaxation time is of the same order as the dynamical time. In this region, the pulsation becomes clearly nonadiabatic and the heat capacity remains sufficiently important, so that driving and damping mechanisms can influence the all pulsation, making the star stable or unstable (see Fig. 4). And finally we have a region of small heat capacity, where very small luminosity imbalances are sufficient to change significantly the entropy of the matter. In this region, $\delta L/L$ becomes quasi-constant (as function of r) and the hypothesis of thermal equilibrium becomes thus perfectly acceptable (see the thermal relaxation times and the graph of $\delta L/L$in Fig. 2).

7 Conclusions

In conclusion, the method we have developed opens the way to interesting bridges between theory and observations. Our new treatment of the outer atmosphere to calculate the eigenfunctions of a nonradially pulsating stars is a fruitful starting point to improve the current methods of analysis of pulsating stars. Up to now, e.g., the nonadiabatic character of the pulsation was treated with an ad hoc parameter in all mode-identification methods based on multicolour photometry (see e.g. Watson 1988; Heynderickx et al. 1994; Garrido 2000).

Also, the available codes for the calculation of theoretical line-profile variations (e.g. Townsend 1997) rely on adiabatic eigenfunctions. It has been claimed in the literature (Balona 1987; Lee et al. 1992) that nonadiabatic effects in the calculation of a line profile are much more important than the pulsational velocity effects in rapid rotators. However, these results are based on an (each time different) arbitrary parametrization of the eigenfunction of the temperature. Moreover, the discrepancies between observed and theoretical line-profile variations is sometimes ascribed to the neglect of temperature effects (e.g. Aerts et al. 1992). Our method allows us to justify or contradict all these suggestions on a much firmer base than was possible before. To find an answer to the question whether or not it is necessary to use nonadiabatic eigenfunctions in the calculation of line-profile variations would be an important step forward in the interpretation of such kind of data and thus also in the application of asteroseismology. We will elaborate on this problem in a forthcoming paper.

Acknowledgements

This work was supported by the F.R.I.A. (Fonds de Recherche pour l'Industrie et l'Agriculture). We express our thanks to R. Scuflaire who helped us greatly in this study and to A. Noels who proposed us to work on this subject. We would also like to thank C. Aerts who encouraged us to write this paper and J. De Ridder who gave us with C. Aerts useful comments about it. And finally, we thank P. Bradley who gave us also useful comments.

References

 

• Aerts, C., De Pauw, M., & Waelkens, C. 1992, A&A, 266, 294 In the text NASA ADS  
• Aerts, C. 1996, A&A, 314, 115 In the text NASA ADS  
• Ando, H., & Osaki, Y. 1975, PASJ, 27, 581 NASA ADS  
• Balona, L. A. 1987, MNRAS, 224, 41 In the text NASA ADS  
• Breger, M., Handler, G., Garrido, R., et al. 1999, A&A, 349, 225 In the text NASA ADS  
• Castor, J. I. 1971, ApJ, 166, 109 In the text NASA ADS  
• Cox, A. N., Morgan, S. M., Rogers, F. J., & Iglesias, C. A. 1992, ApJ, 393, 272 In the text NASA ADS  
• Cox, J. P. 1980, Theory of Stellar Pulsation (Princeton Univ. Press, Princeton) In the text  
• Dziembowski, W. 1977, Acta Astron., 27, 95 In the text NASA ADS  
• Dziembowski, W., & Jerzykiewicz, M. 1996, A&A, 306, 436 In the text NASA ADS  
• Gabriel, M. 1989, A&A, 208, 122 In the text NASA ADS  
• Garrido, R. 2000, in The 6th Vienna Workshop on $\delta\,$Scuti and related stars, ed. M. Montgommery, & M. Breger, PASP Conf. Ser., in press  
• Heynderickx, D., Waelkens, C., & Smeyers, P. 1994, A&AS, 105, 447 NASA ADS  
• Keeley, D. A. 1977, ApJ, 211, 926 NASA ADS  
• Ledoux, P., & Walraven, Th. 1958, Handbuch der Physik, ed. S. Flügge, 51, 353 In the text  
• Lee, U., Jeffery, C. S., & Saio, H. 1992, MNRAS, 254, 185 In the text NASA ADS  
• Pesnell, W. D. 1990, ApJ, 363, 227 In the text NASA ADS  
• Saio, H., & Cox, J. P. 1980, ApJ, 236, 549 In the text NASA ADS  
• Shibahashi, H., & Aerts, C. 2000, ApJL, in press In the text  
• Telting, J. H., & Schrijvers, C. 1997, A&A, 317, 723 In the text NASA ADS  
• Toutain, T., Berthomieu, G., & Provost, J. 1999, A&A, 344, 188 In the text NASA ADS  
• Townsend, R. 1997, MNRAS, 284, 839 In the text NASA ADS  
• Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars, 2nd edition (Univ. Tokyo Press, Tokyo) In the text  
• Watson, R. D. 1988, Ap&SS, 140, 255 NASA ADS  
• Wilkinson, J. H. 1965, The Algebraic Eigenvalue Problem (Oxford Univ. Press, Oxford) In the text