A&A 484, 419-428 (2008)
DOI: 10.1051/0004-6361:20078496

Linear and nonlinear pulsation models with the variable Eddington factor approximation of radiative hydrodynamics

T. Aikawa

Tohoku Gakuin University, Izumi-ku, Sendai, 981-3193, Japan

Received 17 August 2007 / Accepted 14 February 2008

Abstract
Context. Linear and also nonlinear pulsation models with the variable Eddington factor approximation of the radiative transfer are constructed.
Aims. The aim of the study is to apply hydrodynamic models of radial pulsation to the observed variability in some post-AGB stars. It has been shown that the pulsation behavior could be strongly effected by the radiative field of the envelope of these stars because of high luminosity in these low mass supergiant stars. Thus, it is important to treat the radiative field with a higher approximation than the diffusion approximation which has been used successfully in classical Cepheids.
Methods. The moment equations of radiative transfer are integrated into the hydrodynamic equations with a variable Eddington factor. The factor is calculated independently by solving the transfer equation of spherical geometry. The linear eigen value problem of the radial perturbation from hydrostatic equilibrium is solved, and nonlinear simulations of radial pulsation are performed with the same approximation of the radiative transfer. The method is applied to pulsation of low mass supergiant stars and, for comparison, classical Cepheids.
Results. The properties of the strange modes that appear often in low mass supergiant star models e.g. post-AGB stars are seriously affected by different treatments of the radiative field in the linear analyses and nonlinear simulations.

Key words: stars: variables: general - radiative transfer

1 Introduction

The diffusion approximation of radiative transfer has been used for radial pulsation models of classical cepheids and RR Lyrae variables. The approximation is used for linear analysis and also nonlinear simulation (Castor 1971; Christy 1966a,b) for those stars. For much more luminous stars, however, it is necessary to use higher approximations of the radiative transfer.

The variable Eddington factor approximation for radiative transfer of spherical geometry (Mihalas & Mihalas 1984) is one candidate for this extension. Indeed, it was formulated for stellar pulsation models by Davis (1971) and Karp (1975), and has been applied to nonlinear radial pulsation for classical cepheids (Bendt & Davis 1971).

For these models, the Eddington factor, which is assumed as a constant in the the diffusion approximation, should be a function of the optical depth in plane and spherical geometry, and also is frequency dependent for non-gray atmospheres. The Eddington factor is determined by approximations of the original transfer equations (Unno 1976) or by directly solving the original transfer equation (York 1980).

Davis (1972) applied his radiative models to radial pulsation of W Virginis type stars. While there are no noticeable effects due to different treatments of the radiative transfer on Cepheid models, he pointed out that the diffusion approximation tends to dam up radiation in thin zones, not allowing it to flow freely. Shocks, propagating in thin zones have a tendency to be isothermal, radiating more strongly than they should be. His results suggest that these deficiencies related to shock regions are removed by using radiative transfer models. This comparative study strongly demonstrates the effects of different treatments of radiative transfer.

Fokin (1990) developed a nonlinear pulsation code of radiative transfer which was similar to the one by Davis (1972) and applied it to pulsation of low mass supergiant stars, and in particular, a post AGB star, HD 56126 (Jeannin et al. 1996, 1997).

On the other hand, Chistensen-Dalsgaard & Frandsen (1983) applied radiative transfer models with a variable Eddington factor approximation to a linear analysis of solar oscillations. Zalewski (1991, 1992) developed a linear analysis of radial pulsation with radiative transfer on supergiant stars, and applied the method to linear pulsation of the strange modes for post-AGB stars. He pointed out that the linear growth rates of so-called strange modes that appear with highly excited or damped modes of radial pulsation in linear models of post-AGB stars with the diffusion approximation (Wood 1976; Saio et al. 1984; Aikawa 1991, 1993) change when using radiative transfer for the same model parameters. Thus, comparative studies of different treatments of radiative transfer in linear analysis and nonlinear simulations are necessary.

In this paper, we present linear and nonlinear pulsation models using a unified treatment of radiative transfer. The method will be applied to the pulsation of low mass supergiant stars, particularly strange mode pulsation of post-AGB stars. Some post-AGB stars show photometric variability (see, Van Winckel ). Some of them have long-term observations, i.e. 89 Her and HD 161796 (Fernie 1986; Percy et al. 2000), and pulsation may be responsible for this variability (Lebzelter & Hinkle 2002). Aikawa (1991, 1993) showed that radial pulsation models for luminous low mass supergiant stars have strange modes and nonlinear simulations of these models show irregular oscillations with small amplitudes, which are distinguishing features of the variability of post-AGB stars. The existence of these irregular pulsations may be used to distinguish post-AGB stars from intermediate mass supergiant stars in the spectral types from B to G (Takeda et al. 2007). Moreover, pulsation is sometimes a key to understand the nature of post-AGB stars (Corsico et al. 2007) and post-AGB related stars, for instance, FG Sagittae (Jeffery & Schönberner 2006).

While we are concerned with the effects of radiative transfer on the pulsation behavior in low mass luminous stars, we also construct radiative transfer models for Cepheids for comparison.

2 Moment equations of radiative transfer

The basic hydrodynamic equations for radial pulsation and moment equations of radiative transfer are combined with a variable Eddington factor approximation.

We use the following fundamental equations of radiative transfer in spherical geometry (York 1980):

