2 Physical description

 In order to describe the nonlinear behavior of pulsating stars we assume the stellar plasma to consist of three components, namely gas, radiation and turbulent convection. The evolution of each component is governed by respective conservation laws and coupling terms between these equations account for their interaction. The resulting set of coupled equations called the system of convective radiation hydrodynamics (CRHD) consists of 7 nonlinear partial differential equations which have to be solved under appropriate initial and boundary conditions. All equations are written in spherical symmetry, meaning that each physical quantity is considered to be a function of radius r and time t. Consequently we are restricted to radial pulsations.

2.1 Gas dynamics

 The dynamics of the stellar gas is governed by three conservation laws for mass (continuity equation, Eq. (2) in Table 1), momentum (equation of motion, Eq. (3) in Table 1) and internal energy (gas energy equation, Eq. (4) in Table 1). The independent variables are the gas density $\rho$, the gas velocity u and the gas internal energy e. In contrast to most nonlinear pulsation models we use an equation for the internal energy instead of the total energy equation. A detailed discussion of implications and consequences for stellar pulsation models is given in Feuchtinger & Dorfi (1994, 1997). Additionally, as we consider self-gravitating flows, we solve the Poisson equation for the gravitational field. This is done by analytically integrating the Poisson equation in spherical symmetry yielding the gravitational force which is considered as a source term on the right-hand side of the equation of motion. The remaining integration of the inner mass $m\rm _r$ over the density structure is performed by solving Eq. (1) in Table 1, where we employ $m\rm _r$ as the independent variable.

Further details concerning the equations of hydrodynamics can be found e.g. in Mihalas & Mihalas (1984, Chapter 7).

  

Table 1: Equations of convective radiation hydrodynamics in differential form. D/Dt denotes the comoving or substantial derivative with respect to time. Refer to Table 6 for a comprehensive list of symbols

2.2 Radiative transfer

 Time-dependent radiative transfer is considered in the grey approximation by solving the first two frequency-integrated moment equations of the radiation field representing energy conservation (radiation energy equation, Eq. (5) in Table 1) and momentum conservation (radiation momentum equation, Eq. (6) in Table 1). As independent variables we use the zeroth moment of the radiation field J and the first moment of the radiation field H, which are proportional to the radiative energy density and the radiative flux, respectively (e.g. Mihalas & Mihalas 1984, page 337).

A special remark concerns the dynamical terms in the radiation momentum equation. These terms account for time-dependent effects in the radiation field where the light travel time becomes comparable to the dynamical time scale of the problem (e.g. traveling light fronts). Such phenomena usually are not encountered in typical stellar pulsation environments and consequently these terms can be neglected.

In contrast to that the dynamical terms in the radiation energy equation are important whenever the ratio u/c becomes comparable to the ratio $\lambda_{\rm p}/l_{\rm hyd}$ with $\lambda \rm _p$ and $l_{\rm hyd}$denoting the photon mean free path and the typical hydrodynamical length scale of the problem (Mihalas & Mihalas 1984, page 457). In this case (known as the dynamical diffusion regime) the advection of radiation energy by the moving fluid becomes important and such situations are likely to be encountered in stellar pulsation flows. In fact test calculations carried out with and without including the dynamical terms in the radiation energy equation exhibit differences in the pulsation properties like pulsation amplitudes and light curve shapes. Detailed results and further discussions will be presented in a separate paper.

2.3 Convective transfer

 Radiation hydrodynamical models allow a self-consistent description of pulsating envelopes as long as convection is negligible. This means that a radiative pulsating envelope is determined completely by specifying a set of stellar parameters (mass, luminosity, effective temperature and chemical composition). However, it turned out that several properties of classically pulsating stars (Cepheids, RR Lyrae) like the pulsation amplitudes, the cool edge of the instability strip or double mode pulsations (Kolláth et al. 1998) cannot be explained with purely radiative models. In fact energy transfer by turbulent convection plays an essential role in the physics of stellar oscillations and therefore cannot be neglected (e.g. Yecko et al. 1998).

Due to the lack of a self-consistent convection theory the inclusion of turbulent convection into the radiation hydrodynamical model introduces a number of free parameters. In principle these parameters cannot be determined from theoretical considerations, but have to be chosen in order to fit observational constraints. In that sense a self-consistent modelling of stellar pulsations must await considerable progress in turbulence theory.

