A&A 396, 937-948 (2002)
DOI: 10.1051/0004-6361:20021534

On the role of duplicity in the Be phenomenon

I. General considerations and the first attempt at a 3-D gas-dynamical modelling of gas outflow from hot and rapidly rotating OB stars in binaries

P. Harmanec1,2 - D. V. Bisikalo3 - A. A. Boyarchuk3 - O. A. Kuznetsov4


1 - Astronomical Institute of the Charles University, V Holesovickách 2, 180 00 Praha 8 - Troja, Bohemia, Czech Republic
2 - Astronomical Institute of the Academy of Sciences of the Czech Republic, 251 65 Ondrejov, Bohemia, Czech Republic
3 - Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya 48, Moscow 109017, Russia
4 - Keldysh's Institute of Applied Mathematics, Miusskaya Square 4, Moscow 125047, Russia

Received 24 July 2002 / Accepted 14 October 2002

Abstract
This paper begins a new series of studies devoted to a critical re-examination of the role of duplicity for the Be phenomenon and for the variability patterns observed for many Be stars.
Based on both dynamical and energy considerations and a numerical gas-dynamical modelling, a new hypothesis of the formation of Be envelopes in binaries, via an outflow from a rapidly rotating B star in a detached binary, is outlined. It is shown that such an outflow is facilitated by the presence of a companion to the B star and leads to the formation of an envelope but not to any significant mass exchange between the binary components.

Key words: stars: emission-line, Be - stars: binaries: close - stars: binaries: spectroscopic - stars: fundamental parameters


1 A brief critical evaluation of the existing models of Be stars

1.1 The Be phenomenon

More that one tenth of all B stars exhibit time-variable emissions in their Balmer (and some other) line profiles. These stars, called Be stars, represent perhaps the most variable objects among massive stars. Their continuum and line spectra vary on several time scales, ranging from minutes to decades or more. Actually, the longest time scale of their variability is not known and it is conceivable that every B star may become a Be star for an interval during its evolution. The phenomenon also partly extends to hotter Oe stars and cooler Ae stars.

In spite of all the effort of several generations of stellar astronomers, neither the nature of the Be phenomenon nor the physical causes of variations are well understood at present. There is general agreement that the emission lines are formed in extended gaseous envelopes around these stars. Such envelopes have dimensions of one or even two orders of magnitude larger than the stars themselves and re-radiate the stellar radiation in all directions. Evidence has been accumulated to confirm Struve's (1931) original suggestion that the Be envelopes are rotationally flattened disks. An ultimate proof of this statement came from spectro-interferometric and spectro-polarimetric observations of several Be stars which allowed a partial spatial resolution of their envelopes - see Quirrenbach et al. (1993,1997).

However - as far as the understanding the very complex Be phenomenon is concerned - the current situation with many partial models, explaining often only one aspect of it, is not very satisfactory.

1.2 Attempts to explain the Be phenomenon