$\displaystyle \frac{1}{c}\frac{{\partial I^\nu }}{{\partial t}} + \frac{{(1 - \... ...gma _{\rm s}^\nu \int\limits_{ - 1}^1 {P(\mu ,\mu ')I^\nu } (\mu '){\rm d}\mu '$     (1)

where $I^\nu$ is the specific intensity along the ray with the direction cosine $\mu$ and photon frequency $\nu$. The absorption and scattering coefficients are $\sigma _{\rm a}^\nu$ and $\sigma _{\rm s}^\nu$, and the source function $S^\nu$, the phase function for photon scattering $P(\mu ,\mu ')$. Thomson scattering and Rayleigh scattering are included as absorption (Seaton 1993). Scattering due to dust particles should be included for extended atmospheres, but here we ignore this effect. The time dependent term is also ignored as it is of order v/c (Castor 1972).

We then use the method of moments for the radiative transfer equations to combine them with hydrodynamic equations. The resulting moments, to the second order, are defined as:

\begin{eqnarray*}M_0 &=& 2\pi \int\limits_{ - 1}^1 {I_\mu } {\rm d}\mu = E_\nu c... ... 1}^1 {I_\mu } \mu \mu {\rm d}\mu = \overline{\overline P}_\nu c \end{eqnarray*}


where $E_\nu$ is the radiation energy density, $\overline F_\nu$ the net radiation flux, $\overline{\overline P}_\nu$ the radiation pressure tensor, and c is the speed of light and the subscript $\nu$is photon frequency. When these moments are introduced into Eq. (1), by integration, we obtain two equations,

\begin{displaymath}\frac{1}{{r^2 }}\frac{\partial }{{\partial r}}(r^2 \overline ... ...sigma _{\rm T} \left(E_\nu - \frac{{4\pi B_\nu (T)}}{c}\right) \end{displaymath} (2)

and

\begin{displaymath}c\left[\frac{1}{{r^2 }}\frac{\partial }{{\partial r}}\left(r^... ...nu -E_\nu }{{r}} \right] = - \sigma _{\rm T} \overline F_\nu \end{displaymath} (3)

where $\sigma _{\rm T}=(\sigma '_{\rm a} + \sigma _{\rm s})$, and the first term of $\sigma_{\rm T}$ is the absorption coefficient corrected for induced emission, and the second term is due to Thomson scattering and Rayleigh scattering. The source function is reduced to $S = B_\nu(T)$ as LTE, where $B_\nu(T)$ is the Planck function.

Introducing the variable Eddington factor $f_\nu $ as

\begin{displaymath}f_\nu = \overline{\overline P}_\nu /E_\nu \end{displaymath} (4)

we can close Eqs. (2) and (3) as an equation system for two variables, $E_\nu$ and $\overline F_\nu$. The variable Eddington factor is estimated separately by numerically solving the original transfer equation with the given source function and opacities (York 1980), and may vary with limits between 1/3 and 1.

Finally, the Eqs. (2) and (3) are reduced to:

\begin{displaymath}\frac{1}{{r^2 }}\frac{\partial }{{\partial r}}(r^2 \overline ... ...sigma _{\rm T} \left(E_\nu - \frac{{4\pi B_\nu (T)}}{c}\right) \end{displaymath} (5)

and

\begin{displaymath}c\left[ {\frac{{\partial (f_\nu E_\nu )}}{{\partial r}} + \fr... ... - 1)}}{r}E_\nu} \right] = - \sigma _{\rm T} \overline F_\nu. \end{displaymath} (6)

Now, according to Mihalas & Mihalas (1984), we introduce two mean opacities (in units of $\rm cm^2/g$) to reduce the frequency-dependent Eqs. (5) and (6) to frequency-integrated moment equations:

\begin{displaymath}\frac{1}{{\kappa _R }} = \int_0^\infty {(\sigma _{\rm T} /\rh... ...\int\limits_0^\infty {(\partial B_\nu /\partial T)\rm d\nu } \end{displaymath} (7)


\begin{displaymath}\kappa _P = \int\limits_0^\infty {(\sigma _{\rm T} /\rho) B_\nu \rm d\nu } /\int\limits_0^\infty {B_\nu \rm d\nu }. \end{displaymath} (8)

The former opacity, the Rosseland mean, is introduced as the mean opacity which guarantees the correct radiative energy transport in the diffusion regime. The latter, the Planck mean, is introduced to obtain correct values for the total energy emitted or absorbed by the material.

Finally, after dividing by $\rho$ and integrating over the photon frequency, we obtain two equations which may be coupled with hydrodynamic equations:

\begin{displaymath}\frac{{\partial L_{r} }}{{\partial M_{r}}} = \kappa _P \left (4\sigma T^4 - ck_{E}E\right) \end{displaymath} (9)


\begin{displaymath}\frac{{\partial (fE)}}{{\partial M_{r}}} = - \frac{{\kappa _... ...i r^2 )^2c }}k_{F}L_{r} - \frac{1}{{4\pi r^3 \rho }}(3f - 1)E \end{displaymath} (10)

where

\begin{displaymath}L_{r} = 4\pi r^2 \int_0^\infty {\overline {F}_\nu \rm d\nu } \end{displaymath} (11)


\begin{displaymath}E = \int_0^\infty {E_\nu \rm d\nu } \end{displaymath} (12)


\begin{displaymath}f = \int_0^\infty {f_\nu E_\nu \rm d\nu } /\int_0^\infty {E_\nu \rm d\nu } \end{displaymath} (13)


\begin{displaymath}k_{\rm F} = \kappa _F /\kappa _{\rm R} = \int_0^\infty {(\sig... ...nu\rm d\nu } /\kappa _{\rm R} \int_0^\infty {F_\nu \rm d\nu } \end{displaymath} (14)


\begin{displaymath}k_{\rm E} = \kappa _E /\kappa _{\rm P} = \int_0^\infty {(\sig... ... \rm d\nu } /\kappa _{\rm P} \int_0^\infty {E_\nu \rm d\nu }. \end{displaymath} (15)

Here, the mass coordinate Mr is used instead of r, and $\sigma$ is the Stefan-Boltzmann constant which appears with $\int _0^\infty {B_\nu \rm d\nu}=\sigma T^4/\pi$. $k_{\rm F}$ is the ratio of the flux mean to the Rosseland mean, and $k_{\rm E}$ is the ratio of the absorption mean to the Planck mean. The denominators of the two ratios are calculated using the Planck function as the spectral density function. But the other two mean opacities depend on the full nongrey radiation field. For the first approximation, we assume that the flux mean is equal to the Rosseland mean, and the absorption mean is equal to the Planck mean, i.e. $k_{\rm F}=1$ and $k_{\rm E}=1$ (see Mihalas & Mihalas 1984, for discussions of this). f is the frequency-mean Eddington factor.

Under the hydrostatic equilibrium i.e. $L_{r} = L_{\rm T} = \rm const.$ and the Eddington approximation f=1/3, the above equations are reduced to the well-known expression of the diffusion approximation:

\begin{displaymath}L_{r} = - \frac{{\left(4\pi r^2\right)^2 c}}{{\kappa _R }}\fr... ...{\partial M_{r}}}\left(\frac{{4\sigma T^4 }}{{3c}}\right)\cdot \end{displaymath} (16)

The above expression is obtained with the Eddington approximation f=1/3 and hydrostatic equilibrium.

We introduce K = fE, and then we have

\begin{displaymath}\frac{{\partial K}}{{\partial M_{r}}} = - \frac{{\kappa _R }}... ...2 c }}L_{r} - \frac{1}{{4\pi r^3 \rho }}\frac{{(3f - 1)}}{f}K \end{displaymath} (17)


\begin{displaymath}\frac{{\partial L_{r} }}{{\partial M_{r}}} = \kappa _P \left(4\sigma T^4 - \frac{c}{f}K\right). \end{displaymath} (18)

The boundary conditions for these equations are:

(A) at the bottom of the envelope,

\begin{displaymath}L _{r} = L_{\rm T}=\rm const. \end{displaymath} (19)

(B) at the stellar surface with the optical depth $\tau _0$,

\begin{displaymath}K = \frac{{L_{r} }}{{4\pi r^2 c}}(\tau _0 + q_0 ) \end{displaymath} (20)

where K, L r and r are evaluated at the the optical depth $\tau _0$, and q0 is the Hopf function, evaluated at the optical depth $\tau _0$ but we assume a constant, $1/\sqrt{3}$.

Substituting K in Eqs. (18) to (17) with f=1/3, we can estimate the difference between the non-equilibrium and equilibrium diffusion approximation:

$\displaystyle L_{r} = - \frac{{(4\pi r^2 )^2 c}}{{\kappa _R }}\frac{\partial }{... ...ft(\frac{1}{{3c\kappa _P }}\frac{\partial }{{\partial M_{r}}}L_{r} \right)\cdot$     (21)

The first term on the right hand side is the expression of L r in the diffusion approximation, and so the second term, which consists of the second derivative of L r, yields the difference. Several numerical experiments on pulsation behavior show that the differences are very small.

To evaluate the Eddington factor f, we need to solve the radiative transfer equations with the given source function and the opacities (Yorke 1980). Again assuming f=1/3, we obtain radiative hydrodynamics with a non-equilibrium diffusion approximation. The approximation was used in the models by Dorfi & Feuchtinger (1991). In this paper, instead, we numerically solve the radiative transfer equations for a spherical geometry to evaluate the factor. We assume a grey atmosphere with the opacity of a material that may be estimated with the Rosseland mean opacity. We assume that the Eddington factor obtained from the grey atmosphere is equal to f in Eqs. (17) and (18).

3 Radiative hydrodynamics

We start with the following hydrodynamic equations for radial motion of material in the stellar envelope. The spatial coordinate is the mass within a distance r from the center, M r. So we have the Lagrangian expression of the hydrodynamic behavior. The dependent variables in these equations are the radial coordinate r, the flow velocity u, the specific volume V, the temperature T.

\begin{displaymath}\left( {\frac{{\partial r}}{{\partial t}}} \right)_{M_{r} } = u \end{displaymath} (22)


\begin{displaymath}\left( {\frac{{\partial u}}{{\partial t}}} \right)_{M_{r} } =... ... \left( {\frac{{\partial P}}{{\partial M_{r} }}} \right)_{t} \end{displaymath} (23)


\begin{displaymath}\left( {\frac{{\partial E}}{{\partial t}}} \right)_{M_{r} } +... ...eft( {\frac{{\partial L_{r} }}{{\partial M_{r} }}} \right)_t \end{displaymath} (24)


\begin{displaymath}V = \left( {\frac{\partial }{{\partial M_{r} }}} \right)_t \left( {\frac{{4\pi }}{3}r^3 } \right) \end{displaymath} (25)

where G is the gravitation constant, and the internal energy per unit mass, E and the pressure P are related to the temperature, T and the specific volume V by the equation of state, i.e. E=E(T,V) and P=P(T,V). The variations of these quantities are obtained from the reciprocity relations:
$\displaystyle \delta E = C_V \delta T + PV(\chi _{\rm T} - 1)(\delta V/V)$     (26)
$\displaystyle \delta P = (P\chi _{\rm T} /T)\delta T - (P\chi _\rho )(\delta V/V).$     (27)

Through the energy conservation Eq. (24), the material is coupled with the radiative field.

The boundary conditions for these equations are:

at the bottom of the envelope(C)


\begin{displaymath}r={\rm const.}, u = 0 \end{displaymath} (28)

and at the surface (D)


\begin{displaymath}P = P_{r} = (q_0 /c)L_{r} /(4\pi r^2 ) \end{displaymath} (29)

i.e. we assume that the pressure at the outermost layer is the radiative pressure evaluated at the outer boundary.

3.1 Model parameters

In the following, we study radial pulsation of low mass supergiant stars, in particular, pulsation due to the strange mode. We compare this with pulsation in Cepheid models to emphasize the effects on the strange modes. We summarize the model parameters of the two models for a low mass supergiant star and a Cepheid.

For a low mass supergiant (LmSG) star,


\begin{eqnarray*}&& M = 0.8~M_ \odot \nonumber \\ && L = 3000~L_ \odot \nonumbe... ...m e} = 6300~\rm K \nonumber \\ && X = 0.700,Z = 0.002 \nonumber \end{eqnarray*}


and for a Cepheid,

\begin{eqnarray*}&& M = 5.0~M_ \odot \nonumber \\ && L = 4369~L_ \odot \nonumbe... ...} = 5800~{\rm K} \nonumber \\ && X = 0.700,Z = 0.020. \nonumber \end{eqnarray*}


In these models, we ignored the effect of convection, and we used OP opacity (Seaton 1993; Seaton et al. 1994) for Roseland mean and Planck mean opacities. For regions of low density and low temperature, the opacities supplied by Alexander & Ferguson (1994) are used. We use a simplified equation of state: chemical abundances consist of hydrogen (X), helium (Y) and metal abundances (Z). The metal abundances are assumed to be in solar abundance ratios. Then Na, Al are assumed always ionized, and Mg, Si, Fe are treated as a single element and the first ionization is considered. Other elements are ignored.

3.2 Hydrostatic equilibrium

With Lr being $L_{\rm T}$, a constant in the envelope, and from Eq. (18), we obtain

\begin{displaymath}K=4\sigma fT^4 /c. \end{displaymath} (30)

Substituting this into Eq. (17), we obtain an equation for the temperature distribution in the envelope for a given distribution of the Eddington factor, i.e.,
$\displaystyle \frac{{\kappa _R L}}{{(4\pi r^2 )^2 }} = - \frac{1}{q}\frac{\partial }{{\partial M_{r} }}\left(qf4\sigma T^4\right)$     (31)
$\displaystyle {\rm {where }}\ln q = \int\limits_{M_{\rm c} }^{M_{r} } {\frac{{(3f - 1)}}{{4\pi r^3 f\rho }}} {\rm d}m$      

where $M_{\rm c}$ is the mass coordinate of the bottom of the envelope. We can integrate Eq. (31) starting with the boundary condition at the surface (B) in which we assume the surface as the location of the optical depth, $\tau _0 = 0.001$.

Combining this equation with the hydrostatic equilibrium, we can determine the temperature and density distribution in the envelope for a given distribution of the Eddington factor f. Then, we numerically solve the transfer equation for a given source function and opacity distribution. We iterate this process until the iteration reaches a consistent solution starting with the Eddington approximation, i.e., f=1/3.

Figure 1 show the distributions of the factor in the envelope for both the models. Compared with the Cepheid model, the LmSG model has a much wider region where the Eddington factor deviates from f=1/3, and the region reaches just above the top of the hydrogen ionization zone. This feature is important for the behavior of the strange mode pulsation, because the hydrogen ionization zone is responsible for pulsation driving of the mode (Aikawa & Sreenivasan 1985).


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496fig1.ps}} \end{figure} Figure 1: The distribution of the Eddington factor in the hydrostatic equilibrium for the models of low mass supergiant stars (thick solid line) and Cepheids (dotted line). The zone of the hydrogen ionization for the former model is indicated with horizontal line.
Open with DEXTER


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496fig2.ps}} \end{figure} Figure 2: The limb darkening of LmSG (thick solid line) and Cepheid (dotted line) models and the Eddington approximation (solid line).
Open with DEXTER