In the past decades a number of time-dependent one-dimensional models of turbulent convection have been developed. Most of these models are based on the Boussinesq-approximation which assumes the eddy gas to be incompressible except for buoyancy effects. Turbulent motions are considered as fluctuations of the mean hydrodynamic field and this assumption is used to linearize the equations in the fluctuating parts (see Baker 1987 for an overview). In the context of stellar pulsations the majority of investigations use the convection model of Stellingwerf (1982a) to compute nonlinear limit cycle solutions for RR Lyrae and Cepheid pulsations (Stellingwerf 1982b; Gehmeyr 1992; Bono & Stellingwerf 1994; Bono et al. 1997; Bono et al. 1998).

For the nonlinear pulsation model presented in this paper we use the one-equation model of time-dependent turbulent convection proposed by Kuhfuß (1986) in the version of Wuchterl & Feuchtinger (1998). Although the approach is similar to the Stellingwerf-model there are some important differences discussed in Kuhfuß (1987) and Gehmeyr & Winkler (1992). Recently a similar model was used for linear and nonlinear Cepheid pulsations by Yecko et al. (1998) and Kolláth et al. (1998).

The basic equation (Eq. (8) in Table 1) of the Kuhfuß model is a conservation law for the turbulent kinetic energy density $\bar{\omega}$, which serves as the independent variable of the turbulent field. The essential term is the turbulent driving through buoyancy forces ($S_{\bar{\omega}}$). The basic quantity entering $S_{\bar{\omega}}$ is the entropy gradient $\partial s/\partial r$ which is related to the Schwarzschild criterion by Eq. (12) (Table 1). Note that $S_{\bar{\omega}}$ can be evaluated even if the Schwarzschild criterion indicates convective stability. This means that we do not have to neglect damping buoyancy forces as e.g. done in the Stellingwerf model (Stellingwerf 1982a). An important technical advantage of the adopted formulation for $S_{\bar{\omega}}$ is that we do not have to use a seed for the turbulence field.

Another term important for stellar pulsations is the energy dissipation $E\rm _Q$ through turbulent viscosity. A corresponding term accounting for the viscous momentum generation ($U\rm _Q$) enters the equation of motion. While radiative pulsation models either exhibit too large oscillation amplitudes or have to include artificial dissipation effects (e.g. Feuchtinger & Dorfi 1994) the turbulent viscosity provides a physically motivated dissipation, however, adjustable by a free parameter.

For further details about the convection model and its parameters we refer to Wuchterl (1995) and Wuchterl & Feuchtinger (1998).

2.4 Constitutive relations

 For closing the CRHD equations we need several constitutive relations. They specify the microscopic properties of the gas and the radiation field and are discussed in the following sections. Closure relations for the turbulent convection are implicitly included in the assumptions for modelling the different correlations terms and their physical meaning is outlined in Kuhfuß (1987).

2.4.1 Equation of state

The equation of state specifies pressure $P(\rho,e)$, temperature $T(\rho,e)$, isentropic temperature gradient $\nabla_{\rm s}(\rho,e)$ and specific heat at constant pressure $c_{\rm P}(\rho,e)$ of the gas with respect to the independent variables $\rho$ and e. Note that due to considering the radiation field as an extra component we employ the pure gas pressure instead of the total pressure in the equation of motion ($\partial P/\partial r$), while the radiation pressure is treated as an extra term ($4\pi/c \, \kappa_{\!H} \rho H$). This difference is essential in the case of anisotropic radiation fields (variable Eddington factor, see Sect. 2.4.3) where the radiation pressure tensor does not reduce to a simple scalar pressure.

The equation of state is included in tabulated form and it is essential to use suitable interpolation schemes to obtain smooth surfaces and derivatives (Dorfi & Feuchtinger 1995).

2.4.2 Opacity

The stellar opacity used in the form $\kappa = \kappa(\rho,T)$ is also included in tabulated form. As long as the pulsational driving is located at high optical depths it is appropriate to adopt the Rosseland mean of $\kappa$ which yields the correct momentum transfer in the diffusion limit. However, if e.g. extended atmospheres of pulsating AGB stars involving dust formation and radiation pressure driven winds are investigated, it is important to use other proper means of the opacity (Höfner et al. 1998). A detailed discussion about mean opacities and gray radiative transfer can be found in Mihalas & Mihalas (1984, page 355).