Concerning the very origin of the Be envelopes, the number of hypotheses has been growing steadily and there is no general agreement on a single one. These hypotheses are as follows:
1.
Rotational hypothesis by Struve (1931) explains the formation of Be envelopes by rotational instability of underlying, rapidly rotating stars. Its main argument is the observed correlation between the width of the observed Balmer emission lines and v sin iof the respective stars (which has been confirmed by several later studies - cf., e.g., Hanuschik et al. 1988). Struve himself was aware of the fact that the rotational model does not explain variability, namely the long-term V/R variations, observed for a number Be stars at certain epochs. He ingeniously suggested a slow spatial revolution of temporarily elongated envelopes as an explanation. Critics of the rotational model argue that one observes a number of similarly rapid rotators, called Bn stars, which were never observed to have any emission lines (e.g. Baade 1992). Many of the investigators are also convinced that the Be stars are rotating below their break-up speeds at the equator (cf., e.g., Porter 1996). Boyarchuk (1958) pointed out that the rapid rotation itself cannot account for the Be phenomenon and all later quantitative models led to the same conclusion (see, e.g., Limber & Marlborough 1968). Note also that the rotational hypothesis itself does not offer a clue to the observed variations on several time scales.
2.
Outflow models Several attempts were made to explain the origin of Be envelopes as a spheroidal outlow of gaseous material from the underlying star: the stellar-wind model by Gerasimovic (1934,1935), the variable mass flux model by Doazan & Thomas (1987) and Doazan (1987), the "bi-stable'' or "axi-symmetric'' radiation-driven wind model by Lamers & Pauldrach (1991), Araújo et al. (1994) and Stee & Araújo (1994) and the wind-compressed disk model by Bjorkman & Cassinelli (1993) (see also the review by Bjorkman 2000). The latest one was considered as a very promising idea several years ago but the radiation-hydrodynamics simulations by Owocki et al. (1996) showed that the small non-radial components in line forces inhibit formation of an equatorial disk by the wind-compressed model. Also, this model could have problems providing strong enough stellar winds in Be stars of later spectral subclasses. More generally - outflow models do not offer an explanation for the variability of Be stars.
3.
Non-radial pulsations of Be stars were first suggested to explain the observed rapid light and line-profile variations (Baade 1979,1982; Bolton 1982) but the idea has been developed in a complex attempt to explain the Be phenomenon by Baade and his students. In particular, Rivinius et al. (1998a) argued that constructive interference of several pulsational modes can lead to release of a new Be envelope and claimed agreement with the observed emission-line episodes for $\mu$ Cen. They did not present any energy balance considerations of whether such a mechanism could really work, however. It is certainly true that the series of line profiles obtained for a few Be stars are remarkably similar to theoretical line profiles, based on non-radial pulsation modelling. It is conceivable, however, that the description of line-profile variations in terms of spherical harmonic functions would also work to model line-profile variations due to co-rotating structures, located slightly above the stellar surface, which has been advocated, for instance, by Harmanec (1989), Smith et al. (1998) or Balona (2000). Note also that the attempts to derive the macroscopic stellar properties from the modelling of line-profile changes, interpreted as non-radial pulsations, led to too small values of stellar radii, contradicting the estimates based on Hipparcos parallaxes (compare the Maintz et al. 2002 radius of 3.4 $R_\odot $ of the B2e star $\mu$ Cen from the line-profile modelling to Harmanec's 2000 estimate of 5.2-5.4 $R_\odot $ from the V-band brightness outside emission-line episodes and the Hipparcos parallax). Aerts (2000) also warned that the interpretation of Be stars in terms of non-radial pulsations needs further careful tests. Harmanec (2001) pointed out that the formal description of radial-velocity variations of $\mu$ Cen by six non-radial modes did not lead to the gradual decrease of the rms error of the fit, as would be expected for real multiperiodicity. It is fair to say, however, that until such a modelling as was carried out for the non-radial pulsation model will also be carried out for alternative models, the model of non-radial pulsations represents the best available description of the observed rapid line-profile changes of Be stars.
4.
The helmet-type magnetic field model was presented by Underhill (1983,1987) which assumes the presence of organized magnetic fields in Be stars and interprets their envelopes as helmet-type co-rotating structures. This model has been largely ignored by the community (perhaps due to the fact that there was no observational technique allowing the detection of the putative magnetic fields). Given the observational evidence in favour of co-rotating structures slightly above the stellar photospheres (e.g. Harmanec 1989,1999; Harmanec & Tarasov 1990; Smith et al. 1998 or Balona 2000), and the improving sensitivity to detect even weaker magnetic fields in hot stars, this hypothesis probably deserves further critical examination.
5.
Binary models The last class of models are models that assume that the Be phenomenon is somehow related to the duplicity of the Be stars. They are discussed in detail in the next section.

2 The binary models and the role of duplicity for the Be phenomenon

Kríz & Harmanec (1975) and Harmanec & Kríz(1976) formulated a general hypothesis of the binary nature of Be stars, explaining the Be star envelope as an accretion disk created from gas flowing to the Be star via a Roche-lobe overflow from its unrecognized binary companion. They pointed out that with increasing orbital period one observes either a typical Algol interacting binary or a Be binary or a symbiotic star. Since a given B star occupies less and less space when placed in a binary system with longer and longer orbital period, the space available for the formation of an extended accretion disk also gets larger with increasing orbital period. In systems with shorter orbital periods, the Roche-lobe filling and less massive secondaries are usually much less luminous than their B-type primaries. Only for very long orbital periods do the absolute dimensions of the Roche lobe around the secondary become so large that even a very cool secondary has an optical luminosity comparable to, or even larger than the B star and one observes a symbiotic binary. Kríz & Harmanec (1975) were also able to explain rapid rotation of Be stars as a consequence of tangential accretion from the disk which brings some extra angular momentum. Furthermore, they explained at least some of the observed types of variations of Be stars, for instance the phase-locked V/R, RV, line-width and luminosity changes or long-term E/C and V/R variations and offered some ideas as to how rapid changes could be related to co-rotating structures in the accretion disks, caused by resonances.

While their hypothesis certainly represents a serious attempt to address the Be phenomenon in its complexity and very probably explains the nature of some of the actually observed Be binaries (for example AX Mon, RX Cas, SX Cas, KX And, V360 Lac, $\beta$ Lyr etc. - see Harmanec 2001 for a catalogue of known emission-line binaries), it is now well proven that it cannot be accepted as a universal explanation for the origin of Be envelopes. Already Plavec (1976) pointed out that if all Be stars were binaries with Roche-lobe filling secondaries, one should observe more eclipsing binaries among them than what is actually observed. While Harmanec (1987) slightly weakened this objection, there is a stronger one: detailed studies of several known Be binaries ($\varphi~$Per = HD 10516: Poeckert 1981, Gies et al. 1998; V839 Her = 4 Her: Koubský et al. 1997, for instance) clearly demonstrated that the secondaries in those binaries are not Roche-lobe filling objects but very small stars. The same is also true of binaries composed of a Be star and a compact, X-ray companion. Harmanec (1985) came up with the provocative suggestion that even in massive X-ray binaries the mass is flowing from the X-ray star towards the Be primary and presented some observational facts to support such a view. Recalling an earlier suggestion by Kríz (1982), he also argued that the contraction of the originally mass-losing star had to lead to its rotational instability near the equator, leading to another phase of mass transfer from this star to its (now more massive) counterpart. He called it a case PB of mass transfer and mentioned $\varphi$ Per as a system being possibly in such a mass transfer stage. It is obvious, however, that unless somebody can show how to excite X-ray emission from compact stars without allowing them to accrete mass from their optical companions, Harmanec's (1985) idea is not tenable.

Obviously unaware of Kríz's (1982) and Harmanec's (1985) studies, Pols et al. (1991) also investigated the possibility of formation of Be stars as products of case B mass exchange in binaries. Their approach was different, however. They accepted the idea of Kríz & Harmanec (1975) that some Be stars are case B mass-exchanging binaries but argued that the majority of Be stars are remnants of case B mass exchange in intermediate-mass close binaries after the termination of mass transfer. In other words, they postulated that the Be phenomenon occurs due to some still unknown physical mechanism which is only operational in rapidly rotating stars. The role of the mass exchange in their hypothesis is to rejuvenate and spin-up the original secondaries in binaries. They argued that Be stars in mass-exchanging binaries represent only a small fraction of all Be stars. Estimating the lifetimes of different evolutionary stages, they concluded that more than 80% of post mass-transfer Be stars should have a helium-star companion and that there should be 10 times more Be stars with a white-dwarf companion than those with a neutron-star secondary (observable as an X-ray source). They predicted that many new helium-star and white-dwarf companions should be detectable in the XUV spectral region.

The role of duplicity was critically examined by Baade (1992). Using high-S/N IR spectra near 880 nm, he carried out a search for late-type companions of 35 southern Be stars with a completely negative result. He also expressed some doubts about the existence of many binaries with hot compact companions and his conclusion was that the cause of the Be phenomenon cannot be related to their duplicity.

For the following reasons we believe, however, that the role of duplicity was not still investigated well enough and that Baade's view cannot be accepted as the final word:

1.
Perhaps most importantly, the number of known binaries among emission-line stars has been growing steadily - see, e.g., Gies (2000) or Harmanec (2001). As pointed out by Harmanec (2001), duplicity of Be stars can actually serve two different roles: (i) to explain the formation of Be envelopes via some kind of binary interaction, and (ii) to explain some of the variability patterns of Be stars.
2.
Studies by Slettebak (1987) or Harmanec (2000) demonstrated that Be stars are observed among stars of clearly different evolutionary ages. This may indicate that the cause of the Be phenomenon has to be sought in some external mechanism, not primarily in a physical mechanism related to the stars themselves.
3.
Such a view can be also supported by another similar argument. B stars span a huge range of stellar masses and it is well known from the theory of stellar evolution that the time scales of all processes depend strongly on the stellar mass. However, as pointed by Horn et al. (1982) the time scale of the formation of a new Be envelope was found to be very similar for three Be stars of spectral classes B0, B6 and B8. This again seems to indicate an external mechanism for the formation of such envelopes.
It is obvious that a new generation of powerful optical interferometers will soon be able to resolve many close binaries, so far detectable only by spectroscopy. This should allow new stringent tests for various binary scenarios of the Be phenomenon.

3 Physical model

Let us consider a single, rapidly rotating star. When the velocity of rotation on equator $V_{\rm rot}$ will reach the Keplerian (also called break-up) velocity $V_{\rm br}=\sqrt{GM/R}$, the centrifugal and gravitational attractive forces will compensate each other:

\begin{displaymath}F_{\rm cfg}\equiv\frac{V_{\rm rot}^2}{R}=\frac{GM}{R^2}\equiv F_{\rm grav}. \end{displaymath}

In such situations, the presence of a pressure gradient (not counterbalanced by other forces) permits the matter to outflow from the equatorial belt (called the Roche limit). A particle with specific kinetic energy $\simeq $ $V_{\rm br}^2/2$ rotates along closed trajectories and forms an envelope. To escape to infinity the particle would need to aquire additional energy and reach the parabolic velocity $V=\sqrt{2}V_{\rm br}$.


  \begin{figure} \par\resizebox{\hsize}{!}{\includegraphics{m2949f1.eps}}\end{figure} Figure 1: Relative change in the position of the inner Lagrangian point for the systems with non-synchronous rotation of components $x_{L_1^{\rm rot}}/x_{L_1}$ as a function of the degree of non-synchronous rotation $f=\Omega _\star /\Omega $ for two values of the mass ratio: q=0.1 and q=1. The inserted panel shows the dependence $x_{L_1^{\rm rot}}/x_{L_1}$ vs. q for large values of f.
Open with DEXTER

The situation changes dramatically if the same star is a component of a binary system. Let us consider a binary with a spin-orbit synchronization (the velocities of angular rotation of both components are equal to the velocity of angular revolution of the system) and let us use a Cartesian coordinate system that rotates with an angular velocity $\Omega$ and has its origin in the centre of star 1. The X axis is directed toward star 2, Z axis is parallel to the vector of orbital revolution, and Y axis is so oriented to define a right-hand coordinate system. There are several forces acting on a test particle located between the binary components: gravitational attraction of the two binary components, the pressure gradient, and two forces related to the co-rotating frame used: the centrifugal and Coriolis force. The law of motion of such a test particle was first investigated by Roche (1848,1851) in a ballistic approach (i.e. ignoring the pressure force) as a solution of a restricted three-body problem (assuming the mass to be concentrated into two point masses); for a different formulation, see also Hill (1905). The force field (without Coriolis force and the pressure gradient) in cases of a spin-orbit synchronism can be described by the standard Roche potential $\Phi$:


 \begin{displaymath}\begin{array}{l} \Phi=-{\displaystyle\vphantom{x^2_2}GM_1\ove... ...style\vphantom{x^2_2}M_1+M_2}\right)^2+y^2\right)~, \end{array}\end{displaymath} (1)