The distribution of the Eddington factor in the envelope reflects the limb darkening law at the surface. The Eddington approximation f=1/3 gives the well-known limb darkening law:

\begin{displaymath}I(0,\mu )/I(0,1) = \frac{3}{5}\left( {\mu + \frac{2}{3}} \right) \end{displaymath} (32)

where $I(\tau, \mu)$ is the intensity at the optical depth, $\tau$and the direction cosine $\theta$ and $\mu = \cos (\theta)$.

Figure 2 shows the limb-darkening of the present two models. Except for extreme limb regions, the limb darkening laws of both models coincide with that reached using by the Eddington approximation. The deviation from the Eddington approximation is much greater in LmSGstar model than in the Cepheid model, as expected.

In the following $f_{\rm eq}$ denotes the Eddington factor in equilibrium models.

3.3 Linear pulsation

For linear theory, we combine Eqs. (17) and (18) to one equation:

\begin{displaymath}\left( {\frac{{\partial ^2 r}}{{\partial t^2 }}} \right)_{M_{... ...2 \left( {\frac{{\partial P}}{{\partial M_{r} }}} \right)_{t} \end{displaymath} (33)

and the equation of energy conservation is expressed with the change of the entropy (S):

\begin{displaymath}\left( {T\frac{{\partial S}}{{\partial t}}} \right)_{M_{r} } ... ...\frac{{\partial L_{r} }}{{\partial M_{r} }}} \right) _{t}\cdot \end{displaymath} (34)

