A. Fokin ^1,2,4 - Ph. Mathias ² - E. Chapellier ² - D. Gillet ³ - N. Nardetto ²

1 - Institute of Astronomy of the Russian Academy of Sciences, 48 Pjatnitskaya Str., Moscow 109017, Russia
2 - Observatoire de la Côte d'Azur, Dpt. Gemini, UMR 6203, 06304 Nice Cedex 4, France
3 - Observatoire de Haute Provence, CNRS, 04870 Saint Michel l'Observatoire, France
4 - Isaak Newton Institut Moscow Branch, Russia

Abstract
A hydrodynamical model is generated for the high-amplitude $\beta$ Cephei star BW Vulpeculae, and the spectral line profiles are calculated for different pulsational phases. The pulsational characteristics and line profiles are compared with recent observational data obtained during seven consecutive nights in August 2000. We found a generally good agreement in the basic photometric and spectral parameters. Two strong shock waves appear during one period, and the "stillstand'' is due to the gas dynamics between the passages of these shocks. Note that this good agreement suppose a metallicity Z=0.03, while a metallicity Z=0.02 does not lead to the correct amplitudes and shapes of the curves.

Key words: line: profiles - stars: oscillations - stars: variables: general - stars: individuals: BW Vulpeculae

1 Introduction

Among the $\beta$ Cephei stars, BW Vulpeculae (HD 199140, B2 III) exhibits the most extreme variability of light, radial velocity and line profiles. With a period of 0.201 day, the peak-to-peak amplitude of the radial velocity variation amounts to more than 2K=200 km s^-1, the total range of the light variation is approximately 0.2 mag in V, and, finally, the spectra show well-marked line-doubling. A prominent feature of both the radial velocity curve and the light curve is the presence of a bump. This bump occurs around pulsation phase $\varphi=0.8$ in the light curve, whereas it occurs at pulsation phase $\varphi=1.0$ in the radial velocity variation and is usually called "stillstand''. On each side of the stillstand, the velocity curve shows discontinuities due to line-doubling phenomena.

Different hypotheses have been given to explain such a behaviour. Since the engine of the pulsation involves iron lines, Cox et al. (1992) suggest that the large amplitude of the star could be due to a particular metal enrichment. The peculiar evolutionary status of the star has also been put forward. Indeed, since the period grows at a large rate, 2.37 s century^-1(Horvath et al. 1998), the star is thought to be rapidly expanding on a time scale compatible with the shell-hydrogen burning phase, and thus should be more evolved than "classical'' $\beta$ Cepheids which are thought to be at the end of the core-hydrogen burning phase (Sterken & Jerzykiewicz 1990).

The interpretation of these unusual observational phenomena in a $\beta$ Cephei star is not clear yet. Using a linear mode identification technique performed on both spectral and light variations, Aerts et al. (1998) show that the pulsation mode is radial. In this framework, the actual hydrodynamical picture involves two shock waves per pulsation period, which explains the line doubling in the framework of the Schwarzschild mechanism (1954). A possible scenario (Mathias et al. 1998, hereafter MGFC) is that a wave, originating from the inner layers where the $\kappa$ -mechanism due to iron lines acts, brakes before reaching the photosphere (since even the weak metallic lines are double). This shock imparts an outward motion to the atmospheric layers which then follow a ballistic motion until they fall back. At this stage, the lower layers are more decelerated than the upper ones and induce a velocity gradient that becomes so large that it also brakes into a second shock. In this view, the stillstand in the velocity curve represents a relaxation phase of the atmosphere.

However, this global scenario has several variants. For instance, Young et al. (1981) suggest that a stationary layer is generated during the infalling atmospheric motion by the strong increase in temperature and gas pressure, and that the line doubling is not the result of the propagation of a shock front. Smith & Jeffery (2003) use thermodynamical considerations, especially temperature variations, to explain the Van Hoof effect, since no phase-lag is detected between the different optical lines they considered. In MGFC we were able to measure a phase-lag between the Si III and H $\alpha$ variations which was interpreted in terms of the presence of progressive waves (Mathias & Gillet 1993).

