A&A 395, 899-905 (2002)
DOI: 10.1051/0004-6361:20021337

The structure equations of contact binaries

H. Kähler

Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany

Received 3 April 2002 / Accepted 13 September 2002

The structure equations of contact binaries are discussed. An equation for the transfer of mass between the components is derived. Serious uncertainties concern only the transfer of energy. They are expressed as unknown functions. An evolution code for contact binaries is presented. The restrictions imposed by a spherically averaged treatment of the components are discussed. Apart from these restrictions the code can be adapted to any choice of the functions describing the energy transfer. A simple tentative choice for these functions is proposed.

Key words: stars: binaries: close

1 Introduction

Structure and evolution of contact binaries are complex and by no means well understood (for a review cf. Eggleton 1996). Sinse theoretical work on early-type systems is almost absent, the following discussion will be restricted to late-type systems. Basic observational facts (e.g. Rucinski 19851993) are the light curves, the period-colour relation, the broad range in mass, mass ratio and spectral type, the preference for shallow contact, and the absence of mass ratios close to unity. Most systems appear to be unevolved or slightly evolved only, and some systems are known to be unevolved.

Lucy (1968) had the key idea of a common convective envelope in which the entropy is constant and energy is transferred from the primary to the secondary. Thermal equilibrium turned out to be usually not possible. Lucy (1968) found age zero equilibrium models (violating the period-colour relation) only in a limited mass range and for extreme Population I composition. Moss & Whelan (1970) and Whelan (1972) encountered similar difficulties. Indeed, in unevolved systems with typical masses and typical composition thermal equilibrium was found to be in conflict with the equal entropy condition (Kähler 1995). If observed systems are assumed to be evolved, thermal equilibrium, when possible, turned out to be unstable (Kähler 1997a).

All these results are based on the assumption that the energy transfer between the components occurs in the adiabatic part of the common convective envelope (adiabatic transfer). The modification due to superadiabatic transfer is small (Whelan 1972; Biermann & Thomas 1973). Extreme superadiabatic transfer (Moss & Whelan 1973) gives more freedom, but Hazlehurst (1974) showed that this case must be excluded since the heat capacity of the secondary's subphotospheric layers is too small.

Abandoning the restriction to thermal equilibrium, Lucy (1976), Flannery (1976), and Robertson & Eggleton (1977) obtained contact binary solutions evolving in thermal cycles about a state of marginal contact. The solutions are in the required mass range and in agreement with the period-colour relation, but they have bad light curves since for a considerable part of the time the temperature difference between the components is too large. In a series of papers Rahunen and Vilhu encountered the same difficulty, the so-called light curve paradox, also for systems evolving in cycles without loss of contact. They concluded that the TRO (thermal relaxation oscillation) theory is not able to give a correct description of contact binaries (Rahunen 1983). The cycles investigated by Kähler et al. (1986) have indeed again bad light curves. In a recent paper Hazlehurst (2001) confirms the conclusion that the present theory of cyclic contact binaries is not in a position to resolve the light curve paradox.

As long as a proper hydrodynamic treatment of contact binaries is not possible, assumptions are necessary in any discussion. The results on thermal cycles obtained so far are based on few and apparently weak assumptions. The deficiencies of the resulting solutions have led to a number of investigations based on stronger assumptions in order to find contact binary models in thermal equilibrium which are compatible with the observations and sound from a theoretical viewpoint. All these attempts failed. In particular, the contact discontinuity hypothesis (Shu et al. 1976) violates the second law of thermodynamics (Hazlehurst & Refsdal 1978; Papaloizou & Pringle 1979), the assumption of an extended turbulent envelope in the primary (Kähler 1989; Kähler & Fehlberg 1991) is also in conflict with the second law (Hazlehurst 1993), and the assumption of sufficiently rapid internal mass motions in the primary's envelope (Kähler 1995) is in conflict with the observations (Kähler 1997a). For additional arguments against the possibility of thermal equlibrium see Hazlehurst (1993).