Now we convert the equations for linear pulsation to a finite-difference form in the mass coordinate, Mr. The spatial zoning is indicated by a subscript, which is a half-odd integer for zonal quantities such as, P, T, K, S, f, etc., and an integer for interface quantities such as r, M r, L r, etc. The index increases with radius. i=1 is for the inner-most interface and i=(N+1) is for the outer-most interface.
$\displaystyle \frac{{{\rm d}^2 r_i }}{{{\rm d}t^2 }} = - \frac{{G(M _{r})_i }}{{r_i^2 }} - 4\pi r_i^2 \frac{{(P_{i + 1/2} - P_{i - 1 + 1/2} )}}{{Dm2_i }}$     (35)
i = 2,3, ...N,N + 1      


$\displaystyle (L_{r} )_i = - {\frac{{(4\pi r_i^2 )^2 c}}{{(\kappa _R )_i }}} \frac{{(K_{i + 1/2} - K_{i - 1 + 1/2} )}}{{Dm2_i }}$      
$\displaystyle -{\frac{{(4\pi r_i^2 )^2 c}}{{(\kappa _R )_i }}} {\left[ \frac{1}{{4\pi r^3 \rho }}\frac{{(3f - 1)}}{f}K\right]}_i$     (36)
i=2,3,... N      


$\displaystyle K_{i + 1/2} = \frac{{f_{i + 1/2} }}{c}4\sigma T_{i + 1/2}^4 - \fr... ...2} }}\left( {\frac{{(L_{r} )_{i + 1} - (L_{r} )_i }}{{Dm1_{i + 1/2} }}} \right)$     (37)
i=1,2,... ,N      


$\displaystyle \left( {T\frac{{{\rm d}S}}{{{\rm d}t}}} \right)_{i + 1/2} = - \frac{f_{i+1/2}}{{Dm1_{i + 1/2} }}\left\{ {(L_{r} )_{i + 1} - (L_{r} )_i } \right\}$     (38)
i=1,2,... ,N      

where

Dm1 i+1/2 = (M r) i+1-(M r) i (39)


Dm2 i = (Dm1 i+1/2+Dm1 i-1/2)/2. (40)

The boundary conditions for these equations are:

At the inner boundary:

\begin{displaymath}r_1 = \mbox{const. (fixed value)}. \end{displaymath} (41)


\begin{displaymath}(L_{r})_1 = L_{\rm T}= \mbox{const. (fixed value).} \end{displaymath} (42)

At the surface:

\begin{displaymath}(L_{r})_{N+1} = K_{N+1/2}4\pi r^2c/(\tau _0 + q_0) \end{displaymath} (43)


