A&A 445, 647-652 (2006)
DOI: 10.1051/0004-6361:20053288

Numerical stability of mass transfer driven by Roche lobe overflow in close binaries

A. Büning - H. Ritter

Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany

Received 22 April 2005 / Accepted 24 August 2005

Abstract
Numerical computation of the time evolution of the mass transfer rate in a close binary can be and, in particular, has been a computational challenge. Using a simple physical model to calculate the mass transfer rate, we show that for a simple explicit iteration scheme the mass transfer rate is numerically unstable unless the time steps are sufficiently small. In general, more sophisticated explicit algorithms do not provide any significant improvement since this instability is a direct result of time discretization. For a typical binary evolution, computation of the mass transfer rate as a smooth function of time limits the maximum tolerable time step and thereby sets the minimum total computational effort required for an evolutionary computation. By methods of "Controlling Chaos'' it can be shown that a specific implicit iteration scheme, based on Newton's method, is the most promising solution for the problem.

Key words: binaries: close - stars: evolution - stars: mass-loss - methods: numerical

   
1 Introduction

To compute the long-term evolution of a semi-detached binary, the numerically determined mass transfer rate, which is a function of current stellar and binary parameters, must be a sufficiently smooth function of time. The simplest approach to obtain the mass M₂ of the mass losing secondary star is an explicit forward integration in time:

$$% M_2(t_{{n+1}})=M_2(t_{{n}})+\dot{M}_{2}(t_{{n}})~ \Delta t.$$ (1)

Here, $\Delta t = t_{{n}+1} - t_{{n}}$ denotes the length, and t_n and t_n+1 the start and end points of the nth time step. Every change of M₂ within one time step causes changes of other stellar and binary parameters which affect $\dot{M}_2$, and an additional change of M₂ would be necessary to obtain a consistent result for the current time step. This feedback between $\dot{M}_2$ and other stellar and binary parameters, especially the radius R₂ and the critical Roche radius  $R_{{\rm R},2}$ of the donor star, can destabilize the numerically computed mass transfer rate: $\dot{M}_2$ can become non-continuous between successive time steps and can even show strong fluctuations around its secular mean value. These numerical effects have been known for quite a long time (numerical experience by the authors; U. Kolb, K. Schenker, priv. comm.). Not surprisingly, only very few evolutionary calculations that show these instabilities of the mass transfer rate have been published, e.g., Kolb & Ritter (1990, Figs. 3-5), Sarna (1992, Figs. 1 and 9), D'Antona (1994, Fig. 2), Kolb et al. (2000, Figs. 3 and 4), and Schenker et al. (2002, Fig. 5).

Attempts to use a different explicit integration scheme that takes into account not only $\dot{M}_2(t_{{n}})$ but also $\dot{M}_2(t_{{n-j}})$ for certain values of j (e.g., for a particular average over the last m time steps) did not show major improvements. The reason for this behaviour of the mass transfer rate has remained unknown.

The main purpose of this paper is to show analytically why these numerical instabilities exist, what they are, which methods are suitable to suppress them, and which ones are not. Even in the case of the proposed implicit algorithm, calculation of the mass ransfer rate still limits the maximum tolerable time step in a numerical computation and thus sets the minimum total computational effort required for carrying out such a binary evolution. To illustrate this point: by imposing a fixed mass loss rate on a single low-mass main sequence star, up to 10% of the total mass of the star can be removed per time step (eventually after 1 or 2 initial time steps with a lower rate), and the stellar model still converges. But if the mass loss rate is coupled to the binary parameters, even in the case of our proposed implicit iteration scheme, typically no more than a few 10^-3 of the stellar mass can be removed per time step without losing convergence. When using the explicit algorithm (1), the corresponding value is much lower: often only about 10^-5-10^-4 or even less of the total mass can be removed per time step if fluctuations of the mass transfer rate by up to several orders of magnitude are to be avoided.