where M1, M2 are the two point masses, A is binary separation, $\Omega$ is the angular velocity of orbital motion, and x, y, z are Cartesian coordinates in the adopted frame[*].


  \begin{figure} \par\includegraphics[width=6.8cm]{m2949f2a.eps} \includegraphics[width=6.8cm]{m2949f2b.eps} \includegraphics[width=6.8cm]{m2949f2c.eps}\end{figure} Figure 2: Roche equipotentials (dashed lines) and equipotentials of "asynchronous'' potential $\Psi $ (solid lines) for q=0.1 and different values of f. An arrow shows the linear velocity in $L_1^{\rm rot}$ point which amounts to $0.58~A\Omega$ for the upper panel, $1.88~A\Omega$ for the middle one, and $4.45~A\Omega$ for the bottom panel (the scale of the absolute value of velocity is different for each panel).
Open with DEXTER

The presence of additional forces (absent in the case of a single star) results in violation of equilibrium in the inner Lagrangian point L1. In particular, the pressure gradient cannot be counterbalanced there by the gradient of the Roche potential. Hence, as soon as the star expands and fills the cricital Roche lobe, matter begins to flow towards the binary companion in the vicinity of L1 point but not from the entire equatorial zone. The position of inner Lagrangian point L1 can be derived from equation $\nabla\Phi=0$, which can be rewritten after Kopal (1959):

 \begin{displaymath}\left(\frac{x_{L_1}}{A}\right)^{-2}-\frac{x_{L_1}}{A} =q\left... ...L_1}}{A}\right)^{-2} -\left(1-\frac{x_{L_1}}{A}\right)\right), \end{displaymath} (2)