\begin{displaymath}P_{N+1+1/2} = (q_0/c)(L_{r})_{N+1}/\left(4\pi r_{N+1}^2\right). \end{displaymath} (44)

We assume in the last expression that the pressure at the outside of the outer-most interface is equal to the radiation pressure evaluated at the location.

For a linear analysis we assume that we are given a model that is in hydrostatic equilibrium according to the difference Eqs. (35)-(38) with time-derivatives set to zero. The infinitesimal deviation of any variables from their values in the hydrostatic equilibrium will be indicated by the prefix $\delta $. We also need  $f_{\rm osc}$, the Eddington factor for the perturbation, which might be evaluated by the perturbation on radiative intensity $\delta I _{r}$. We linearize Eqs. (35)-(38) by expanding all the functions appearing in them in a Taylor series about the equilibrium model, retaining only terms of zero or first order in the perturbation. Thus we replace the first order variables as follows:

\begin{displaymath}x_i = (Dm2_i)^{1/2}\delta r_i \end{displaymath} (45)


\begin{displaymath}l_i = \delta (L_{r})_i/(L_{r})_i \end{displaymath} (46)


\begin{displaymath}k_{i+1/2} = \delta K_{i+1/2}/K_{i+1/2} \end{displaymath} (47)


\begin{displaymath}y_{i+1/2} = T_{i+1/2}\delta S_{i+1/2}. \end{displaymath} (48)

Then we have (4N-1) equations with (4N-1) variables, which are arranged with the following vector:
$\displaystyle \vec{V}=(k_{1+1/2},y_{1+1/2},x_2,l_2,..x_i,l_i,k_{i+1/2},y_{i+1/2},..$      
..xN,lN,kN+1/2,yN+1/2,xN+1)T.     (49)

We search for the normal modes of radial pulsation by assuming that the time dependence is exponential. Therefore we have a factor $\exp (i\omega t)$ for all the perturbations and replace ${\rm d}/{\rm d}t$ by $i\omega $. The eigen frequency $\omega $ is complex in general.

For these variables, we have the following equations:

                          $\displaystyle \omega ^2 x_i$ = G11,i xi - 1 + G12,i xi + G13,i xi + 1  
    + G21,i yi - 1 + 1/2+G22,i yi + 1/2 (50)
    i = 2,3,... ,N,N + 1  


                         li = $\displaystyle - \frac{(4\pi r^2 )^2 c}{(L_{r} )_i (\kappa _R )_ i}\delta F_i + ... ... c}{(L_{r} )_i (\kappa _R )_ i } \bigl[ { BL1_{1,i} x_{i - 1} + BL1_{2,i} x_i }$  
    $\displaystyle +{ BL1_{3,i} x_{i + 1} + BL2_{1,i} y_{i - 1+ 1/2} + BL2_{2,i} y_{i + 1/2} } \bigr]$ (51)
    i = 2,3,... N  


                         ki + 1/2 = $\displaystyle \frac{f_{i + 1/2} }{K_{i + 1/2} c}4(4\sigma )T_{i + 1/2}^4 \biggr[ {\frac{y_{i + 1/2} }{(c_V T)_{i + 1/2} } + (\Gamma _3 - 1)_{i + 1/2}}$  
    $\displaystyle \times {({\rm d}R1_{i + 1/2} x_i + {\rm d}R2_{i + 1/2}x_{i + 1} )} \biggr] + \frac{1}{cK_{i + 1/2} (\kappa _{\rm P})_{i +1/2} }$  
    $\displaystyle \times \left[ {\frac{(L_{r} )_{i + 1} }{Dm1_{i + 1/2} }l_{i + 1} - \frac{(L_{r} )_i }{Dm1_{i + 1/2} }l_i } \right]$ (52)
    i = 1,2,... ,N  


$\displaystyle -i\omega y_{i + 1/2} = \frac{{(L_{r} )_i }}{{Dm1_{i + 1/2} }}l_i - \frac{{(L_{r} )_{i + 1} }}{{Dm1_{i + 1/2} }}l_{i + 1}$     (53)
i = 1,2,... ,N.      

In Eq. (46),
$\displaystyle F = \frac{{(K_{i + 1/2} - K_{i - 1 + 1/2} )}}{{{\rm d}m2_i }} + \left[ {\frac{1}{{4\pi r^3 \rho }}\frac{{(3f - 1)}}{f}K} \right]_i$     (54)

and $\delta F $ is the linearly perturbed part of F.

Then we have the eigen value problem with a matrix, M, and the vector, $\vec{V}$:

\begin{displaymath}\vec{M(\omega )}\cdot \vec{V} = 0 \end{displaymath} (55)

where $\vec{M(\omega )}$ is a band matrix in which the elements are functions of quantities in hydrostatic equilibrium and the eigen frequency, $\omega $. The band width of this band matrix is 10, if we arrange linearized Eqs. (50)-(53) in order.

We solve this eigen value problem by iterations. For Cepheids, the eigen values estimated by the quasi-adiabatic approximation will be appropriate for the initial guess of the eigen values. For a luminous low mass supergiant star, however, non-adiabatic effects are so strong that the initial guess of eigen values is useless. Therefore, we adopt the following strategy.

We introduce a control parameter $\epsilon $ in Eq. (38):

              $\displaystyle \left( {T\frac{{{\rm d}S}}{{{\rm d}t}}} \right)_{i + 1/2}$ = $\displaystyle - \epsilon \frac{1}{{Dm1_{i + 1/2} }}\left[ {(L_{r} )_{i + 1} - (L_{r} )_i } \right]$ (56)
    i=1,2,... ,N.  

The parameter, which controls the adiabaticity of perturbations, varies between 0 and 1. The special case, $\epsilon = 0$, corresponds to adiabatic perturbation, and $\epsilon = 1$ corresponds to fully non-adiabatic perturbation.


  \begin{figure} \par\resizebox{8.5cm}{!}{\includegraphics{8496fig3.ps}} \end{figure} Figure 3: The pulsation periods (days) and the growth rates as a function of $\epsilon $ for the Cepheid model.
Open with DEXTER

Equation (55) then is rewritten logically as

\begin{displaymath}\vec{M(\omega ,\epsilon )}\cdot \vec{V} = 0. \end{displaymath} (57)