Collecting these results, thermal equilibrium appears to be as impossible as thermal disequilibrium. Nevertheless contact binaries are observed. We are thus led to the conclusion that the assumptions made so far in the theoretical discussion are not weak but too restrictive. This casts a doubt upon results which have been regarded as reliable for decades. In particular, we expect that if the restrictions are removed the light curve paradox for thermal cycles can be resolved and/or thermal equilibrium can be shown to be possible. The purpose of the present paper is to collect the structure equations of contact binaries on the basis of truly weak assumptions, to discuss the remaining uncertainties, and to present an evolution code for contact binaries which is flexible enough to be an appropriate tool when investigating the theoretical problems and trying to resolve the conflict between theory and observations.

2 The structure equations of contact binaries

We adopt a self-consistent spherically averaged treatment of the components (Kähler 19861997b) and assume uniform rotation and conservation of total mass M=M1+M2 and angular momentum J. Each component is assumed to be in hydrostatic equation, but the Roche equipotential condition is not assumed to be satisfied. The system may be either in contact or semi-detached or detached.


Table 1: Coefficients cjk describing the area of the neck in systems with $q\ge 0.1$.

c1k c2k c3k c4k c5k

$ 0.1875129 \times 10^{-0}$ $ 0.3983312 \times 10^{-1}$ $ 0.8526208 \times 10^{-2}$ $ 0.1033800 \times 10^{-2}$ $ 0.9590364 \times 10^{-3}$
1 $ 0.7118728 \times 10^{-1}$ $ 0.2322256 \times 10^{-1}$ $ 0.8475562 \times 10^{-2}$ $ 0.6828108 \times 10^{-3}$ $ 0.1776656 \times 10^{-2}$
2 $-0.2894722 \times 10^{-1}$ $-0.4160200 \times 10^{-2}$ $ 0.4443636 \times 10^{-3}$ $-0.6483157 \times 10^{-3}$ $ 0.8321587 \times 10^{-3}$
3 $ 0.1142797 \times 10^{-1}$ $ 0.1214730 \times 10^{-2}$ $-0.3399640 \times 10^{-3}$ $-0.2098207 \times 10^{-3}$ $-0.4379926 \times 10^{-4}$
4 $-0.4940387 \times 10^{-2}$ $-0.5475169 \times 10^{-3}$ $ 0.1630188 \times 10^{-3}$ $ 0.1086021 \times 10^{-3}$ $-0.2663696 \times 10^{-4}$
5 $ 0.1916213 \times 10^{-2}$ $ 0.2852279 \times 10^{-3}$ $-0.5841164 \times 10^{-4}$ $-0.4100397 \times 10^{-4}$ $ 0.1298322 \times 10^{-4}$
6 $-0.9130680 \times 10^{-3}$ $-0.1526878 \times 10^{-3}$ $ 0.1392672 \times 10^{-4}$ $ 0.1545703 \times 10^{-4}$ $-0.4117447 \times 10^{-5}$
7 $ 0.1851706 \times 10^{-2}$ $ 0.6550784 \times 10^{-4}$ $-0.9329447 \times 10^{-6}$ $ 0.2139876 \times 10^{-5}$ $ 0.1465290 \times 10^{-5}$
8 $-0.1245003 \times 10^{-2}$ $-0.2541965 \times 10^{-4}$ $-0.1047323 \times 10^{-5}$ $-0.3300985 \times 10^{-5}$ $-0.5642312 \times 10^{-6}$
9 $-0.1259457 \times 10^{-2}$ $ 0.3325350 \times 10^{-4}$ $ 0.2875518 \times 10^{-5}$ $-0.9470286 \times 10^{-5}$ $-0.5865679 \times 10^{-6}$
10 $ 0.9270583 \times 10^{-3}$ $-0.2953845 \times 10^{-4}$ $-0.3925625 \times 10^{-5}$ $ 0.5986984 \times 10^{-5}$ $ 0.7838701 \times 10^{-6}$
11 $ 0.8772300 \times 10^{-3}$ $ 0.5431091 \times 10^{-5}$ $ 0.2145996 \times 10^{-5}$ $ 0.5477786 \times 10^{-5}$ $-0.1043399 \times 10^{-6}$
12 $-0.6223828 \times 10^{-3}$ $ 0.2346910 \times 10^{-5}$ $-0.3850357 \times 10^{-6}$ $-0.3856733 \times 10^{-5}$ $-0.9468877 \times 10^{-7}$