All the ideas mentioned so far were purely based on good observational material, but a model is desirable to clarify what happens in the star. Up to now, the only attempt to model the observations was performed by Moskalik & Buchler (1994). They used a nonlinear pulsation model where the dynamics was governed by a unique outward propagating shock originating at the bottom of the He II ionization zone. In this view, the consecutive strong compression provokes a sudden jump of the Rosseland mean opacity which contributes to the formation of an apparent discontinuity in the observed radial velocities. However, this results in a stillstand which is at a value of about -100 km s^-1 in the rest frame of the star, whereas its observed value is around the stellar $\gamma$ -velocity at -9.2 km s^-1 (MGFC).

The main objective of this paper is to interpret observations of BW Vulpeculae using an auto-coherent pulsation model which has already been successfully used for different classes of radial pulsators, from RR Lyrae (Fokin & Gillet 1997) to RV Tauri (Fokin 2001) and post-AGB (Jeannin et al. 1997). The paper is organized as follows: in Sect. 2, we briefly describe our new high-time-resolution spectra. The pulsation model is applied to these new data in Sect. 3, while Sect. 4 deals with a comparison of line profile variations derived from observations and those computed from the model. Some concluding remarks are given in Sect. 5.

2 Observations

Spectra were obtained at the Observatoire de Haute-Provence with the 1.52 m telescope using the AURELIE spectrograph during 7 consecutive nights, from August 14 to 21, 2000. The spectral resolution was around 25 000 over a 120 Å spectral range centered on the Si III triplet at 4552, 4567 and 4574 Å. This relatively low spectral resolution allowed a very good temporal sampling: with a mean exposure time of 2 min, more than 100 spectra per pulsation period were obtained. The measured signal-to-noise is between 100 and 150. Reductions were performed using the standard IRAF package.

Our observations are presented in phase according to the pulsation period. By convention, the pulsation phase $\varphi=0$ corresponds to maximum luminosity. Since we have no simultaneous photometric observations, dates of maximum luminosity were computed following the most recent ephemeris provided by Horvath et al. (1998). Because this ephemeris gives dates concerning light minima, we added 0.1116 d (0.555 P) to retrieve the usual phase convention (Sterken et al. 1987).

3 Nonlinear model

3.1 Model description

The basic stellar parameters for BW Vul are still uncertain. According to different authors, the mass is between 11 and $14~M_{\odot}$ , luminosity between $\log L/L_{\odot}=4.146$ and 4.431 and $\log T_{\rm eff}$ between 4.33 and 4.386 (Aerts et al. 1998; Lesh & Aizenman 1978; Heynderickx 1992; Moskalik & Buchler 1994). We tried different sets of parameters, and finally chose a model close to the second turn-over point on the $11~M_{\odot}$ evolutionary track of Dziembowski & Pamyatnykh (1993). The parameters of this 150-zone model are: $M=11~M_{\odot}$ , $\log L=4.176$ , $\log T_{\rm eff}=4.362$ , which is close to the BW Vul model published by Moskalik & Buchler (1994).

Our model was calculated with the radiative Lagrangian code by Fokin (1990) which uses variable Eddington factors and a time-dependent transfer equation. The inner boundary was fixed to $T=4.4\times 10^7$ K, corresponding to about 5% of the photometric radius and the envelope contained 83% of the stellar mass. We used the OPAL92 opacity tables, and studied both Z=0.02 and 0.03 metallicity. We started the calculations with an initial velocity profile of 10 km s^-1at the surface.

The $\beta$ Cepheid models are characterized by an extremely slow growth rate, so usually the limit cycle is achieved after some 10⁵ pulsational cycles.

The model with Z=0.02 (the metallicity used by Moskalik & Buchler 1994) reached its fundamental limit cycle with the period of 0.211 days and bolometric and radial velocity amplitudes $\Delta m=0.15$ mag and 2K=40 km s^-1, respectively. The relative radial amplitude at the surface is $\Delta R/R=2.5\%$ . Its pulsation is sinusoidal and synchronic in most of the atmosphere. This model represents a typical $\beta$ Cepheid star, having very small amplitude, but has little in common with BW Vul.

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig1.ps} \end{figure}$	Figure 1: Theoretical bolometric light curve ( upper diagram) and displacement of different mass zones ( lower panel) for a BW Vul model with $M=11~M_{\odot}$ , $\log L/L_{\odot}=4.176$ , $\log T_{\rm eff}=4.362$ and Z=0.03.
Open with DEXTER