To solve this eigen value problem, we at first set $\epsilon $to 0. For this special case the eigen value problem is complete and this means we will derive a complete set of eigen values. We then gradually increase the value of $\epsilon $ to 1. We therefore obtain eigen values that have counterparts in adiabatic perturbation.

For the Eddington factor, we assume at first $f=f_{\rm eq}$ in the linear pulsation. $f_{\rm eq}$ is the Eddington factor for the hydrostatic equilibrium. Figure 3 shows the pulsation period and the growth rates for lower modes as a function of $\epsilon $ in the Cepheid model. The difference of the present results to those obtained by the Eddington approximation is quite small.

Figure 4 shows the results for the LmSG star model. Compared with Fig. 5 obtained for the same model with the diffusion approximation, the positive value of the growth rate of the strange mode, the line of which is labeled 3 in Figs. 4 and 5, is remarkably reduced in this approximation.

As the next problem, we evaluate f in the presence of the linearly perturbed quantities. We start with the following equation for perturbed intensity $\delta I$ which is derived from the Lagrangian perturbed material quantities:

\begin{displaymath}\mu \frac{{\partial (\delta I)}}{{\partial m_{r} }} + \frac{1... ...- \delta \left[ {\frac{{\kappa (I - S)}}{{4\pi r^2 }}} \right] \end{displaymath} (58)

where mr is the mass coordinate, and the operation $\delta $ is the Lagrangian perturbation of the bracket. We ignore the terms associated with $ \frac{{\partial I_0}}{{\partial \mu }}$.
  \begin{figure} \par\resizebox{8.5cm}{!}{\includegraphics{8496fig4.ps}} \end{figure} Figure 4: The pulsation periods (days) and the growth rates as a function of $\epsilon $ for the LmSG star model.
Open with DEXTER


  \begin{figure} \par\resizebox{8.5cm}{!}{\includegraphics{8496fig5.ps}} \end{figure} Figure 5: The pulsation periods (days) and the growth rates as a function of $\epsilon $ for the LmSG star model obtained by using the diffusion approximation.
Open with DEXTER

Evaluating the bracket, followed by the r-coordinate, we have the equation for complex $\delta I$:


$\displaystyle \mu \frac{{\partial (\delta I)}}{{\partial r}} + \frac{{(1 - \mu ... ...(I - S)\frac{{\delta \kappa }}{\kappa } - 2(I - S)\frac{{\delta r}}{r}} \right]$     (59)

where I and S are the mean intensity and source function of the hydrostatic equilibrium model. Thus, we solve the radiative transfer equation of spherical geometry with the complex source function. The boundary conditions of this are the same as used for the equilibrium model.

From $\delta I$, we can evaluate the Eddington factor which is complex in general. We denote this  $f_{\rm osc}$. Then we solve the linear eigen value problem using the Eddington factor  $f = f_{\rm osc}$ to obtain the eigen value and the perturbations. To arrive at consistent solutions, we again require iterations. We solve the linear problem with an assumed  $f_{\rm osc}$. The starting form of the assumed  $f_{\rm osc}$ is $f=f_{\rm eq}$ and then we evaluate f with the perturbed quantities until a consistent solution is reached. For some cases, the iteration unexpectedly skips between different modes. Then we take the following strategy:

\begin{displaymath}f^{n+1}_{\rm osc} = (1- \varepsilon)f^{n-1}_{\rm osc} + \varepsilon f^{n}_{\rm osc} \end{displaymath} (60)

where $\varepsilon$ is a small factor, about 0.05. Starting with $f^{0}_{\rm osc} = f_{\rm eq}$ and  $f^{1}_{\rm osc}$ which is obtained from the perturbed quantities with $f=f_{\rm eq}$, we continue the iteration until $f^{n+1}_{\rm osc}=f^{n}_{\rm osc}$ within the small error. By this method, we successfully obtain consistent solutions of the eigen value problem.

Figure 6 shows the paths of the convergence of the iteration for the low mass supergiant star model. We obtain the eigen values up to the 10th modes. Except for the strange mode, the effects on the periods and the growth rates of using $f = f_{\rm osc}$ are small. For the strange mode, the present approximation has serious effects on the pulsation period as well as the growth rate. We demonstrate in Fig. 7 the distribution of $f_{\rm osc}$ for the strange mode in the LmSG model. At the hydrogen ionization zone, it deviates strongly from $f_{\rm eq}$, the value of the equilibrium model. We summarized the results of the linear models using various Eddington factor approximations in Table 1. Compared with the ordinary modes, the strange mode has remarkably strong changes for the growth rates and also the periods.


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496fig6.ps}} \end{figure} Figure 6: The paths of the convergence of iteration for the model with $f = f_{\rm osc}$ for the LmSG model. The paths of the modes, up to the 10th overtone, are indicated with the label of the mode. The ordinate is the real part of the eigen value, and the abscissa is the imaginary part.
Open with DEXTER


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496fig7.ps}} \end{figure} Figure 7: The distribution of the Eddington factor, $f_{\rm osc}$ (complex), for various modes for the LmSG model. The real and imaginary parts of the factor are denoted by solid and dotted lines, respectively.
Open with DEXTER

Table 1: The periods and the growth rates of the LmSG model.

3.4 Nonlinear pulsation

We start nonlinear pulsation models with expressions (17), (18) and (22)-(25). As described in Eqs. (35)-(38) as the finite difference expression, these equations may be converted directly to finite-difference forms. However, dynamic re-zoning, during which ionization zones and shock regions are confined to fine zoning regions may reduce the effects of artifacts in the light curves. So we introduce re-zoning as a diffusion equation on mesh points. The equations now become:

\begin{displaymath}\left( {\frac{{\partial r}}{{\partial t}}} \right)_{x} = u +... ...r} \left( {\frac{{\partial r}}{{\partial M_{r} }}} \right)_{t} \end{displaymath} (61)


\begin{displaymath}\left( {\frac{{\partial u}}{{\partial t}}} \right)_{x} = - \... ...M\left( {\frac{{\partial u}}{{\partial M_{r} }}} \right)_{t} \end{displaymath} (62)


\begin{displaymath}L_{r} = - \frac{{(4\pi r^2 )^2 c}}{{\kappa _R }}\left( {\fra... ...} + \frac{1}{{4\pi r^3 \rho }}\frac{{(3f - 1)}}{f}K} \right) \end{displaymath} (63)