2.4.3 Eddington factor

 For the closure of the two radiation moment equations an equation of state for the radiation field by means of the Eddington factor $f_{\rm edd} = J/K$ is specified. For the case of classical stellar pulsations (e.g. Cepheids, RR Lyrae) where the stellar atmosphere does not influence the pulsational behavior (Feuchtinger & Dorfi 1994), the Eddington approximation ($f_{\rm edd} = 1/3$)which accounts for anisotropies in the radiation field up to the first order, can be used. However, if the atmosphere is taken into account, as e.g. in Feuchtinger & Dorfi (1994) or Höfner & Dorfi (1997), it is important to employ a variable Eddington factor. In that case the radial structure of the Eddington factor is calculated by solving the time-independent grey radiative transfer equation (e.g. Balluch 1988).

2.5 Source function

For the source function S(T) appearing in the energy coupling between matter and radiation we assume local thermodynamical equilibrium and neglect scattering. Consequently S(T) is given by the Stefan-Boltzman law:

$$S(T) = \frac{\sigma}{\pi} \, T^4.$$ (1)

2.6 The model equations

In Table 1 the complete set of model equations for gas dynamics, radiative transfer and turbulent convection in differential form are compiled. With regard to the numerical method of solution using a conservative volume discretization on an adaptive mesh (cf. Sect. 3), these equations have to be reformulated in conservation form and for a time-dependent coordinate system. This is done by integration over a time-dependent volume leading to the integral form of the CRHD equations as given in Table 2. Note that due to numerical reasons we solve an equation for the total internal energy ($\rho e + 4\pi/c\,J + \rho \bar{\omega}$) instead of the gas internal energy equation. The corresponding total internal energy equation (Eq. (4) in Table 1) can be derived by adding the respective energy equations of gas, radiation and convection (Eqs. (4), (5) and (8) in Table 1).

  

Table 2: Equations of convective radiation hydrodynamics in integral form. d/dt denotes the time derivative with respect to a moving coordinate system which is essential for discretization on an adaptive grid. V=V(t) denotes a time-dependent volume and $\partial V$ its surface. Refer to Table 6 for a comprehensive list of symbols and to the text for further details

2.7 Boundary conditions

 For modelling stellar pulsation the CRHD equations are solved for a spherical shell representing the region of the star which is dominated and influenced by the oscillations.

Within the inner boundary located at $R_{\rm int}$ we specify a rigid core which influences the model by its mass $M_{\rm int}$ and luminosity $L_{\rm int}$.Depending on the considered type of star the inner radius $R_{\rm int}$ is constrained by two demands: First it has to be located above the zones where nuclear burning takes place, as we do not consider these effects. Secondly it has to be located deep enough in order to ensure that the eigenfunctions of the pulsations vanish at the inner boundary. Consequently for the mechanical inner boundary we assume the gas velocity to vanish (u = 0). An estimate for $R_{\rm int}$ can be obtained from linear stability analysis. However, the independence of the nonlinear solution on $R_{\rm int}$ has to be checked. For the turbulent convection at the inner boundary we assume the turbulent kinetic energy $\bar{\omega}$ to vanish which constrains the core to be radiative.

The outer boundary condition depends on the problem under consideration. If we are interested in classical stellar pulsations (e.g. Cepheids, RR Lyrae) which are strongly bound objects, oscillating weakly non-adiabatic and exhibiting negligible mass loss rates we use a Lagrangian outer boundary located above the photosphere. For that case the motion of the boundary is determined from the equation of motion. This requires that the medium outside the boundary is specified by an external pressure $P_{\rm ext}$ which we assume to be constant. Note that $P_{\rm ext}$ has to be sufficiently low compared with the pressure at the photosphere (Feuchtinger & Dorfi 1994). Physically this means, that pulsational waves are almost perfectly reflected at the steep pressure decline of the atmosphere which is a good approximation for these stars. Consequently the penetration of waves into the stellar atmosphere are not considered which has no effect on the pulsational behavior (Feuchtinger & Dorfi 1995). Again it has to be checked that the nonlinear solution is independent of $P_{\rm ext}$. Usually we adopt a value of $P_{\rm phot}/P_{\rm ext} = 30$ with $P_{\rm phot}$ denoting the photospheric pressure.

If, on the other hand, pulsations of extended envelopes like AGB stars are to be investigated (which as a matter of fact is a really tough problem), because of the high mass loss rates expected in those stars, a Lagrangian boundary is not suitable. Pulsations are highly nonadiabatic meaning that transmission of mechanical energy into the atmosphere cannot be neglected. Consequently the outer boundary condition has to account both for a considerable expansion of the atmosphere and for the mass outflow generated by the pulsations. For further discussion we refer to Feuchtinger et al. (1993).