The model with Z=0.03, on the contrary, has reached the limit cycle with very large amplitudes, $\Delta m_{\rm bol}=0.7$ mag and 2K=260 km s^-1, with the period P=0.217 days. The relative radial amplitude at the surface is $\Delta R/R=12\%$ . This model is presented in Figs. 1-6. In Fig. 1 one can clearly see the bump of the light curve at phase 0.8, as well as complicated motions in the upper atmosphere with shock waves. According to Barry et al. (1984), the estimated bolometric magnitude should be about 0.75 mag, which is close to our theoretical value of 0.7 mag. Unfortunately, the bolometric light curve cannot be directly compared with the observational curve since the observations give only the photometry in a few selected spectral bands. On the other hand, the theoretical estimation of the visual light curve is not very certain. Notably, it is difficult to estimate $T_{\rm eff}$ with sufficient accuracy (i.e. the temperature at $\tau = 2/3$ in the chosen continuum) in a Lagrangian mesh because of strong variations of $\tau$ and T over only a few mass zones.

We note that, as further test calculations have shown, the amplitudes and the character of the pulsation are not sensitive to small variations of L, $T_{\rm eff}$ and M. For instance, we tried several models located on the evolutionary tracks for M from 10 to 14 $M_{\odot}$ . We varied $T_{\rm eff}$ by 16 $\%$ and L from 10 000 to 18 600 $L_{\odot}$ , and obtained only slight variations in $\Delta m$ (from 0.5 to 0.7) and 2K (from 195 to 270 km s^-1). All these models have the characteristic bump in their light curves and a stillstand in the velocity curves. The period of the models varied as well, so our choice of the BW Vul model was motivated by the closeness of the theoretical to the observed periods.

As we show below, the bump and the stillstand are the results of a passage of two strong shocks formed close to the region of instability ( $T \approx 2.5\times 10^5$ K). The Linear Non Adiabatic (LNA) analysis reveals that there is no low-mode resonance in the BW Vul model up to the third overtone, so the Cepheid-like explanation of the bump is not relevant.

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig2.ps} \end{figure}$	Figure 2: Velocity curves for different mass zones of the BW Vul model Z=0.03. The curves are shifted relative to each other for clarity. The scale is the same for all zones, 30 km s^-1 between two tick marks.
Open with DEXTER

In Fig. 2 we represent the velocity curves for all mass zones in the outer atmosphere. The scale is the same for all zones, and the curves are shifted to each other for clarity. Two main shocks are clearly seen, at about $\varphi=0.8$ and 1.0. To confirm them, we plot the maxima of the compression rate for the most important compression/shock waves versus the mass zone in Fig. 3. We also indicate the mass position of the zone T=250 000 K where the Z-peak $\kappa$ -mechanism acts, the photosphere, and both boundaries of the He ionization zone. Figure 4 shows the maxima of the compression rate vs. radius. Here we must note that between phases 0.7 and 0.85 the first shock is receding, as noted by MGFC. Indeed, the velocity of the falling atmosphere during this time largely exceeds the velocity of the shock, so the first shock is captured by the falling matter until phase 0.85-0.9, when the shock becomes upraising in radius. On the contrary, the second shock is from the beginning upraising in radius. We recall that the shocks are always propagating outwards $\it {in~ mass}$ .

We also remark that the outer boundary of the He ionization zone (at about T=40 000 K) remains strangely flat during the whole phase of "stillstand'' between the two shock waves (Fig. 4).

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig3.ps} \end{figure}$	Figure 3: Shock propagation through the mass grid. Points mark the maxima of compression corresponding to the most important compression/shock waves. The lower solid curve below corresponds to the middle of the Z-peak region. The two upper solid curves limit the He ionization region. The dashed curve indicates the level of the photosphere.
Open with DEXTER

3.2 Shock dynamics

The general picture of the shock dynamics in the BW Vul model is as follows. In the expansion phase soon after $\varphi=0.1$ the atmosphere decelerates and then (after phase 0.45) falls down - with almost constant deceleration, $\approx$ 13 m s^-2, which is about 4 times less than the mean gravity in the model atmosphere. During the contraction phase, the compression of the gas is not homogeneous. At $\varphi = 0.6$ in the region of T=100 000-250 000 K the rate of compression exceeds that in the regions above.

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig5.ps} \end{figure}$