\begin{displaymath}K = \frac{f}{c}4\sigma T^4 - \frac{f}{{c\kappa _P }}\left( {\frac{{\partial L_{r} }}{{\partial M_{r} }}} \right) \end{displaymath} (64)


\begin{displaymath}\left( {\frac{{\partial E}}{{\partial t}}} \right)_{x} + P\le... ...\frac{{\partial V}}{{\partial M_{r} }}} \right)_{t} } \right] \end{displaymath} (65)

where x is a new coordinate which is a function of M r and t. We shall use the coordinate to describe the zonal interfaces that are not fixed in mass, and $\dot M$ denotes $\left( {\frac{{\partial M_{r} }}{{\partial t}}} \right)_{x} $ which is calculated with algorithms for dynamic re-zoning.

The boundary conditions for these equations are the same in linear pulsation:

At the inner boundary:

\begin{displaymath}r_1 = \mbox{const. (fixed value)} \end{displaymath} (66)


u1 = 0 (67)


\begin{displaymath}(L_{r})_1 = L_{\rm T}= \mbox{const. (fixed value).} \end{displaymath} (68)

At the surface:

\begin{displaymath}(L_{r})_{N+1} = K_{N+1/2}4\pi r^2c/(\tau _0 + q_0) \end{displaymath} (69)


\begin{displaymath}P_{N+1+1/2} = (q_0/c)(L_{r})_{N+1}/\left(4\pi r_{N+1}^2\right) \end{displaymath} (70)

where again we assume that the pressure at the outer boundary is the radiation pressure at the boundary.

The finite difference scheme that has been used most successfully is presented in the following equations:

$\displaystyle r_i^{n + 1} - r_i^n = \frac{{(\Delta t)^{n + 1/2} }}{2}(u_i^{n + ... ...\left(r_{i + 1}^{n + 1} + r_{i + 1}^n - r_{i - 1}^{n + 1} - r_{i - 1}^n \right)$     (71)


                                       $\displaystyle u_i^{n + 1} - u_i^n = \frac{{GM_i^{n + 1/2} (\Delta t)^{n + 1/2} }}{{r_i^{n + 1} r_i^n }}$  
    $\displaystyle \qquad +\frac{{2\pi (\Delta t)^{n + 1/2} }}{{3Dm2_i^{n + 1} }}\left[ {(r_i^{n + 1} )^2 + r_i^{n + 1} r_i^n + (r_i^n )^2 } \right]$  
    $\displaystyle \qquad \times\left[ {(p + q)_{i - 1/2}^{n + 1} + (p + q)_{i - 1/2}^n - (p +q)_{i + 1/2}^{n + 1} - (p + q)_{i + 1/2}^n } \right]$  
    $\displaystyle \qquad +\frac{{\dot M_i^{n + 1} (\Delta t)^{n + 1/2} }}{{4Dm2_i^{n + 1} }} (u_{i + 1}^{n + 1} + u_{i + 1}^n - u_{i - 1}^{n + 1} - u_{i - 1}^n )$ (72)


$\displaystyle (L_{r})_i^{n + 1} = \left( {\frac{{(4\pi (r_i^{n + 1} )^2 )^2 c}}... ... {\frac{1}{{4\pi r^3 \rho }} \frac{{(3f - 1)}}{f}K} \right)_i^{n + 1} } \right]$     (73)


\begin{displaymath}K_{i + 1/2}^{n + 1} = \frac{{f_{i + 1/2}^{n + 1} }}{c}4\sigm... ... 1} - (L_{r} )_i^{n + 1} }}{{Dm1_{i + 1/2}^{n + 1} }}} \right) \end{displaymath} (74)


                                        $\displaystyle E_{i + 1/2}^{n + 1} - E_{i + 1/2}^n + \frac{1}{2}\left[ {(p + q)_{i + 1/2}^{n + 1} + (p + q)_{i + /1/2}^n } \right]$  
    $\displaystyle \qquad \times (V_{i + 1/2}^{n + 1} - V_{i + 1/2}^n ) = \frac{{\le... ...{r} )_{i + 1}^{n + 1} } \right)}}{{Dm1_{i + 1/2}^{n + 1} }}(\Delta t)^{n + 1/2}$  
    $\displaystyle \qquad + \frac{1}{2}\dot M_{i + 1/2}^{n + 1/2} (\Delta t)^{n + 1/... ...} + \frac{{E_{i + 1/2}^{n + 1} - E_{i - 1 + 1/2}^{n + 1} }}{{Dm2_i^{n + 1} }} }$  
    $\displaystyle \qquad{+ (p + q)_{i + 1/2}^{n + 1} \left( {\frac{{V_{i + 1 + 1/2}... ...+ 1/2}^{n + 1} - V_{i - 1 + 1/2}^{n + 1} }}{{Dm2_i^{n + 1} }}} \right)} \Biggr]$ (75)

where we introduce the artificial viscosity q to stabilize shock waves. The formula used is:

\begin{displaymath}q_{i + 1/2}^n = \frac{{C_{\rm q} }}{{V_{i + 1/2}^n }}\left[ {... ..._{\rm qcut} \sqrt {p_{i + 1/2}^n V_{i + 1/2}^n } )} \right]^2 \end{displaymath} (76)

where $C_{\rm q} $ is a parameter to control the amount of viscosity and  $C_{\rm qcut}$ is another parameter which prevents unnecessary damping in the deep interior of the envelope (Stellingwerf 1975).

Thus we have an implicit difference scheme which guarantees numerical stability without having to satisfy the Courant-Friedrichs-Lewy condition. The nonlinear equations are solved, as usual, by use of the Newton-Raphson technique. We expand the dependent vector variable as

\begin{displaymath}X_i^{n + 1} \to X_i^{n + 1} + \Delta X_i^{n + 1}. \end{displaymath} (77)