2.1 The differential equations

Let $\omega$ be the angular velocity, mi the mass variable, gi the effective gravity in the component i, and $\sigma_i$the source (when positive) or sink (when negative) of energy per unit of mass caused by the interaction of the components. Using standard notation otherwise, the basic differential equations are
$\displaystyle \frac{\partial r_i}{\partial m_i}$=$\displaystyle \frac{1}{4\pi r_i^2 \rho_i},$ (1)
$\displaystyle \frac{\partial P_i}{\partial m_i}$=$\displaystyle -\frac{1}{4\pi r_i^2}\;g_i,\quad
g_i=\frac{{\rm G}m_i}{r_i^2}-\frac{2}{3}\omega^2 r_i,$ (2)
$\displaystyle \frac{\partial l_i}{\partial m_i}$=$\displaystyle \varepsilon_i-T_i \dot{s}_i+\sigma_i,$ (3)
$\displaystyle \frac{\partial T_i}{\partial m_i}$=$\displaystyle \frac{T_i}{P_i}\frac{\partial P_i}{\partial m_i} \nabla_i.$ (4)

Concerning the temperature gradient, let $\nabla_{{\rm ad},i}$ be the adiabatic gradient,

\nabla_{{\rm rad},i}=\frac{3}{16 \pi {\rm acG}}
\frac{\kappa_i l_i P_i}{m_i T_i^4}
\end{displaymath} (5)

the radiative gradient, and $\nabla_{{\rm t},i}$ the gradient in the presence of turbulence. Note that turbulence can be caused by convective instability ( $\nabla_{{\rm rad},i}>\nabla_{{\rm ad},i}$) as well as by circulation currents. The gradient $\nabla_{{\rm t},i}$ can be obtained from a mixing-length theory if $\nabla_{{\rm rad},i}>\nabla_{{\rm ad},i}$, and from a generalised mixing-length theory (Kähler 1989) otherwise. Let $\alpha_+$ be the mixing length and $\alpha_-$ the generalised mixing length, both in units of the pressure scale height. With these definitions we have