Figure 5: Time evolution of the profiles of luminosity ( upper diagram), gas pressure ( middle) and velocity ( bottom) versus the number of the mass zone between phases 0.69-0.87, corresponding to the generation of shock 1. Phase 0.69 refers to the thick solid curve, and phase 0.87 to the thin three-dotted curve. Mass zone number 60 approximately corresponds to the Z-peak zone having T=250 000 K, while mass zone 150 corresponds to the top of the atmosphere.

Open with DEXTER

At the same time the luminosity from the inner zones starts increasing rapidly, but it is effectively absorbed in the outer region of the Z-peak zone (T=200 000 K). From phase 0.8 (the beginning of the "stillstand'') until 0.95 this absorption is especially strong (see Figs. 5-6). This absorption creates an over-pressure above the Z-peak zone ( $\kappa$ -mechanism). After approximately phase 0.75 the compression wave, caused by this over-pressure, starts propagating outwards and shortly transforms into a shock wave, hereafter called shock 1. This is clearly seen in Fig. 5.

Soon this shock enters the zone of He ionization, where it reaches a velocity amplitude of about 100 km s^-1and a density compression rate of about 13. Due to the perturbation of the temperature and density of the gas, the opacity in the wake of this shock increases by a factor of 2. Consequently, the strong radiative absorption in the region of the Z-peak and in the wake of shock 1 is the cause of the observed bump in the light curve. On the other hand, as we show below, there is no strict stillstand in both the observed and theoretical velocity curves.

Shock 1 then increases in amplitude, up to 140 km s^-1, and reaches a compression rate of about 100. It rapidly passes through the atmosphere and escapes. During the photometric bump (until the escape of shock 1) the absorption in the Z-peak zone and above it still continues, also due to the increasing opacity in the wake of shock 1.

After the escape of shock 1 the outer envelope starts expanding, while the inner shells are still in compression. The expansion of the outer atmosphere is slow. The absorption in the Z-peak region is still about 11 000 $L_{\odot}$ , but the luminosity from the inner region increases, so the total stellar luminosity starts increasing after the short bump. However, the accumulation of thermal energy due to absorption in this zone continues.

Near phase 0.95 the compression starts expanding outwards from the inner zones (see Fig. 6). Soon after, a new shock is formed in the helium ionization zone ("shock 2'').

The outer atmospheric layers are now falling, accelerated by the gravity. The absorption in the Z-peak region stops, but the luminosity from the inner region also becomes less, so the surface luminosity soon starts decreasing.

Just before phase 1.1 the inner (postshock) envelope expands with a velocity of 110 km s^-1, while the outmost (preshock) atmospheric layers fall with a velocity of 50 km s^-1. At $\varphi=1.1$ the second shock arrives at the surface, and the whole envelope starts expanding again (main expansion).

An important result from the above analysis is that both shock waves have their origin in the region lying well below the photosphere (the Z-peak zone for shock 1, and the Helium ionization zone for shock 2), and are seemingly due to the $\kappa$ -mechanism in the Z-peak zone.

4 Line profiles: theory vs. observations

After the model was generated, we calculated a series of snap-shots of the atmospheric structure (about 50 per pulsational period) to study the theoretical line profiles. The line transfer problem for each model atmosphere was solved with the code of Fokin (1991) under the LTE assumption. For all phases we assumed the same microturbulent velocity of 1 km s^-1, and we considered a projected rotation velocity of $v~\sin i=24$ km s^-1 (Stankov et al. 2003). The line profiles have been convolved with the relevant instrumental profile.

In Figs. 7 and 8 we present the comparison of the theoretical and observed profiles of the Si III 4553 Å line for almost a full pulsational cycle.

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig7.ps} \end{figure}$	Figure 7: The sequences of calculated (solid curves) and observed (dots) profiles of the Si III 4553 Å line for phases 0.239-0.664.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.1cm,clip]{0418-fig9.ps} \end{figure}$	Figure 9: Theoretical (solid) and observational (points) Full Widths at Half Magnitude of the Si III 4553 Å line. The comparison is presented for all seven consecutive nights, with the dates indicated in each diagram.
Open with DEXTER