The resulting linear system for $\Delta X_i^{n + 1}$ has a banded matrix. The band width depends on the arrangement of variables. We adopt the arrangement of variables as follows:
$\displaystyle \Delta K_{1 + 1/2}^{n + 1} ,\Delta T_{1 + 1/2}^{n + 1} ,\Delta r_... ...\Delta L_2^{n + 1} ,\Delta K_{2 + 1/2}^{n + 1} ,\Delta T_{2 + 1/2}^{n + 1} ,...$      
$\displaystyle ...,\Delta r_i^{n + 1} ,\Delta u_i^{n + 1} ,\Delta L_i^{n + 1} ,\Delta K_{i + 1/2}^{n + 1} ,\Delta T_{i + 1/2}^{n + 1} ,...$      
$\displaystyle ...,\Delta r_N^{n + 1} ,\Delta u_N^{n + 1} ,\Delta L_N^{n + 1} , ... ...elta T_{N + 1/2}^{n + 1} , \Delta r_{N + 1}^{n + 1} ,\Delta u_{N + 1}^{n + 1} .$     (78)

We have (5N-1) variables with the same number of equations and the band width is 16.

To determine $\dot M$, we follow Castor et al. (1977) and Aikawa & Simon (1985) to keep a certain feature of the calculation stationary with respect to the mesh. We concentrate on ionization zones, because the zone should remain in the fine zoning during pulsation. Otherwise we will have many artifacts in the light curves that we should compare with observation data. Rezoning is defined by the equation:

\begin{displaymath}\left( {\frac{{\partial M_{r} }}{{\partial t}}} \right)_x = ... ...ft[ {Y\frac{{\partial M_{r} }}{{\partial x}}} \right]_{t} + A. \end{displaymath} (79)

In this equation B, Y, and A are adjustable functions that best meet the above objectives.

This equation expresses the diffusive motion of the mesh in the mass coordinate, and is expressed in the finite difference form:


$\displaystyle M_i^{n + 1} = M_i^n + B_i^n \Delta t^{n + 1}\left[ Y_{i + /12}^n ... ...{i - 1/2}^n (M_i^{n + 1} - M_{i - 1}^{n + 1} ) \right] +A_i^n \Delta t^{n + 1}.$     (80)

This tri-diagonal matrix for Min + 1 is solved quickly with a given Min, Ain, Bin, and Yi+1/2n. Thus we obtain the instantaneous mass coordinate at the intermediate time between n and n+1 time steps.

The simulation of the initial value problem is performed using the values of the static equilibrium except for the velocity profile. Moreover, this time we use fixed values of the parameters for the artificial viscosity, $C_{\rm q} = 4.0, C_{\rm qcut} = 0.02$. Then the simulation continues until we obtain stationary states of pulsation, which are considered as observed states of pulsating stars.

The non-linear simulations are performed at first with $f=f_{\rm eq}$, i.e, the Eddington factor in the non-linear models equal to the one in the hydrostatic equilibrium. Figure 8 shows light curves at the photosphere for the Cepheid models. The simulation is started with velocity profile of the fundamental mode with a scale factor of -10 km s-1 at the surface, and the model shows a limit cycle oscillation after a run of the time interval of about 10 000 days. The results obtained by using the diffusion approximation are also included. The light curves obtained without re-zoning, or with a variable Eddington factor approximation using snapshot values, $f_{\rm dyn}$ computed during the simulation, give only minor changes to the results shown in Fig. 8.

  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496fig8.ps}} \end{figure} Figure 8: The light curves at the photosphere of the Cepheid model by the variable Eddington factor approximation (sold line) and the diffusion approximation (dotted line).
Open with DEXTER

For the low mass supergiant star model, we first perform the simulations with the Eddington factor  $f_{\rm eq}$. We also used snapshot values of the factor evaluated at each time step during the simulation. We denote the factor as  $f_{\rm dyn}$ for this case. We believe that the simulation using  $f_{\rm dyn}$ may correspond to the linear model with  $f_{\rm osc}$. Figure 9 shows the model light curves obtained with different assumptions on the Eddington factor. It is noted that the differences are quite effective in changing the light curve of pulsation due to the strange mode. Figure 10 shows the periodogram of the light curves. It is shown that the light curve obtained by using the dynamic variable Eddington factor approximation is sinusoidal with a main period of 3.9 days and there are no higher components beside the primary component. On the other hand, the other two light curves are more complicated. For the model with $f=f_{\rm eq}$, the period of the main component is 6.5 days and there are some peaks in their higher harmonics. For the mode with f=1/3, it seems chaotic pulsation previals, although there are some peaks at the periods of about 6 days.


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496fig9.ps}} \end{figure} Figure 9: The light curves at the photosphere of the LmSG star model obtained using the variable Eddington factor approximations: the diffusion approximation ( top), $f_{\rm eq}$ ( middle) and $f_{\rm dyn}$ ( bottom).
Open with DEXTER


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{8496figa.ps}} \end{figure} Figure 10: The periodograms of the model light curves. The light curve obtained using the dynamic variable Eddington factor approximation ( $f_{\rm dyn}$) is shown as a solid line. Compared with those obtained by the static variable Eddington factor approximation ( $f_{\rm eq}$) (dashed line), and by the diffusion approximation(dotted line), the light curve obtained by  $f_{\rm dyn}$ is quite simple.
Open with DEXTER

4 Conclusions

1.
Linear and nonlinear radial pulsation models are constructed with a unified treatment of radiative transfer. They are applied to the pulsation of low mass supergiant stars and also a Cepheid.
2.
Higher approximations than the diffusion approximation for the radiative transfer cause only small changes in the Cepheid model.
3.
The models are applied to the strange modes which appear in the pulsation of post-AGB stars and are assumed to be responsible for the photometric variability of these stars. The pulsation behaviours are strongly affected by different treatments of the radiative transfer in the pulsation model. This is important for astroseismological studies using pulsation in the post-AGB stars.
4.
According to Aikawa (1991, 1993), the limit cycles appear in models with small values of luminosity, and chaotic pulsation due to the strange mode will appear in more luminous models. Thus, pulsation behaviours may be a function of luminosity. This factor, combined with the results if Figs. 9 and 10, suggest that chaotic pulsations, which are commonly observed in post-AGB stars, will require much higher luminosities in the  $f_{\rm dyn}$ models than in the diffusion models.

Acknowledgements
Part of this work was supported by the Japanese a Grant-in-Aid for Scientific Research of the Ministry of Education, Culture, Sports, science and Technology project number 14540229.

References

 
Copyright ESO 2008