\nabla_i=\left\{ \begin{array}{ll}
\nabla_{{\rm rad},i}\qua...
...{{\rm t},i} & \mbox{in turbulent layers}.
\end{array} \right.
\end{displaymath} (6)

Convectively stable turbulent layers are expected in the common turbulent envelope of early-type contact binaries, and possibly also in late-type systems. Concerning the extent of these layers some assumption is necessary, which will be labeled by the number $i_{\rm t}$. We may assume that convectively stable turbulent layers are absent ( $i_{\rm t}=0$). Alternatively we may assume that the extent of these layers is determined by the condition that layers with energy sources and sinks caused by circulation and accretion are turbulent ( $i_{\rm t}=1$). Other assumptions (to be labeled with $i_{\rm t}>1$) are also possible.
\end{figure} Figure 1: The function $\beta (q)$ describing the area of the neck in systems in shallow contact (see text).
Open with DEXTER

Since mass M and angular momentum J are conserved, the global balance of energy (Kähler 1997b) requires that

\sum_{i=1}^2 \int_0^{M_i} \sigma_i ~{\rm d}m_i=L_{\rm acc},
\end{displaymath} (7)


L_{\rm acc}=(\Psi_1-\Psi_2+h_1-h_2)\dot{M_2}
\end{displaymath} (8)

is the accretion luminosity. The quantities hi and $\Psi_i$ denote, respectively, specific enthalpy and potential at the surface of the component i. In the presence of mass exchange ( $\dot{M}_2\neq 0$) let $i={\rm g}$ in the gainer (the mass gaining component), and $i={\rm l}$ in the loser. Since the accretion luminosity is caused by energy sources in the gainer's outermost layers, we can write

\sigma_i=\sigma_{i,{\rm circ}}+\delta_{i{\rm g}}\sigma_{\rm acc}
\end{displaymath} (9)


\int_0^{M_{\rm g}} \sigma_{\rm acc}~{\rm d}m_{\rm g}=L_{\rm acc},
\end{displaymath} (10)

and the luminosity transferred by circulation currents from the primary to the secondary is

\Lambda=-\int_0^{M_1}\sigma_{1,{\rm circ}}~{\rm d}m_1
= \int_0^{M_2}\sigma_{2,{\rm circ}}~{\rm d}m_2.
\end{displaymath} (11)

2.2 The boundary conditions

The central boundary conditions are

m_i=0: \quad r_i=l_i=0.
\end{displaymath} (12)

The outer boundary conditions are in a close approximation
$\displaystyle m_i= m_{{\rm f},i}:\quad P_i$ = $\displaystyle \tilde{P}_i(M_i,L_i,R_i,\omega)$ (13)
Ti = $\displaystyle \tilde{T}_i(M_i,L_i,R_i,\omega)$ (14)
ri = $\displaystyle \tilde{r}_i(M_i,L_i,R_i,\omega)$ (15)
li = Li, (16)

where the functions $\tilde{P},\tilde{T},\tilde{r}$ are defined by outer layer integrations, assuming li=Li, and the fitting mass

m_{{\rm f},i}=x_{\rm f}M_i,\qquad x_{\rm f}=1-10^{-11}
\end{displaymath} (17)

is very close to the surface, just below the photosphere. This choice makes it possible to treat systems with energy sources or sinks in the subphotospheric layers.

2.3 Conservation of mass and angular momentum

Total mass M, angular velocity $\omega$ and the separation A of the components are connected by Kepler's law

\omega^2={\rm G}M/A^3.
\end{displaymath} (18)

Conservation of mass requires

M1+M2=M, (19)

and conservation of angular momentum requires

\end{displaymath} (20)

with the moments of inertia

\Theta_0=\mu A^2, \quad \mu=M_1M_2/M
\end{displaymath} (21)


\Theta_i=\frac{2}{3}\int_0^{M_i}r_i^2~{\rm d}m_i \quad(i=1,2).
\end{displaymath} (22)

2.4 The Roche potential and the effective potential

The components are treated as spherical stars. In addition to the Roche potential we have therefore in each component an effective central potential. In a self-consistent treatment this potential is well-defined (see Kähler 1986). The surface potential of the component i is

\Psi_i=-\frac{{\rm G}M}{A}\left( \frac{\lambda_i^2}{3}
+\frac{\mu_i}{\lambda_i} + \frac{(1-\mu_i)(3-\mu_i)}{2}\right),
\end{displaymath} (23)

where $\mu_i$ and $\lambda_i$ are, respectively, normalised mass and radius of the component defined by

\mu_i=M_i/M, \quad \lambda_i=R_i/A.
\end{displaymath} (24)

Note that $\Psi_i$ is the average value of the Roche potential over the spherical surface. The degree of contact Fi of the component i is defined by

F_i=(\Psi_i-{\Psi_{\rm L}}_1)/({\Psi_{\rm L}}_2-{\Psi_{\rm L}}_1),
\end{displaymath} (25)

where ${\Psi_{\rm L}}_n$ denotes the Roche potential at the Lagrangian point n. The effective central potential $\psi_i(r_i)$ in the component i is the solution of the differential equation

\frac{{\rm d}\psi_i}{{\rm d}r_i}=g_i
\end{displaymath} (26)

with the property $\psi_i(R_i)=\Psi_i$. In the outermost layers of the components the effective central potential almost coincides with the spherically averaged Roche potential. In the interior the two potentials differ and only the effective potential is useful. We shall also need the normalised potential $\varphi_i$ defined by

\psi_i=-\frac{{\rm G}M}{A}\varphi_i.
\end{displaymath} (27)

2.5 The extent of the neck

If a system in contact satisfies the Roche equipotential condition $\Psi_1=\Psi_2$, the area B of the neck between the components can be written in the form

B=A2 b(F,q), (28)

where b is a function of mass ratio q=M2/M1 and degree of contact F. A close approximation for $q\ge 0.1$ is

b=\sum_{j=1}^5c_j~F^j, \quad c_j=\sum_{k=0}^{12}c_{jk}~y^k, \quad
\end{displaymath} (29)

with coefficients listed in Table 1. For systems in shallow contact ($F\ll 1$) we have

\end{displaymath} (30)

with the function $\beta (q)$ shown in Fig. 1.

The area of the neck can be used can be used to derive constraints for the exchange of mass and energy. We shall need the function Bidefined by

B_i=\left\{ \begin{array}{ll}
A^2~b(F_i,q)\quad & \mbox{if ~ $F_i>0$ ,} \\
0 & \mbox{otherwise}.
\end{array} \right.
\end{displaymath} (31)

2.6 The transfer of mass

If the system is semi-detached, matter in the neck is streaming from the loser to the gainer, approximately with the velocity of sound, and through an area which is smaller than $B_{\rm l}$ (cf. Fig. 1.6 of Pringle 1985). Accordingly we have

\dot{M}_{\rm g}<B_{\rm l}\rho_{\rm l}c_{\rm s,l},
\end{displaymath} (32)

where $\rho_{\rm l}$ and $c_{\rm s,l}$ denote, respectively, density and sound velocity at the loser's inner critical surface. Since $c_{\rm s}\approx\sqrt{P/\rho}$ and since $B_{\rm g}=0$ in semi-detached systems we can write

\dot{M}_{\rm g}=f_{\rm M}(B_{\rm l}-B_{\rm g})\sqrt{\rho_{\rm l}P_{\rm l}},\quad f_{\rm M}<1.
\end{displaymath} (33)

Following Paczynski & Sienkiewicz (1972), Pringle (1985) showed that

\dot{M}_g=f_{\rm P}(\rho c_{\rm s}^3)_{\rm l} \omega^{-2}
\end{displaymath} (34)

with $f_{\rm P}$ of order unity. As expected, Eqs. (33) and (34) turned out to be roughly equivalent. In most cases investigated equivalence requires $f_{\rm M}/f_{\rm P}\sim 0.3\ldots0.4$.

We need an equation for the rate of mass exchange which depends continuously on the structure of the system and which is physically plausible in all cases. Equation (33) has these properties. If matter is streaming almost freely from the loser to the gainer, the rate can be written as the product of density, sound velocity and effective area in the neck and thus in the form of Eq. (33). This is the case not only in semi-detached systems but also in systems in contact as long as $B_{\rm g}\ll B_{\rm l}$. In detached systems without mass exchange the equation is trivially valid.

The general case of systems in contact is more complex since matter is streaming from the loser to the gainer (at the rate $\mu_{\rm g}$, say) as well as from the gainer to the loser (at the rate $\mu_{\rm l}$). Each rate separately can be approximated by a product of characteristic values for density, velocity (usually smaller than the sound velocity) and effective area. For the net rate $\dot{M}_{\rm g}=\mu_{\rm g}-\mu_{\rm l}$however an expression of this form is not obvious, except for the case (already mentioned) that $B_{\rm g}\ll B_{\rm l}$ ans thus $\mu_{\rm l}\ll \mu_{\rm g}$. Nevertheless Eq. (33) is plausible as an approximation since providing fulfilment of the Roche equipotential condition in a close approximation, with the possible exception of short phases during a cycle just before contact is broken or just after contact has been reestablished. As long as $B_{\rm g}\ll B_{\rm l}$ in these phases, Eq. (33) is certainly adequate. For these reasons Eq. (33), with $f_{\rm M}$ taken in the range $0.1\ldots0.3$, is expected to be reasonable in any system.

2.7 The transfer of energy

With the rate of mass exchange also the accretion luminosity  $L_{\rm acc}$ is known from Eq. (8). The luminosity $\Lambda$ transferred by circulation currents however is still unknown.

It is convenient to use the variable xi=mi/Mi and to write the instantaneous distribution of the energy sources and sinks in the form

\sigma_{\rm acc} = L_{\rm acc}~\vartheta_{\rm acc}(x_{\rm g}),
\end{displaymath} (35)

\begin{displaymath}\sigma_{i,{\rm circ}}= (-1)^i\Lambda~\vartheta_i(x_i)\quad (i=1,2)
\end{displaymath} (36)

with normalised functions, i.e.

\int_0^{x_{\rm f}}\vartheta_{\rm acc}(x_{\rm g})~
{\rm d}x_{\rm g}=\int_0^{x_{\rm f}}\vartheta_i(x_i)~{\rm d}x_i=1.
\end{displaymath} (37)

So far it has almost always been assumed that the energy sources and sinks occur only in the outer layers. For this reason it is useful to introduce a simple and smooth normalised function which is zero in the interior and positive in the outer layers. The function $\vartheta(a,x)$ defined by

\vartheta=\left\{ \begin{array}{ll}
...x>a$ }\\
0 & \quad \mbox{otherwise} \\
\end{array} \right.
\end{displaymath} (38)


0<a<x_{\rm f},\quad V=1/(x_{\rm f}-a)
\end{displaymath} (39)

has these properties. Since the sources caused by accretion occur in the gainer's surface layers, we may assume that

\vartheta_{\rm acc}(x_{\rm g})=\vartheta(x_{\rm acc},x_{\rm g}),
\end{displaymath} (40)

where $x_{\rm acc}$ is a number very close to $x_{\rm f}$.

Sources and sinks caused by circulation currents occur only in systems in contact. Since the fitting mass is very close to unity, the contact condition

R_i>R_{{\rm c},i}\quad (i=1,2)
\end{displaymath} (41)

(where $R_{{\rm c},i}$ is the inner critical radius of the component i) can be approximated by

r_i(x_{\rm f})>R_{{\rm c},i}\quad (i=1,2).
\end{displaymath} (42)

It is usually assumed that energy sources/sinks occur only in the secondary's/primary's outer layers. If this assumption is reasonable, we can approximately write

\end{displaymath} (43)

where ai is a number somewhat smaller than $x_{\rm f}$. The assumption is however possibly too restrictive. For example, circulation currents deep in the interior of the components as proposed by Webbink (1977) are possibly associated with energy sources and sinks. As a further example, Martin & Davey (1995) presented arguments for the existence of additional energy sources/sinks in the primary's/secondary's layers near the critical surface. If these arguments are valid, the functions $\vartheta_i(x_i)$ change sign and additional parameters may be necessary. In any case we may assume that the function $\vartheta_i(x_i)$ is approximately described by few parameters $a_i,b_i\ldots$

The uncertainties concerning the energy transfer by circulation currents have been expressed as uncertainties in the quantities $\Lambda$ and $a_i,b_i\ldots(i=1,2)$. We shall assume that these quantities are (in way so far unknown) determined by the instantaneous spherically averaged structure of both components and thus by the set of functions

r_i(m_i),~P_i(m_i),~l_i(m_i),~T_i(m_i)\qquad (i=1,2),
\end{displaymath} (44)

which will be comprised in the symbol $\Gamma$. This assumption gives the largest freedom within the framework of a spherically averaged description of the components. Symbolically the functionals $\Lambda,a_i,\ldots$ will be written as functions

\Lambda(\Gamma),~a_i(\Gamma),~b_i(\Gamma),\ldots \quad (i=1,2).
\end{displaymath} (45)

To be more specific on this point, if the physical parameters (mass, angular momentum and composition) and the functions (44) are given, then from their range of definition also M1 and M2 are known. Since the functions are not independent but coupled by differential equations, also the angular velocity $\omega$ is known from the equation of hydrostatic equilibrium, and the separation A follows from Kepler's law. The instantaneous structure of the system (within the framework of a spherically averaged description of the components) is therefore completely known. Accordingly, $\Gamma$ represents the instantaneous structure. If $\Gamma$ is given, the rate of mass exchange and the effects of accretion are known. The evolution of the system (i.e. the time derivative $\dot{\Gamma}$) however is only determined if in addition the energy sources and sinks associated with circulation currents are given. Assuming that $\Lambda,a_i,b_i\ldots$ are still unspecified functions of $\Gamma$ we have the largest possible freedom. An example for functions of this type is the "transport equation''

\Lambda=Kd^m\Delta s^n
\end{displaymath} (46)

proposed by Hazlehurst & Refsdal (1980) and adopted by Rahunen & Vilhu (see Rahunen 1983), where d is the depth of contact (determined by the radii and the critical radii of the components), $\Delta s$ is the entropy difference between the two outer convective zones, and K,m,n are parameters.

3 Uncertainties in the structure equations

The equations derived in the preceding subsections contain the parameters $\alpha_+,\alpha_-,i_{\rm t},f_{\rm M},x_{\rm acc}$and the unknown functions listed in Eq. (45). The presence of these parameters and functions reflects the difficulties of the hydrodynamical problem. Concerning the choice of the parameters there is little freedom. We decided to use as standard values

\alpha_+=1.5,\;\alpha_-=1,\;f_{\rm M}=0.1,\;x_{\rm acc}=1-10^{-10}
\end{displaymath} (47)

and to compare results obtained with different values of $i_{\rm t}$. The uncertainties in the functions however are serious. We have assumed that the sources and sinks of energy caused by circulation are determined by the spherically averaged structure of the components. Is this assumption reasonable from the viewpoint of hydrodynamics? In the general case (i.e. when considering the evolution of a system starting from an arbitrarily given initial configuration) the answer is uncertain since the sources and sinks depend on details of the internal mass motions which escape notice in a spherically averaged description.

When studying the thermal evolution of unevolved systems, however, we are interested only in the long-time behaviour, which is expected to be either a stable limit-cycle (i.e. a thermal cycle) or a fixed point (i.e. a configuration in stable thermal equilibrium). More complex attractors cannot be excluded, but previous results (Kähler et al. 1986) suggest that they are not realised. Under these circumstances the problem is less difficult and a spherically averaged treatment of the components is reasonable.

To see this, consider first a given system evolving in relaxed thermal cycles. All properties (e.g. $\Lambda$) are periodic functions of the time t, and t can certainly be eliminated by inverting functions describing the spherically averaged structure of the components as a function of t. For example, if the primary's degree of contact F1(t) is a monotonic function in a certain time interval, then t=t(F1) and thus $\Lambda=\Lambda(F_1)$ in this interval. Accordingly, for a given system evolving in cycles the functions (45) can be defined in a variety of ways. We are interested in a physical definition which is (approximately) valid not only for the particular system under consideration but for all systems evolving in cycles as well for all systems in (stable or unstable) thermal equilibrium. Such a definition can hopefully be obtained in the future from 3-D hydrodynamic calculations. In the meantime we have to make plausible assumptions and to introduce free parameters. Equation (46) is an example.

Suppose next that thermal equilibrium can be achieved. In an approximate discussion the uncertainties in the parameters listed in Eq. (47) can be neglected. The structure of a configuration is then determined by the values of $\Lambda,~a_i\ldots \; (i=1,2)$, i.e. by few numbers, not functions. The investigation of possible equilibrium configurations is therefore much less difficult than the investigation of thermal cycles. In fact, it is possible to investigate all candidates for reasonable equilibrium configurations. We are interested in configurations which are dynamically as well as thermally stable. The dynamical stability problem is well determined. More difficult is the thermal stability problem since it involves the change of the quantities $\Lambda,~a_i\ldots$ in the course of thermal perturbations, and thus, if these quantities can be treated as functions of $\Gamma$, the derivatives of these functions.

4 A simple model

Some choice of the functions in Eq. (45) being necessary, here we define a simple model for late-type contact binaries involving only few parameters in addition to those in Eq. (47). We assume that the energy flux in the neck is not larger than

\rho c_{\rm s} T(s_1-s_2)\approx (\rho P)^{\frac{1}{2}}T(s_1-s_2),
\end{displaymath} (48)

with all quantities taken at the critical surface and quantities without index averaged between the components, and that therefore

\Lambda=f_{\rm E} B_2
(\rho_1 P_1 \rho_2 P_2)^{\frac{1}{4}}(T_1 T_2)^{\frac{1}{2}}(s_1-s_2),
\end{displaymath} (49)

where $f_{\rm E}$ is a parameter not larger than unity. Concerning the extent of the energy sources and sinks, let $x_i=x_{{\rm c},i}$ at the critical radius of a system in contact. We assume that Eq. (43) is valid and that

a_i = x_{\rm f}-\alpha_i(x_{\rm f}-x_{{\rm c},i}),
\end{displaymath} (50)

where $\alpha_i$ is a parameter. In other words, we make the usual assumption that energy sources/sinks occur only in the secondary's/primary's envelope, and assume that the extent of the sources or sinks is correlated with the extent of the layers above the inner critical surface.

In view of the persistent difficulties in the theoretical treatment of contact binaries we cannot expect that this model gives results which are fully compatible with the observations as well as with theoretical constraints. We do, however, expect that the deficiencies of the results can be used to reduce the inconsistencies and to improve the model. In other words, the model is intended as a starting point in the attempt to find a consistent description of contact binaries.

5 The evolution code

The equations collected above have been cast into an evolution code for contact binaries. The structure equations for both components are solved simultaneously. The code can be adapted to any choice of the functions in Eq. (45).

The non-Lagrangian mass variables xi=mi/Mi are used. Since the components exchange mass, this choice minimises the changes in pressure and temperature at the meshpoints near the surface during a time step.

As expected, numerical stability turned out to be sometimes a serious problem. While in quiet evolutionary phases time steps of more than 103 years can be used, in some difficult phases steps of 10-4 years are required. We found it necessary to treat the difference equations as proposed by Sugimoto (1970).

A problem with the temperature gradient concerns only the secondary's envelope. In some phases during a cycle the extent of the convective layers may not converge in the iteration process. Some layers oscillate between radiative transport and convection. This difficulty is removed if the extent of radiative and of turbulent layers is artificially kept fixed after a certain number of iterations.

Concerning the choice of an initial configuration, one of the equations in the Henyey scheme deals with the exchange of mass between the components. If this equation is replaced by an equation prescribing the mass ratio, and if artificial energy sources/sinks are introduced in the secondary's/primary's interior, systems in thermal equilibrium with any given mass ratio can easily be constructed. The energy sources and sinks can be adjusted to obtain a configuration in contact satisfying the Roche equipotential condition. Switching then off the artificial energy sources and sinks we obtain a contact configuration in thermal disequilibrium which, if dynamically stable, can be used as initial configuration.

6 Concluding remarks

We have collected the equations for zero age contact systems, assuming conservation of angular momentum. (In our opinion it is premature to study the effects of nuclear evolution and loss of angular momentum as long as the structure of unevolved systems is not well understood.) Are these equations based on truly weak assumptions? The continuity equation is beyond doubt. The equation of hydrostatic equilibrium (2) is a very close appoximation in the interior, and a close approximation also in the outer layers if (i) the velocity of the internal mass motions (in a rotating frame) is not larger than the the sound velocity and thus much smaller than the orbital velocity and (ii) the turbulent pressure is unimportant. If the velocity of the internal mass motions cannot be neglected, modifications of the equation of hydrostatic equilibrium and of the Roche equipotential condition will be necessary as described by Kähler (1995), and this requires a modification of Eq. (33). There are also uncertainties in the transport equation which can be estimated comparing results for different assumptions (i.e. different values of $i_{\rm t}$). Serious uncertainties concern only the energy balance. Some quantities decribing the energy transfer between the components are unknown. When investigating thermal cycles, these quantities can be treated as (so far unknown) functions depending on properties of the spherically averaged structure of the components. An evolution code has been presented which can be adapted to any choice of these functions. With this freedom and with the possibility to introduce other modifications (concerning the temperature gradient and departures from hydrostatic equilibrium in the outer layers) the treatment of contact binaries is free from restrictive assumptions.

The transfer of mass and the exchange of energy between the components have been treated separately. In reality these processes are coupled since determined by the internal mass motions. A combined treatment will be possible only when we get more insight into these motions. A simple model has been proposed as a starting point for more detailed investigations, and in an accompanying paper first results are presented.


Copyright ESO 2002