where q=M2/M1 denotes the mass ratio.

In case of asynchronous rotation of star 1,[*] we should also include centropedal acceleration and Coriolis force into the equilibrium conditions. The presence of these two terms is related to the motion of the stellar matter in the adopted co-rotating frame. In such a case, the force field at the stellar surface is given by asynchronous Roche potential (see, e.g., Plavec 1958; Kruszewski 1963; or Limber 1963):

 \begin{displaymath}\Psi=\Phi-\hbox{$~^{1}\!/_{2}$ }(\Omega_\star^2-\Omega^2) (x^2+y^2) \end{displaymath} (3)

where $\Omega_\star$ is the angular velocity of rotation of the star in question. Similar to the synchronous case it is possible to introduce a concept of the inner Lagrangian point for asynchronous rotation, $L_1^{\rm rot}$, i.e. a point where the pressure gradient ceases to be counterbalanced by other forces and where the matter begins to flow towards the companion when reaching this point. The position of this point is given by the condition $\nabla\Psi=0$ which leads to the following equation (see, e.g., Pratt & Strittmatter 1976)

 \begin{displaymath}\left(\frac{x_{L_1^{\rm rot}}}{A}\right)^{-2}-f^2\frac{x_{L_1... ...t)^{-2} -\left(1-f^2\frac{x_{L_1^{\rm rot}}}{A}\right)\right), \end{displaymath} (4)

which is similar to Eq. (2) valid for synchronous rotation. In this equation $f=\Omega _\star /\Omega $ and the sign of f does not affect the solution, i.e., the sense of the stellar rotation in the laboratory coordinate system does not affect the position of $L_1^{\rm rot}$. Therefore - without loss of generality - we consider only the cases of $f\ge0$, i.e. cases when the directions of stellar rotation and binary revolution are the same. The ratio $x_{L_1^{\rm rot}}/x_{L_1}$ as a function of q and f is shown in Fig. 1. It is obvious that for stellar rotation rates slower than the orbital revolution (f<1), the "non-synchronous'' Roche lobe is larger than the standard Roche lobe, achieving maximum for f=0. When the star rotates faster than the binary revolves, the "non-synchronous'' Roche lobe becomes smaller than the standard one. Note that formally $x_{L_1^{\rm rot}}/x_{L_1}\to0$ as $f\to\infty$. In real binaries, however, the position of $x_{L_1^{\rm rot}}$ is actually limited by the break-up rotation velocity of the star in question, so that $x_{L_1^{\rm rot}}/A$ cannot be smaller than $(\Omega_{\rm br}/\Omega)^{-2/3}(q+1)^{-1/3}$.


  \begin{figure} \par\resizebox{\hsize}{!}{\includegraphics{m2949f3.eps}}\end{figure} Figure 3: The ratio of linear velocity in $L_1^{\rm rot}$ point (in observer's coordinate system) to critical velocity $\eta=\Omega_\star\cdot x_{L_1^{\rm rot}}/\protect\sqrt{GM_1/x_{L_1^{\rm rot}}}\cdot100\%$ vs. q for different values of asynchronicity parameter $f=\Omega _\star /\Omega $.
Open with DEXTER


  \begin{figure} \par\resizebox{10cm}{!}{\includegraphics{m2949f4.eps}}\end{figure} Figure 4: Extra energy $\Delta E$ (in units of $A^2\Omega ^2$) needed for a particle leaving $L_1^{\rm rot}$ to reach the L1 point, plotted vs. asynchronicity parameter $f=\Omega _\star /\Omega $ for values of binary mass ratio of q=0.1 and q=1. Dash-pointed line shows value $\Delta E=0$. The inserted panel shows the dependence of $\Delta E$ vs. f for large values of f.
Open with DEXTER

This is illustrated in Fig. 2 which shows the equatorial plane of a binary with q=1. The isoline of the Roche potential $\Phi$ passing through L1 (i.e. the XY projection of the classical Roche limit) is shown by a dashed line, while the isoline of potential $\Psi $, passing through $L_1^{\rm rot}$, is drawn by a solid line. The "asynchronous Roche lobes'' are shown for three values of "asynchronicity'' parameter: f=2, f=10, and f=100 (see panels a), b), and c) of Fig. 2, respectively). The vectors of linear velocity of $L_1^{\rm rot}$ point ( ${\vec V}_{\rm rot}=({\vec\Omega}_\star-{\vec\Omega})\times{\vec r}_{L_1^{\rm rot}}$) in the adopted coordinate system are also shown.


 \begin{figure} \par\includegraphics[width=8.8cm]{m2949f5a.eps} \includegraphics[width=8.8cm]{r2949f5a.eps} \end{figure} Figure 5a: Top panel: bird's-eye view of density isosurfaces on level $\rho =5\times 10^{-5}\cdot \rho (L_1^{\rm rot})$ for the moment of time $t=0.514P_{\rm orb}$. Middle panel: the slice of formed envelope by XY plane (density distribution and velocity vectors in equatorial plane). Vector in top right corner corresponds to velocity 500 km s-1. Bottom panel: the slice of the envelope by XZ (density distribution in frontal plane). Coordinates in all three panels are expressed in $R_\odot $.


 \begin{figure} \par\includegraphics[width=8.8cm]{m2949f5b.eps} \includegraphics[width=8.8cm]{r2949f5b.eps} \end{figure} Figure 5b: The same as Fig. 5a but for the time $t=0.545P_{\rm orb}$.


 \begin{figure} \par\includegraphics[width=8.8cm]{m2949f5c.eps} \includegraphics[width=8.8cm]{r2949f5c.eps} \end{figure} Figure 5c: The same as Fig. 5a but for the time $t=0.586P_{\rm orb}$.


 \begin{figure}\par\includegraphics[width=8.8cm]{m2949f5d.eps} \includegraphics[width=8.8cm]{r2949f5d.eps} \end{figure} Figure 5d: The same as Fig. 5a but for the time $t=0.608P_{\rm orb}$.


 \begin{figure}\par\includegraphics[width=8.8cm]{m2949f5e.eps} \includegraphics[width=8.8cm]{r2949f5e.eps} \end{figure} Figure 5e: The same as Fig. 5a but for the time $t=0.639P_{\rm orb}$.


 \begin{figure}\par\includegraphics[width=8.8cm]{m2949f5f.eps} \includegraphics[width=8.8cm]{r2949f5f.eps} \end{figure} Figure 5f: The same as Fig. 5a but for the time $t=0.660P_{\rm orb}$.