In addition to the mechanical outer boundary conditions a condition for the radiation field must be specified. Here we assume the radiative flux to have only an outward component leading to

$$J = \frac{1}{\bar \mu} \, H$$ (2)

with $\bar \mu$ denoting a quantity accounting for the geometry of the radiation field at the outer boundary. If we assume Eddington's approximation (cf. 2.4.3), $\bar \mu$ is given by $\bar \mu = 1/2$while in the case of a variable Eddington factor $\bar \mu$ has to be determined from the radiation transport equation (Balluch 1988).

For the convective outer boundary we assume the gradient of the turbulent kinetic energy $\partial\bar{\omega}/\partial r$ to vanish which in principle allows a convective flux across the external radius. However, the outer boundary usually is located somewhere in the atmosphere where normally no convection occurs.

In Table 3 the complete set of boundary conditions for the case of classical stellar pulsations is summarized.

  

Table 3: Boundary conditions for the case of classical stellar pulsations (see text for details)

  

Table 4: Hydrostatic and local limit of the convective radiation hydrodynamics equations given in Table 1. See text for details

2.8 Static equations and initial conditions

 In order to compute the nonlinear limit cycle evolution of a particular pulsating envelope we start from a hydrostatic configuration referred to as the initial model. Due to the implicit numerical method we use to solve the nonlinear equations (cf. Sect. 3) the structure of the initial model must be an exact solution of the full nonlinear equations. Numerically speaking the initial model has to lie within the convergence radius of the implicit solution scheme. Consequently we use the hydrostatic and local limit of the nonlinear equations given in Table 1 to compute the structure of the initial model. This limit is obtained by omitting all time derivatives, setting the gas velocity identically to zero ($u \equiv 0$), and neglecting turbulent pressure and overshooting. The remaining equations compiled in Table 4 yield a system of ordinary differential equations similar to the equations of stellar structure. Together with the constitutive relations given in Sect. 2.4 this system in principle can be solved by a forward integration.

However, for the case of pulsating envelopes without considering the stellar atmosphere the static equations can be further simplified by the following assumptions:

First we assume Eddington's approximation. Secondly we set J = S which is true as long as layers are convectively stable. In this case the equations given in Table 4 reduce to the equations of stellar structure (e.g. in Kippenhahn & Weigert 1990). If now, during the forward integration, a layer becomes convectively unstable (as determined by the Schwarzschild criterion) we replace the radiative temperature gradient in the radiation diffusion equation by the actual convective gradient, which can be determined from the static turbulent kinetic energy equation (Eq. (6) in Table 4). Analogous to the cubic equation of the mixing length theory a nonlinear equation for the actual temperature gradient can be derived (Wuchterl & Feuchtinger 1998).

The above outlined method represents only an approximative solution of the static system compiled in Table 4. In particular at the transition regions between radiative and convective zones where both radiative and convective fluxes exhibit a nonvanishing divergence, the assumption J = S is not fulfilled. However, in the majority of applications the solution lies well within the accuracy range necessary for the implicit solution method.

Due to practical reasons we start the integration at the photosphere where the initial values for mass, temperature, and luminosity are set to the stellar parameters of the model. Additionally we need the photosphere pressure which at first is set to an estimate. Beginning from the photosphere in a first step we perform an outward integration up to the desired external radius (cf. Sect. 2.7). At this point the outer boundary condition for the radiation field has to be matched and this relation is used to determine the photosphere pressure by iteration. After the atmosphere structure is computed in this way we continue by performing an inward integration from the photosphere down to the desired radius of the inner boundary condition.

In order to guarantee the initial model computed by forward integration to be hydrostatic with respect to the nonlinear equations, in a further step the nonlinear code is used to relax the model on the complete CRHD system. This is done by using the same numerical parameters as for the fully dynamical computation. A successful relaxation is done if the time step of the nonlinear evolution becomes large compared to the thermal time scale of the model, which normaly establishes after 10 - 20 steps.

Even though the dynamical evolution is started with time steps lower than the dynamical time scale of the model, due to numerical dissipation a hydrostatic, but vibrationally unstable model usually remains hydrostatic, if the growth rates of the oscillations are small. In this case the pulsations are initiated by a velocity perturbation applied to the hydrostatic initial model. Depending to the problem under consideration the perturbation is taken stochastic or corresponding to the eigenfunction of the expected pulsation mode (Stobie 1969).