This paper is organized as follows: first, we discuss in Sect. 2 the necessary input physics before we show in Sect. 3 that the mass transfer rate in our model is physically stable. Subsequently, in Sect. 4.1 we prove mathematically that for a simple explicit algorithm the time-discretized mass transfer rate becomes unstable if $\Delta t$ is greater than a critical value. We also give an estimate of how many time steps are required to calculate the binary evolution. Next, we show in Sect. 4.2 that for increasing the mass transfer rate undergoes a Feigenbaum scenario Feigenbaum 1978,1980; for a more accessible review see, e.g., Thompson & Stewart 1986 and finally becomes chaotic. In Sect. 4.3 we discuss how the onset of chaos can be suppressed in terms of "Controlling Chaos'' and we present an implicit iteration scheme for the integration of the mass transfer rate. Finally, in Sect. 5 we discuss a number of points which have to be considered for a practical implementation.

   
2 Input physics

In this paper we use the same nomenclature as in Büning & Ritter (2004), hereafter called Paper I. For the mass transfer rate we use

$$% \dot{M}_{2}=-\dot{M}_{0}~ \exp \left( \frac{\Delta R}{H_{{\rm P}}}\right),$$ (2)

where

$$% \Delta R:=R_{2}-R_{{\rm R,2}}$$ (3)

denotes the difference between the radius R₂ and the Roche radius  $R_{{\rm R},2}$ of the donor, $H_{{\rm P}}$ the photospheric pressure scale height, and $\dot{M}_0 > 0$ a weakly varying function of several system parameters (for details, see Ritter 1988).

Instead of Eq. (2) other relations between  $\dot{M}_2$ and $\Delta R$ have also been used in the literature. For example Tout et al. (1997) have adopted a power-law dependence $\dot{M}_{2} \propto (\Delta R/R_2)^3$. However, a mass transfer prescription other than Eq. (2) results mainly in a different characteristic scale length

$$% H := \left( \frac{{{\rm d}}\ln \left(-\dot{M}_2\right)}{{{\rm d}}\Delta R} \right)^{-1}$$ (4)

at the secular mass transfer rate $\overline{\dot{M}}_2$ which is given by

$$% \tau_{{\rm M}} := - \frac{M_2}{~\overline{\dot{M}}_2} = (\zeta_{{\rm s}} - \zeta_{{\rm R}}) \tau_{{\rm d}}'.$$ (5)

Here, $\tau_{{\rm M}}$ is the secular time scale of the mass loss, $\tau _{\rm d}'$ the driving time scale including thermal relaxation of the donor, $\zeta_{{\rm s}}$ and  $\zeta_{{\rm R}}$ are respectively the adiabatic mass radius exponent of the donor and the mass radius exponent of the Roche radius, where $\zeta_{\rm s} > \zeta_{\rm R}$ is required for physical stability of mass transfer (for details, see Paper I). For the mass transfer prescription (2), we have $H \equiv H_{{\rm P}}$.

   
3 The time-continuous system

The mass transfer rate $\dot{M}_2$ is uniquely invertibly coupled to $\Delta R$ by Eq. (2). For simplicity, we will consider the time evolution of $\Delta R$ instead of $\dot{M}_2$ which is given by Eq. (36) of Paper I:

$$% \frac{{{\rm d}}}{{{\rm d}}t} \Delta R = (\zeta_{{\rm s}} - ... ...(\Delta R)}{M_2} + \frac{R_2}{\tau_{{\rm d}}'} =: F(\Delta R).$$ (6)

Since for our stability considerations we neither take into account changes in the system parameters ( $\zeta_{{\rm s}}$, $\zeta_{{\rm R}}, \tau_{{\rm d}}'$) nor irradiation of the donor as in Paper I, our model is completely described by this 1-dim. autonomous differential equation which simplifies the linear stability analysis significantly.

At the stationary value $\overline{\Delta R}$ which is the only fixed point (FP) of Eq. (6) and which is equivalent to the stationary mass transfer rate, using Eq. (2) we get

$$% DF(\overline{\Delta R}) := \frac{{{\rm d}}F(\overline{\Delt... ...frac{R_{2}}{H_{{\rm P}}} \frac{~\overline{\dot{M}}_{2}}{M_{2}}$$ (7)

(cf. Eq. (44) of Paper I). Since $\dot{M}_2 < 0$ and $\zeta_{\rm s} > \zeta_{\rm R}$, DF is negative, not only for $\overline{\Delta R}$, but even for all $\Delta R \in {\rm I\!R}$. This means that the stationary mass transfer rate, i.e., $\overline{\Delta R}$ is stable and that all solutions  $\Delta R(t)$ of Eq. (6) converge to  $\overline{\Delta R}$. The convergence occurs on a time scale of

$$% \tau \approx \frac{H_{{\rm P}}}{R_2} \tau_{\rm d}'$$ (8)

since

$$% \Delta R(t) = \overline{\Delta R} + \left[ \Delta R(0) - \overline{\Delta R} \right]~ {\rm e}^{-DF~t}$$ (9)

in the linearized system. This has also been discussed in more detail by D'Antona et al. (1989).

   
4 The time-discretized system

   
4.1 The explicit algorithm

The simplest way to obtain the donor mass numerically as a function of time is an explicit forward integration of $\dot{M}_2$ as given by Eq. (1). In this time-discretized system the evolution of $\Delta R$ is given by an iteration equation which follows directly from Eq. (6):

 

$$% \Delta R_{{n+1}} \Delta R_{{n}} + \left[ \frac{R_{2}}{\tau_{{\rm d}}'} + \left( \z... ...rm R}} \right) \frac{R_{2}}{M_{2}} \dot{M}_{2}(\Delta R_{n}) \right] ~ \Delta t$$ =

$$\Phi \left( \Delta R_{{n}} \right).$$ =: (10)

Here, $\Delta R_{n}$ denotes the value of $\Delta R$ at t = t_n which is in practice a (numerical) approximation for the "real'' value of $\Delta R$ in the time-continuous system.

Obviously, Eqs. (6) and (10), i.e., the time-continuous and the time-discretized systems have the same FP since $F(\Delta R) = 0$ if and only if $\Phi(\Delta R) = \Delta R$. But, since $\Phi(\Delta R)$ is a map while $F(\Delta R)$ is a vector field, the condition for stability is different: according to the Hartmann-Grobmann theorem (e.g., Guckenheimer & Holmes 1983), a FP of a map is stable if the absolute values of all eigenvalues of the Jacobi matrix of at are less than unity, and is unstable if at least one eigenvalue has an absolute value greater than unity.

In the case of $\Phi(\Delta R)$ from Eq. (10) this means: $\overline{\Delta R}$ is stable if

$$% \left\vert D\Phi \right\vert := \left\vert \frac{{{\rm d}}\... ...frac{~ \overline{\dot{M}}_{2}}{M_{2}} \Delta t\right\vert < 1,$$ (11)

and unstable if the absolute value is greater than unity. As can be easily shown by using Eq. (5), Eq. (11) is equivalent to

$$% \Delta t < 2 \frac{H_{{\rm P}}}{R_{2}} \tau_{{\rm d}}' =: \Delta t_{{\rm max}}.$$ (12)

This means that in the neighbourhood of $\overline{\Delta R}$ all solutions of the explicit iteration scheme (1) converge to the FP if $\Delta t$ is less than the critical time step length  $\Delta t_{{\rm max}}$. But for $\Delta t > \Delta t_{{\rm max}}$ the FP is unstable and the solutions of the time-discretized system (10) diverge although the solutions of the time-continuous system (6) converge to the FP.

If the initial value of the iteration is $\Delta R_0$, then in the linearized system the orbit of $\Delta R_0$, i.e., $\{\Delta R_0, \Delta R_1, \Delta R_2, \ldots \}$ is given by

$$% \Delta R_{{n}+1} = \Phi(\Delta R_{n}) \approx \overline{\De... ...\right)^{n+1}~ \left(\Delta R_0 - \overline{\Delta R} \right),$$ (13)

which can be shown by induction over n.

The following cases are possible (without proof):

1.

$\Delta t < \frac{H_{\rm P}}{R_2} \tau_{\rm d}'$: the orbit converges directly to the FP; for $\Delta t \rightarrow 0$ the continuous limit is reached.

2.

$\Delta t = \frac{H_{\rm P}}{R_2} \tau_{\rm d}'$: the orbit converges after one iteration to the FP.

3.

$\frac{H_{\rm P}}{R_2} \tau_{\rm d}' < \Delta t < 2 \frac{H_{\rm P}}{R_2} \tau_{\rm d}'$: since $D\Phi(\overline{\Delta R})$ becomes negative, the orbit converges alternatingly to the FP.

4.

$2 \frac{H_{\rm P}}{R_2} \tau_{\rm d}' < \Delta t$: the orbit diverges.

Therefore, any binary evolutionary code which uses an explicit iteration scheme like Eq. (1) encounters this numerical instability. This is the case even if we assume that the binary evolutionary code solves the equations of stellar structure with infinite accurracy because this instability is a result of the time-discretization itself which is basically unavoidable for numerical computations.

It is now possible to estimate how many time steps an explicit method like Eq. (1) needs at least for the computation of a mass transfer phase during which the maximum time step length  $\Delta t_{{\rm max}}$ is used. From Eqs. (5) and (12) it follows that at most the mass fraction

$$% \frac{\Delta M_{\rm max}}{M_2} = \left\vert \frac{ \overlin... ... \frac{2}{\zeta_{\rm s} - \zeta_{\rm R}} \frac{H_{\rm P}}{R_2}$$ (14)

can be removed from the donor per time step. After n time steps the initial mass M₂⁽⁰⁾ has been reduced to

$$% M_2({t_{n}}) = \left(1 - \frac{2}{\zeta_{\rm s} - \zeta_{\rm R}} \frac{H_{\rm P}}{R_2} \right)^{n} M_2({t_0}).$$ (15)

For low-mass main sequence (MS) stars, $\zeta_{\rm s} \ga -\frac{1}{3}$(Hjellming 1989; Hjellming & Webbink 1987); for sufficiently small mass ratio M₂ / M₁ of the binary, $\zeta_{\rm R} \rightarrow -\frac{5}{3}$in the analytical approximation of Paczynski (1971), and therefore $\zeta_{\rm s} - \zeta_{\rm R} \approx 1$. As an example: for $\frac{H_{\rm P}}{R_2} = 10^{-4}$ which is a typical value for low-mass MS stars and the low-mass limits $\zeta_{\rm s} = -\frac{1}{3}$and $\zeta_{\rm R} = -\frac{5}{3}$ at least 4600 time steps are necessary to reduce the donor mass by a factor of two according to Eq. (15). This is the most optimistic limit where $\Delta t = \Delta t_{\rm max}$. For thermally unstable systems, which can even approach the onset of dynamical instability Hjellming & Webbink 1987; Schenker et al. 2002; and Beer & Podsiadlowski 2002, easily 10⁴ or 10⁵ time steps are necessary to reduce the donor mass by a factor of two.

The main reason for this increased computational demand in the case of thermally unstable mass transfer can be understood as follows:

The donor star of a binary in which mass transfer is thermally unstable has a deep radiative envelope. Unperturbed stars with a deep radiative envelope can be described (to some extent) by a polytropic stellar structure with a polytropic index $n \la 3$ and an adiabatic index $\gamma = 5/3$ (monoatomic ideal gas), see e.g. Hjellming & Webbink (1987), Table 1. Because for polytropes of index n, $\zeta_s = (1-n)/(3-n)$, for radiative stars (with $n \la 3$) $\zeta_s$ is a very large, positive number. However, this holds only for unperturbed stars in thermal equilibrium and only in the limit of infinitesimally small mass loss (see Hjellming & Webbink (1987), last paragraph of Sect. II). On the other hand, for finite mass loss $\zeta_s$ decreases rapidly with the amount of mass lost. Yet for such stars, at least during the initial phases of thermal timescale mass transfer, $\zeta_s$ is still much larger than unity, i.e. typically of order 10-100 (as shown in Hjellming & Webbink 1987, their Figs. 3 and 4). As a consequence, $\zeta_s{-}\zeta_{\rm R}$ is then also of the order of 10-100, and  $\Delta M_{\rm max}$ as given in Eq. (14) is smaller and the minimum number of required time steps increases by that factor.

As can be seen from Eq. (15), an artificial increase of $H_{{\rm P}}$ can significantly reduce the required number of time steps. This has been used in the past by numerous authors to speed up the computations (e.g., Hameury 1991). This approach yields the correct mass transfer rate only when mass transfer is close to stationary. However, it cannot be used if the turn-on and turn-off of mass transfer is important, as is the case for irradiation-induced mass transfer cycles (cf. Paper I).

   
4.2 The Feigenbaum scenario

The linear stability analysis discussed in Sect. 4.1 describes only the system dynamics near the FP, i.e., the local dynamics; we will now briefly discuss the global dynamics.

When going to dimensionless quantities

$$% x := \frac{\Delta R - \overline{\Delta R}}{H_{\rm P}},\quad \delta := \frac{R_2}{H_{\rm P}} \frac{\Delta t}{\tau_{\rm d}'}$$ (16)

and using Eqs. (2) and (5), Eq. (10) is equivalent to

$$% x_{{n}+1} = x_{n} + \left[ 1 - \exp(x_{n}) \right] ~\delta =: x_{n} + G(x_{n}) =: f_\delta(x_{n}).$$ (17)

The family of functions $f_\delta$ is topologically conjugated to a family of $\cal S$-unimodal functions: thus, $f_\delta$ is qualitatively similar to the well-known logistic map:

$$% g_\delta(x) = \delta x~(1 - x)$$ (18)

and undergoes a Feigenbaum scenario for increasing $\delta$ which is called the control or chaos parameter.

For $\delta = 2$, corresponding to $\Delta t = \Delta t_{\rm max}$, the FP becomes unstable and bifurcates into an unstable FP and a stable cycle of period 2. For increasing values of $\delta$, the 2-cycle also becomes unstable, and a stable 4-cycle appears, and so on. At a critical value of $\delta \approx 2.7$, the system dynamics become chaotic, and beyond this value the chaotic dynamics are permanently interrupted by finite windows of regular dynamics which correspond to the existence of a stable cycle of any period. Figure 1 shows 50 iterations of Eq. (17) of one single orbit for 5000 different values of $\delta$. The FP $\bar{x} \equiv 0$ corresponds to  $\overline{\Delta R}$. At about $\delta \approx 3.1$, a stable orbit of period 3 appears, whose existence is a mathematical proof for the existence of chaotic dynamics in the system according to the famous "Period three implies chaos'' by Li & Yorke (1975). x is a logarithm of $\dot{M}_2$ and $\delta$ is linear in $\Delta t$. Thus, the transitional region from stable to chaotic dynamics is rather small.

	Figure 1: Feigenbaum diagram of map (17). 50 iterations of one orbit for 5000 different values of $\delta$are shown.
Open with DEXTER

Since chaotic or "random'' values in the mass transfer rate are not desirable, how can the chaotic dynamics be suppressed? As mentioned above, an artificial increase of $H_{{\rm P}}$ will increase $\Delta t_{{\rm max}}$ and therefore shift the onset of chaos to higher values of $\delta$. In general, the choice of a mass transfer prescription which differs from Eq. (2) will only shift the onset of chaos to different values of $\delta$ but it cannot suppress it. At least, if $\dot{M}_2$ grows as $x^\alpha$, $f_\delta$basically keeps its behaviour (without proof) for  $\alpha \geq 2$ and this is the case for all physical models for computing mass transfer.

   
4.3 Controlling Chaos

It is possible to suppress chaotic dynamics and to stabilize a given FP by introducing small pertubations in terms of "Controlling Chaos''. There are two distinct approaches: the first one varies one system parameter by a small amount  $\varepsilon_{n}$ for each time step to enforce stability of the FP (Ott et al. 1990), but this method requires the a priori knowledge of the FP, and this is not the case for our system. The second method adds a small correction  $\varepsilon_{n}$ to the new iteration value, i.e.,

$$% x_{{n}+1} := f(x_{n}) + \varepsilon_{n}.$$ (19)

The simplest approach uses $\varepsilon_{n} = K~(x_{n} - x_{n-1})$with a suitable constant K (Pyragas 1992). This so-called delayed dynamical feedback involves x_n-1 and x_n, i.e., it uses information from previous iteration steps. This can extend the stability limit to larger values of $\delta$ depending on K, but as numerical experience shows, this is by far not sufficient in our case. The inclusion of all prior iteration values (Socolar et al. 1994) like $\varepsilon_{n} = K~(x_{n} - x_{n-1}) + R \varepsilon_{n-1}$ with a suitable constant R achieves significantly better results. For $R\rightarrow 1$ the stability of the FP can be extended to arbitrarily large values of $\delta$ but at the expense of an arbitrarily small basin of attraction which makes this approch difficult for practical application.

A nonlinear approach of the form

$$% x_{n+1} = f(x_{n-1}) + K~(f(x_{n}) - f(x_{n-1})) + R \varepsilon_{n}$$ (20)

as has been proposed by de Sousa Vieira & Lichtenberg (1996) yields a significant improvement of the basin of attraction. For the specific choice of R=K and $K = (Df(\bar{x}) - 1)^{-1} DF(\bar{x})$, the FP becomes superstable, i.e., the FP iteration (20) converges quadratically in a sufficiently small neighbourhood of the FP.

Since our map (17) is of the form $f_\delta(x) = x + G_\delta(x)$, Eq. (20) is equivalent to

$$% x_{n+1} = x_{n} + DG_\delta^{-1}(\bar{x}) G_\delta(x_{n}).$$ (21)

Because the position of the FP is a priori unknown, the best available approximation to  $DG_\delta(\bar{x})$ is given by  $DG_\delta(x_{n})$. Then Eq. (21) turns into a standard Newton's method for  $G_\delta(x)$. Therefore, we conclude that Newton's method is most likely the best available method to compute $\bar{x}$, and  $\overline{\dot{M}}_2$ with a reasonable computational effort.

   
5 The implicit algorithm

While the explicit iteration scheme is a simple time integration scheme which requires one iteration per time step, the proposed implicit iteration scheme is a fixed point search which requires several iterations per time step until the fixed point is found with sufficiently high accuracy.

Fortunately, the equations of stellar structure are typically solved by Newton's method, more explicitly by the so-called Henyey method (Kippenhahn et al. 1967; Kippenhahn & Weigert 1990). Thus, the best solution is to include the mass transfer rate or the stellar mass itself as an additional variable in the Henyey method and to perform the fixed point search simultaneously with the solution of the equations of stellar structure (Paper I). This has also been done independently by Benvenuto & de Vito (2003). For details of our implementation, see Büning (2003).

Figure 2 shows the mass transfer rate obtained with our binary evolutionary code as a function of time for one specific binary system. Since our evolutionary code has not been designed for the analysis of chaotic dynamics, we have kept $\Delta t$ constant and instead have used the fact that the timescale of thermal relaxation varies, which is the dominant term in  $\tau _{\rm d}'$ for the binary system in question. When mass transfer starts, the thermal relaxation of the donor increases (i.e., $\tau _{\rm d}'$ decreases, cf. Eq. (16)), and, therefore, $\delta$ also increases. At some point, thermal relaxation reaches a maximum and decreases afterwards, and so does $\delta$.

Figure 2: Mass transfer rate $\dot{M}_2 [M_\odot/ {\rm yr}]$ as a function of time $t[{\rm yr}]$ for an early case A mass transfer in a binary system with $M_1 = M_2 = 2~M_\odot$ and using a constant time step length of $5000~{\rm yr}$. The variation of the chaos parameter $\delta$ is caused by the variation of thermal relaxation. The system dynamics are most unstable when $\tau _{\rm d}'$ is minimal. Every dot corresponds to one single time step in the calculation with the explicit scheme (1); the dashed line shows the result of the same computation, but using the implicit method.

Open with DEXTER

When using the explicit algorithm (1), the resulting mass transfer rate undergoes a period doubling bifurcation before it exhibits chaotic dynamics. Within the chaotic region a periodic window of period 5 appears. After thermal relaxation and hence also $\delta$ has reached its maximum, the Feigenbaum scenario evolves backwards. Since the decrease of $\delta$ occurs slower than its increase before, even stable cycles of period 8 and 4 appear before the secular mass transfer rate finally becomes stable. In contrast, the result obtained with our implicit algorithm for the same system and the same system parameters is shown by the dashed line.

Although the implicit algorithm stabilizes the FP and even provides a quadratic convergence of the iteration, chaotic dynamics are still present outside of a neighbourhood of the FP. Therefore, a good initial value for the iteration is necessary. We used a linear extrapolation of $\dot{M}_0$ and $H_{{\rm P}}$ to determine the initial value for $\dot{M}_2$ by Eq. (2). First, we perform typically 2-4 iterations and keep the preestimated value for $\dot{M}_2$constant until the solution for the other stellar parameters has almost converged. Otherwise, even small fluctuations of R₂ might push the next iteration value of $\dot{M}_2$ out of the basin of attraction of the FP and prevent convergence. Then, we finish with mostly 2-6 iterations using the full implicit algorithm to determine the correct mass transfer rate.

Furthermore, to calculate the mass transfer rate with a numerical accuracy of better than $1\%$, the stellar radius has to be determined with a very high numerical accurracy. According to Eq. (2), for a scale height of $H_{\rm P} \approx 10^{-4}~ R_2$, which is typical for low-mass MS stars, R₂ has to be determined with a relative accuracy of the order of 10^-6 in order to get $\dot{M}_2$ with a relative accuracy of the order of 10^-2. To compute the radius with such a high accuracy the stellar physics and especially the equation of state must be a very smooth function of its variables. Every small discontinuity, especially in the outer layers of the star, can cause small jumps in the stellar radius which appear, magnified by the factor  $R_2/H_{\rm P}$, as significant jumps in the mass transfer rate.

   
6 Summary and conclusions

We have used a simple analytical model, a 1d autonomous ordinary differential equation to describe the time evolution of the mass transfer rate. We have shown that, while the FP in the time-continuous system, i.e., the "physical'' mass transfer rate is stable, the FP in the time-discretized system, i.e., the "numerical'' mass transfer rate becomes unstable if the length of the time step $\Delta t$exceeds a critical value  $\Delta t_{{\rm max}}$ given by Eq. (12). We have estimated that even in the ideal case where $\Delta t = \Delta t_{\rm max}$ at least several thousand time steps are necessary to reduce the donor mass in a low-mass binary system by a factor of two. For systems with thermally unstable mass transfer, it is even worse.

We outline a mathematical proof that the iteration equation for the time-discretized system shows a behaviour similar to that of the logistic map and, for $\Delta t > \Delta t_{{\rm max}}$, undergoes a series of period doublings which finally leads to chaotic dynamics, i.e., to apparently random values of the computed mass transfer rate. The choice of a different explicit prescription to calculate the mass transfer rate results only in a shift of the critical time step length  $\Delta t_{{\rm max}}$ which, in turn, depends on the characteristic scale length H at the FP  $\overline{\Delta R}$, i.e., at the secular mass transfer rate.

In terms of "Controlling Chaos'' we have briefly discussed several methods to stabilize the FP. Various modified iteration schemes which are equivalent to different explicit iteration schemes have been discussed in the literature, but they all do not show sufficient stabilization. Therefore, we suggest that using a different explicit iteration scheme may shift  $\Delta t_{{\rm max}}$ to higher values but will not solve the problem. An implicit scheme, Newton's method, is the most promising solution because it stabilizes the FP formally for $\Delta t \rightarrow \infty$, although for this to be the case a sufficiently good initial value for the iteration is required.

In practice, the implicit algorithm reduces the number of required time steps by at least a factor of 10. Another advantage is that the implicit algorithm either yields the "correct'' mass transfer rate or does not converge at all, whereas the explicit algorithm provides a result in every case, even if it is random. Therefore, in our binary evolutionary calculations we reject results of the last time step if convergence is not reached and restart it with a smaller $\Delta t$.

Acknowledgements

We thank A. Weiss and H. Schlattl for providing their stellar evolutionary code, and H. Schlattl for valuable support and helpful discussions about numerics. We also thank U. Kolb for providing unpublished details about his numerical calculations.

References

• Arnold, V. I. 1983, in Grundlehren der mathematischen Wissenschaften: A Series of Comprehensive Studies in Mathematics, 250, Geometrical Methods in the Theory of Ordinary Differential Equations, ed. M. Berger, B. Eckmann, & S. R. S. Varadhan (New York, Berlin: Springer-Verlag)
• Beer, M. E., & Podsiadlowski, P. 2002, MNRAS, 335, 358 [NASA ADS] [CrossRef] (In the text)
• Benvenuto, O. G., & de Vito, M. A. 2003, MNRAS, 342, 50 [NASA ADS] [CrossRef] (In the text)
• Boccaletti, S., Grebogi, C., Lai, Y.-C., Mancini, H., & Maza, D. 2000, Phys. Rep., 329, 103 [NASA ADS] [CrossRef]
• Büning, A. 2003, Ph.D., Ludwig-Maximilians-Universität München (In the text)
• Büning, A., & Ritter, H. 2004, A&A, 423, 281 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Collet, P., & Eckmann, J.-P. 1980, in Progress in Physics, 1, Iterated Maps on the Interval as Dynamical Systems, ed. A. Jaffe, & D. Ruelle (Boston: Birkhäuser Verlag)
• Cvitanovic, P., Ed. 1984, Universality in Chaos (Bristol: Adam Hilger Ltd) (In the text)
• D'Antona, F. 1994, Mem. Soc. Astron. It., 65, 285 [NASA ADS] (In the text)
• D'Antona, F., Mazzitelli, I., & Ritter, H. 1989, A&A, 225, 391 [NASA ADS] (In the text)
• de Sousa Vieira, M., & Lichtenberg, A. J. 1996, Phys. Rev. E, 54, 1200 [NASA ADS] [CrossRef] (In the text)
• Feigenbaum, M. J. 1978, J. Stat. Phys., 19, 25 [CrossRef]
• Feigenbaum, M. J. 1980, Los Alamos Science, 1, 4, see also reprint in Cvitanovic (1984)
• Guckenheimer, J., & Holmes, P. 1983, in Applied Mathematical Sciences, 42, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (New York: Springer-Verlag) (In the text)
• Hameury, J.-M. 1991, A&A, 243, 419 [NASA ADS] (In the text)
• Hjellming, M. S. 1989, Ph.D., Illinois University at Urbana-Champaign
• Hjellming, M. S., & Webbink, R. F. 1987, ApJ, 318, 794 [NASA ADS] [CrossRef]
• Jackson, E. A. 1989, Perspectives of Nonlinear Dynamics (Cambridge: Cambridge University Press)
• Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution (Berlin: Springer-Verlag)
• Kippenhahn, R., Weigert, A., & Hofmeister, E. 1967, in Methods in Computational Physics, Astrophysics, ed. B. Alder, S. Fernbach, & M. Rosenberg (New York: Academic Press), 7, 129
• Kolb, U., Davies, M. B., King, A., & Ritter, H. 2000, MNRAS, 317, 438 [NASA ADS] [CrossRef] (In the text)
• Kolb, U., & Ritter, H. 1990, A&A, 236, 385 [NASA ADS] (In the text)
• Li, T.-Y., & Yorke, J. A. 1975, Am. Math. Monthly, 82, 985 (In the text)
• Ott, E., Grebogi, C., & Yorke, J. A. 1990, Phys. Rev. Lett., 64, 1196 [NASA ADS] [CrossRef] (In the text)
• Paczynski, B. 1971, ARA&A, 9, 183 [NASA ADS] (In the text)
• Pyragas, K. 1992, Phys. Lett. A, 170, 421 [NASA ADS] [CrossRef] (In the text)
• Ritter, H. 1988, A&A, 202, 93 [NASA ADS] (In the text)
• Sarna, M. J. 1992, MNRAS, 259, 17 [NASA ADS] (In the text)
• Schenker, K., King, A. R., Kolb, U., Wynn, G. A., & Zhang, Z. 2002, MNRAS, 337, 1105 [NASA ADS] [CrossRef] (In the text)
• Socolar, J. E., Sukow, D. W., & Gauthier, D. J. 1994, Phys. Rev. E, 50, 3245 [NASA ADS] [CrossRef] (In the text)
• Thompson, J. M. T., & Stewart, H. B. 1986, Nonlinear Dynamics and Chaos (Chichester, New York: John Wiley and Sons) (In the text)
• Tout, C. A., Aarseth, S. J., Pols, O. R., & Eggleton, P. P. 1997, MNRAS, 291, 732 [NASA ADS] (In the text)
• Wiggins, S. 1990, in Texts in Applied Mathematics, 2, Introduction to Applied Nonlinear Dynamical Systems and Chaos (New York: Springer-Verlag)