As discussed above, the matter can outflow from the stellar surface as it reaches $L_1^{\rm rot}$ point. This fact changes the limiting value of the break-up velocity when the outflow begins. In Fig. 3, the values of the linear velocity at $L_1^{\rm rot}$ point are plotted as a function of binary mass ratio q for different values of asynchronicity parameter f. All velocities are expressed in the units of the critical velocity $V_{\rm br}$, derived for a single star of the same properties. The results presented in Fig. 3 show that in a number of cases even a small additional increase in the rotational velocity can lead to an outflow of matter in the vicinity of $L_1^{\rm rot}$ point. However, for values typical for Be stars, say $q\sim0.1$ and $f\sim100$, the outflow occurs for rotational velocities only slightly smaller than the break-up velocity of the respective star. It is important to realize, however, that when  $V_{\rm rot}$ gets close to $V_{\rm br}$ for Be stars, which are members of binary systems, the outflow of matter occurs only in the vicinity of $L_1^{\rm rot}$, not from the the whole equatorial belt of the star (the Roche limit for single stars).

The presence of a companion to the Be star results in another notable change in the mechanism of the outflow from the surface of a rotating star. For a single star, a particle can escape to infinity if it has the parabolic velocity. For a binary star, the particle leaving via $L_1^{\rm rot}$ point can escape from the system if its energy is large enough to reach the L1 point (this energy is smaller than that needed to reach the escape velocity from a single star) since after it reaches L1, it can be captured by the gravitational field of the companion. Particles with energies insufficient to reach the L1 point will move along closed trajectories around the star from which they escaped. The extra energy $\Delta E$, needed to get a particle with velocity $V_{\rm rot}=(\Omega_\star-\Omega)\cdot x_{L_1^{\rm rot}}$ to the vicinity of L1 point, is plotted in Fig. 4 as a function of f for two mass ratios, q=1 and q=0.1, and is expressed in the units of the characteristic energy of the system $E_{\rm sys}=A^2\Omega^2$. It is obvious that the sum of potential energy of a particle in  $L_1^{\rm rot}$ point plus its kinetic energy given by the stellar rotation is much smaller than the potential energy in L1 point

\begin{displaymath}\Phi(L_1^{\rm rot})+\hbox{$~^{1}\!/_{2}$ }(\Omega_\star-\Omega)^2\cdot x_{L_1^{\rm rot}}^2< \Phi(L_1). \end{displaymath}

The value of extra energy needed to reach L1 is a few orders of magnitude larger than the characteristic energy of the system $E_{\rm sys}$. Considering that the effective temperatures of B stars range roughly from 10 000 to 30 000 K, it is clear that the thermal energy cannot change the overall energy balance significantly. One is, therefore, led to the conclusion that the outlow from rapidly rotating B stars in binaries via the $L_1^{\rm rot}$ point should lead to the formation of roughly Keplerian equatorial disks around such stars but not to a significant mass transfer towards their companions.

It is necessary to point out, however, that the above analysis of the energy balance is not exhaustive. For instance, the numerical investigations carried out by Narita et al. (1994) show that if the rotational velocity is close to $V_{\rm br}$ and viscosity is considered, a disk in a binary system may evolve from an accretion disk to an outflowing one. At the same time, their numerical simulations showed that the transfer of the angular momentum is far more pronounced than the mass loss and that the mass loss by the disk is significant only on the evolutionary time scale (i.e. $\sim$106 years). Unfortunately, it is impossible to carry out 3-D gas-dynamical simulations for such long time intervals with present-day computers. One only has to expect that a viscous smearing will not significantly influence the solution obtained only over a time interval comparable to the orbital period of the binary (say, less than a year).

It summary, the outlow of matter via the $L_1^{\rm rot}$ point seems to represent the most probable scenario of the formation of the Be envelope for a rapidly rotating B star which is a member of a binary system.

4 The numerical model

4.1 Parameters of binary