In Figs. 9-11 we present the detailed comparison of different features of the predicted profiles versus the observational ones obtained on 7 consecutive nights in August 2000. Note that the theoretical curve is the same in each diagram, and is compared with the observed curves for different nights. These three diagrams show, respectively, the FWHM (Fig. 9), the residual flux (RF: Fig. 10) in the minima, and the radial velocities measured at the minima of the principal absorption components (Fig. 11).

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig10.ps} \end{figure}$	Figure 10: Same as in Fig. 9 for the residual flux of the principal component. We note that the theory-observation discrepancy in flux is normally less that 0.1, or less than 15%.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{0418-fig11.ps} \end{figure}$	Figure 11: Same as in Fig. 9 for the radial velocities measured from the minima of the principal absorption component.
Open with DEXTER

Note that the observational curves are noticeably variable from cycle to cycle. If we fix some phases, we shall see that the discrepancy between the observed and predicted curves also varies significantly from night to night. This behaviour cannot be explained by our model, which is strictly periodic.

The comparison of the FWHM curves (Fig. 9) is good, and for most phases and nights even excellent. This is due to the fact that the FWHM is less sensitive to opacity or temperature variations.

In contrast with the FWHM curves, the RF diagrams (Fig. 10) show greater differences between theory and observations. A good agreement is obtained at phases 0.07-0.4 and 0.75, while between the phases of the shock passages the difference reaches more than 0.1. We suggest that the temperature of the gas after the first shock passage is not correct in the model. This can probably be improved after correction of the radiative cooling rate in the shock wake. Nevertheless, we are optimistic about these diagrams since the relative discrepancy of the flux is less than 15%.

Finally, the radial velocity curves (Fig. 11) show very good agreement with the observations for almost all the nights. We note that the "stillstand'' is rather an idealization, especially clearly seen in the curve of August 16: during this period the velocity varies significantly, a fact confirmed by our model.

5 Discussion and conclusion

Our nonlinear model represents reasonably well the main observed features of BW Vulpeculae, i.e. period, amplitudes, stillstand and bump, light and velocity curves. We also confirm that two shocks are generated consecutively in the stellar envelope - one at each phase of the observed velocity discontinuities (and not one shock as previously stated by Moskalik & Buchler 1994). These shocks are extremely rapid, so the corresponding line-doublings are very short (less than 0.02 P). The physical origin of these shocks is not very clear yet, but we suggest that the main mechanism is a strong radiative absorption in the zone of the "Z-peak''. Also, the characteristic asynchronous motions of the upper and lower envelope regions can contribute to the shock generation. We stress that both shocks are generated below the photosphere.

It is noteworthy that the first shock revealed by the model is at first receding in radius, while the second one is always rising, in agreement with what was described in MGFC.

We found that a metallicity Z=0.02 is too low, and that the observed amplitudes can be reproduced only with Z=0.03. In the limit cycle regime the maxima of the kinetic energy $\log E_k$ for these two metallicities are 42.47 and 43.26, respectively. The linear analysis shows that the only driving zone in the models is that of the "Z-peak'' of opacity at $\log T=5.2{-}5.4$ . Just above and below there are two regions of positive dissipation, with the deep one being more effective. LNA calculations show that the model with Z=0.03 is unstable in the F-mode, while the Z=0.02 model is only marginally unstable. We found that the lower damping zone becomes more effective as Z decreases from 0.03 to 0.02. Apart from being a direct effect of the metal abundance, it seems that this is also due to the fact that the work integrand in the lower dissipation zone increases by an equal amount. The higher relative amplitude in this zone can provoke higher gradients and, consequently, the growth of positive dissipation, including the nonlinear regime. Unfortunately, linear analysis alone cannot explain definitively the low energetics of the Z=0.02 model, and can give only a rough idea.

We suggest that other members of the $\beta$ Cepheid group, which have much smaller amplitudes than BW Vul, must have lower metallicity, as is shown in our second model, which is identical except for Z=0.02.

Hydrodynamic models for $\beta$ Cephei variables

I. BW Vulpeculae revisited