For a numerical investigation, we have chosen a binary with parameters that are typical for binary Be stars. Since Be stars are most abundant around the spectral class B2, we have chosen parameters corresponding to a normal main-sequence B2 star according to the empirical calibration by Harmanec (1988) and adopted a mass ratio of 0.1 which is also typical for known Be binaries. In particular, we used the mass of Be star (star 1) $M_1=8.6M_\odot$, mass of star 2 $M_2=0.86M_\odot$, binary mass ratio q=M2/M1=0.1, binary separation $A=121R_\odot$, and orbital period $P=50^{\mbox{d}}$. Inner Lagrangian point is then located at the distance of $x_{L_1}=87R_\odot$ from the centre of star 1. The equatorial radius of the Be star $R_1=4.6R_\odot$ was chosen. In accordance with the above considerations we assumed the Be component to be large enough to reach the `asynchronous' Roche lobe, i.e. that $x_{L_1^{\rm rot}}=R_1$. Using the adopted values of $x_{L_1^{\rm rot}}/A$ and q, one gets the asynchronicity parameter of Be star $f=\Omega _\star /\Omega $ via Eq. (4). In this case, f is equal to 128.6 (implying that the period of rotation of the Be star is 0 $^{\rm d}\!\!.$39). The value of linear velocity at $L_1^{\rm rot}$ point in the adopted coordinate system is equal to $V_{\rm rot}=(\Omega_\star-\Omega)\cdot x_{L_1^{\rm rot}}=594$ km s-1. Note that the value of the critical velocity of a single star with the same mass and radius is equal to $V_{\rm br}=\sqrt{GM_1/R_1}=599$ km s-1. This means that the value of velocity in $L_1^{\rm rot}$ point equals 99.2% of the critical velocity.

As pointed out above, the matter located in $L_1^{\rm rot}$ can reach the L1 point only when it gets some extra energy $\Delta E$. Figure 4 shows that for the adopted values of q and f it would be necessary to add energy $\Delta E=10.5\cdot A^2\Omega^2$ to reach L1. For the binary in question, the thermal energy needed to provide such an energy excess would require temperatures of at least 1.86 millions of K. In the numerical model, we actually adopt the value of effective temperature corresponding to a normal star of a similar mass, 22 900 K. One then gets $\Delta E=0.013\cdot A^2\Omega^2$, and the matter outflowing from the $L_1^{\rm rot}$ point cannot get to distances larger than $R=9.1R_\odot\approx2R_1$. Hence, one has to expect formation of an envelope extending to a few stellar radii.

4.2 Gas-dynamical equations

To describe the gas flow, we have used 3-D gas-dynamical equations in cylindrical coordinates. We have modified the original conservative form of equations in cylindrical coordinates to obtain a system similar to the system of gas-dynamical equations in Cartesian coordinates. This approach permits us to treat the flow near the axis more accurately (see, e.g., Pogorelov et al. 2000). The corresponding equations are:

\begin{displaymath}\partial_t\rho +\partial_r\rho v_r +\frac{1}{r}\partial_\varphi\rho v_\varphi +\partial_z\rho v_z= -\frac{\rho v_r}{r}, \end{displaymath}


\begin{displaymath}\begin{array}{l} \partial_t\rho v_r +\partial_r\left(\vphanto... ...\frac{v_\varphi^2-v_r^2}{r} +2\rho\Omega v_\varphi, \end{array}\end{displaymath}


\begin{displaymath}\begin{array}{l} \partial_t\rho v_\varphi +\partial_r\rho v_r... ...Phi -\rho\frac{2v_r v_\varphi}{r} -2\rho\Omega v_r, \end{array}\end{displaymath}


\begin{displaymath}\begin{array}{l} \partial_t\rho v_z +\partial_r\rho v_r v_z +... ...\right)=-\rho\partial_z\Phi -\rho\frac{v_r v_z}{r}, \end{array}\end{displaymath}


\begin{displaymath}\begin{array}{l} \partial_t\rho E +\partial_r\rho v_r h +\fra... ...i -\rho v_z\partial_z\Phi -\rho\frac{v_r h}{r}\cdot \end{array}\end{displaymath}

Here $\partial_t\equiv\partial/\partial t$, $\partial_r\equiv\partial/\partial r$, $\partial_\varphi\equiv\partial/\partial\varphi$, $\partial_z\equiv\partial/\partial z$, $\rho$ is the density, ${\vec v}=(v_r, v_\varphi, v_z)$ the velocity vector, P pressure, $E=\varepsilon+\hbox{$~^{1}\!/_{2}$ }\vert{\vec v}\vert^2$ the full specific energy, $\varepsilon$ specific internal energy, and $h=\varepsilon+P/\rho+\hbox{$~^{1}\!/_{2}$ }\vert{\vec v}\vert^2$ specific full enthalpy. Gas-dynamical equations are written in the adopted co-rotating frame (i.e. in the frame where the centres of stars are in rest), so the Coriolis force is included in momentum equations. As usual, the system of gas-dynamical equations is closed using the equation of state. The equation of state of a perfect gas $P=(\gamma-1)\rho\varepsilon$ was adopted, $\gamma$ being the adiabatic index. In calculations, we adopted the value of $\gamma$ very close to 1 to mimic the energy losses due to radiative cooling (see, e.g., Sawada et al. 1986; Bisikalo et al. 1998). The Roche potential in cylindrical coordinates has the following form:

\begin{displaymath}\begin{array}{l} \Phi=-\frac{GM_1}{\sqrt{r^2+z^2}} -\frac{GM_... ...t(\frac{M_2}{M_1+M_2}\right)^2Ar\cos\varphi\right). \end{array}\end{displaymath}

4.3 The finite-difference model

To solve the system of gas-dynamical equations we used a monotonic Roe's scheme (Roe 1986) of first-order approximation with Osher's flux limiters (Chakravarthy & Osher 1985) that increases the order of approximation and leaves the scheme monotonous.

Gas flow was simulated in a cylinder $r\le1.1\cdot x_{L_1}=95R_\odot$, $0\le z\le 50R_\odot$ (because of symmetry with respect to the equatorial plane, calculations could only be conducted in the upper half-space). Non-uniform finite-difference grids (denser near the Be star and the equatorial plane) containing $50\times25\times30$ gridpoints on r, z, and $\varphi$, respectively, were used.

As for the initial condition, we adopted a rarefied gas with $\rho_0=10^{-6}$, P0=10-4, and ${\bf v_0}=0$.

The boundary conditions were defined as follows: in the gridpoints that correspond to $L_1^{\rm rot}$ we adopted the condition of injection of matter: $\rho(L_1^{\rm rot})=1$, $T(L_1^{\rm rot})=22~900~\mbox{K}$, which corresponds to the sound velocity $c(L_1^{\rm rot})=14~\mbox{km}~\mbox{s}^{-1}$, $v_r(L_1^{\rm rot})=c(L_1^{\rm rot})$, $v_\varphi(L_1^{\rm rot})=594~\mbox{km}~\mbox{s}^{-1}$, $v_z(L_1^{\rm rot})=c(L_1^{\rm rot})$. Note that an arbitrary value of the boundary density $\rho_0$ can be chosen since the system of equations can be scaled with respect to $\rho$ and P. To derive the true values of density in a specific system with a known mass loss rate, the calculated densities must simply be changed in accordance with the scale determined from the ratio of the true and model mass-loss rate. The boundary conditions were derived by solving the Riemann problem between the gas parameters ( $\rho_0,~{\vec v}_0,~P_0$) in $L_1^{\rm rot}$ point and the parameters in the computation gridpoint closest to it (see, e.g., Sawada et al. 1986; Sawada & Matsuda 1992, or Bisikalo et al. 1998). A full absorption of matter was assumed for the rest of the Be-star surface and for the outer boundary of the computational domain. We have verified that the outer boundary has virtually no effect on the results of computations. The boundary condition at the stellar surface is more important and less clear. However - considering the strong gravitational pull of the massive Be star and centrifugal acceleration, insufficient to cause a mass outflow - we believe that the assumption of the full absorption of matter is legitimate and a physically correct one.

The 3-D gas-dynamical calculations represent a very time-consuming task even for the most powerful present-day computers. That is why for this first investigation we followed the mass outflow only for the minimum interval of time needed. The calculations started as an injection into vacuum and converged to a more or less stable solution after about one third of the orbital period. They were continued to about one full orbital cycle after the emerging pattern of variations remained stable within the adopted accuracy. It is conceivable that future modelling, which we plan to continue until truly steady-state situation, will reveal a somewhat different structure of the envelope. However, it cannot alter our principal finding that the duplicity of a B star can lead to the formation of Be envelope even in a detached binary system.

5 Results of calculations

The six panels of Fig. 5 show the bird-eye view of density isosurfaces in the vicinity of the Be star at the level $\rho =5\times 10^{-5}\cdot \rho (L_1^{\rm rot})$. These isosurfaces are shown for six moments of time: $t=0.514P_{\rm orb}$, $0.545P_{\rm orb}$, $0.586P_{\rm orb}$, $0.608P_{\rm orb}$, $0.639P_{\rm orb}$, and $0.660P_{\rm orb}$, respectively. Projections of the envelope into XY plane (density distribution and velocity vectors in equatorial plane) and XZ plane (density distribution in the frontal plane) are also depicted in those six panels. An analysis of the results shows that the outflow of matter from the Be star in the vicinity of $L_1^{\rm rot}$ point results in the formation of an envelope with a fast retrograde apsidal motion. The mean angular velocity of apsidal motion can be derived from data of Fig. 6 where the time evolution of two angles, the angle between X axis and the direction to the centre of mass of the envelope (dashed line), and the angle between X axis and the direction to the centre of mass of the outer part of the envelope (solid line; values of density $\rho \in [10^{-3}, 10^{-2}]$ are shown).[*] It follows from the data of Fig. 6 that $\Omega_{\rm aps}=-6\Omega$, and $P_{\rm aps}=\hbox{$~^{1}\!/_{6}$ }P_{\rm orb}$, respectively. Moreover, it is clearly seen that the velocity of the apsidal motion is variable; the motion gets slower in the interval $\alpha\in[0.4\pi,0.6\pi]$. This region corresponds to the zone of interaction between the envelope and the outflowing stream of matter from $L_1^{\rm rot}$. It seems that this interaction causes a standstill of the apsidal motion. The analysis of Fig. 6 also shows that the centre of mass of the whole envelope oscillates within the interval $\alpha\in[0.35\pi,0.65\pi]$ while the centre of mass of the outer layers makes a full revolution within the interval $\alpha\in[-\pi,\pi]$. This finding seems to indicate the presence of a strong differential rotation of the envelope. It is interesting to note that the elongated shape of the envelope implies varying velocity of rotation, the velocity being larger than the Keplerian one in one part of the orbit and lower than Keplerian in the rest of the orbit. Figure 7 illustrates this phenomenon and depicts the ratio of angular velocity to the Keplerian one along the isolines $\rho =9\times 10^{-5}$ and $\rho =5\times 10^{-3}$ for the time $t=0.514P_{\rm orb}$ (see also Fig. 5a). The length of each bar characterizes the deviation of angular velocity from the Keplerian one $\nu =\vert V-V_{\rm Kep}\vert/V_{\rm Kep}$ (the bar in upper right corner corresponds to $\nu =1$) and its direction - inward or outward - denotes the sub-Keplerian or super-Keplerian flow regime, respectively.


  \begin{figure} \par\resizebox{\hsize}{!}{\includegraphics{m2949f6.eps}}\end{figure} Figure 6: Time variation of the angle between X axis and the direction to the centre of mass of the envelope (dashed line), and the angle between X axis and the direction to the centre of mass of the outer part of the envelope (with values of density $\rho \in [10^{-3}, 10^{-2}]$ - solid line). Six asterisks correspond to six moments of time presented in Fig. 5.
Open with DEXTER


  \begin{figure} \par\resizebox{\hsize}{!}{\includegraphics{m2949f7.eps}}\end{figure} Figure 7: The ratio of angular velocity to Keplerian one along the isolines $\rho =9\times 10^{-5}$ and $\rho =5\times 10^{-3}$ for the time $t=0.514P_{\rm orb}$. The length of each bar characterizes how much the angular velocity deviates from the Keplerian one $\nu =\vert V-V_{\rm Kep}\vert/V_{\rm Kep}$ (the bar in the upper right corner corresponds to $\nu =1$) and its direction - inward or outward - characterizes the sub-Keplerian or super-Keplerian flow regime, respectively. The dashed circle represents the surface of star 1. Coordinates x and y are expressed in $R_\odot $.
Open with DEXTER

The drawings of the envelope in Fig. 5 are given for six moments of time that cover the entire $P_{\rm aps}$. Our analysis provides, therefore, a complete picture of the flow and gives an estimate of the main parameters of the envelope as well. It is seen that for some time ( $t=0.514P_{\rm orb}$, $t=0.545P_{\rm orb}$, and $t=0.586P_{\rm orb}$ (which corresponds to $\alpha_{\rm out}\in[\pi,0.3\pi]$ for the vector pointed to the centre of mass of the outer layers of the envelope), the envelope has a torus-like shape, the thickness of the envelope being $h\sim R_1$, i.e. exceeding the polar radius which equals $\sim$0.65R1. In the rest of the time ( $t=0.608P_{\rm orb}$, $t=0.639P_{\rm orb}$, and $t=0.660P_{\rm orb}$ (corresponding to $\alpha_{\rm out}\in[0.3\pi,-\pi]$), one can see the elongation of the envelope. Its shape becomes nearly disk-like with a characteristic thickness $h\sim\hbox{$~^{1}\!/_{2}$ }R_1$. It is obvious that the changes in the envelope shape result from both the presence of the binary companion and the interaction of the gas in the envelope (during its apsidal motion) with the stream of gas leaving $L_1^{\rm rot}$. The characteristic linear sizes of the envelope in the equatorial plane are the following: for times when it has a torus-like shape it is $\sim$3R1 (on the level $\rho=5\times10^{-5}\times\rho(L_1^{\rm rot})$), while for the time when it has disk-like shape, its size increases to $\sim$4.5R1[*].

A test calculation for a binary composed from the same two stars but with an orbital period of 15 days also led to rapid retrograde revolution of the outflowing disk, though with a less regular pattern.

6 Discussion and future prospects

Gas-dynamical simulations of flow structure in binary Be stars has proven the existence of another possible mechanism for the formation of a Be-type envelope in binary Be stars: outflow of matter from the vicinity of "asynchronous'' inner Lagrangian point $L_1^{\rm rot}$ of a rapidly rotating B star. Note that for the values of macroscopic parameters, typical for Be binaries, the energy of the gas is insufficient to reach the inner Lagrangian point L1 of the binary, therefore there is no mass transfer between the binary components involved in the process. The outflowing gas only forms an envelope around the B star.

The main characteristics of such an envelope, according to our first numerical simulations, are:

1.
The envelope has a complex shape; is is not axially symmetric but elongated in one direction.
2.
Rapid cyclic changes of the envelope geometry between the disk-like and torus-like shape occur on a time scale almost an order of magnitude shorter than the binary orbital period.
3.
The envelope undergoes a fast retrograde apsidal motion, again on a time scale of a fraction of the orbital period, and the speed of this apsidal advance differs for the inner and outer parts of the envelope, indicating differential rotation inside the envelope.
Opponents of this interpretation can object that statistical studies showed that the Be stars are rotating well below their break-up velocities at the equator, which makes the mechanism proposed here inoperational. We do believe that the question of how close the true rotational velocities to the critical ones are needs to be carefully investigated for particular Be stars, not only on statistical grounds. Note, for instance, that Harmanec (2002) found that the observed v sin iof $\gamma$ Cas seems to agree with the expected break-up velocity for the best currently available estimates of the basic physical properties of this Be binary.

Undoubtedly, the complex variations of the envelope shape found here would strongly influence the observational appearance of Be stars in detached binaries. However - being aware of the limitations of our present gas-dynamical simulations - we are not attempting to compare our results with the observed variations for real detached Be binaries and postpone that for another study in this series, based on longer series of gas-dynamical modelling.

Here, we only offer a few thoughts on the potential of the new model:

1.
Our modelling should be complemented by calculations of emission-line profiles for the variable Be envelopes originating from the outflow discussed here. Yet, it seems that one should observe rapid changes in the intensity, width and V/R ratio of emission profiles, probably phase shifted for different lines, originating in different parts of the envelope. Even more stringent constraints may come from the analyses of forthcoming interferometric observations with VLTI focal instruments, such as a detection of phase-locked variations of the intensity maps at different wavelengths, both in the spectral lines and in continuum.
2.
We believe that the new model has the potential to explain the observed irregular emission-line episodes of Be stars. If a B star in a binary system is close to the limit of its rotational instability, even a small accidental disturbance can lead to the outflow of matter via the $L_1^{\rm rot}$ point. This can lead to a change in the total angular momentum of the system which, in turn, can either facilitate the outlow of matter or to stop it. One can assume that there should be a longer period of quiescence after a large emission episode than after a minor one.
3.
An interesting area for investigation is the case when the parameter of asynchronous rotation f is small, say 2 or so, since our analysis showed that this leads to conditions favourable for mass transfer between the binary components in systems which are still detached from the Roche lobe. Existence of such systems would substantially alter our current ideas about the process of mass transfer in binaries.
In the following papers of this series, we intend to address the above ideas via detailed modelling.

Acknowledgements
We thank the referee, Dr. Philippe Stee, for his constructive suggestions which helped us to improve and clarify some parts of the text. The use of the computerized bibliography from the NASA Astronomical Data System is also gratefully acknowledged. The study of PH was partly supported from the research plans J13/98 113200004 and AV 0Z1 003909 and from project K2043105. PH also acknowledges the support from the collaborative program KONTAKT ME402(2000) and CONACyT. D.B., A.B. and O.K. were partly supported by RFBR via grants 02-02-16088 and 00-15-96722, by FP "Astronomy", by program of RAS "Nonstationary Stars" and by INTAS via grant 00-491.

References

 

Copyright ESO 2002