# Decoding of the light changes in eclipsing Wolf-Rayet binaries. I. A non-classical approach to the solution of light curves

C. Perrier, J. Breysacher, G. Rauw

► To cite this version:

C. Perrier, J. Breysacher, G. Rauw. Decoding of the light changes in eclipsing Wolf-Rayet binaries. I. A non-classical approach to the solution of light curves. Astronomy & Astrophysics - A&A, 2009, 503, pp.963. 10.1051/0004-6361/200911707. hal-00479926

HAL Id: hal-00479926

https://hal.science/hal-00479926v1

Submitted on 28 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Decoding of the light changes in eclipsing Wolf-Rayet binaries

## I. A non-classical approach to the solution of light curves

C. Perrier<fnref n="1" />, J. Breysacher<fnref n="2" />, and G. Rauw<fnref n="3" />,\*

Received 23 January 2009 / Accepted 2 June 2009

## ABSTRACT

**Aims.** We present a technique to determine the orbital and physical parameters of eclipsing eccentric Wolf-Rayet + O-star binaries, where one eclipse is produced by the absorption of the O-star light by the stellar wind of the W-R star.

**Methods.** Our method is based on the use of the empirical moments of the light curve that are integral transforms evaluated from the observed light curves. The optical depth along the line of sight and the limb darkening of the W-R star are modelled by simple mathematical functions, and we derive analytical expressions for the moments of the light curve as a function of the orbital parameters and the key parameters of the transparency and limb-darkening functions. These analytical expressions are then inverted in order to derive the values of the orbital inclination, the stellar radii, the fractional luminosities, and the parameters of the wind transparency and limb-darkening laws.

**Results.** The method is applied to the SMC W-R eclipsing binary HD 5980, a remarkable object that underwent an LBV-like event in August 1994. The analysis refers to the pre-outburst observational data. A synthetic light curve based on the elements derived for the system allows a quality assessment of the results obtained.

**Key words.** stars: early-type – stars: mass-loss – binaries: eclipsing – stars: individual: HD 5980

## 1. Introduction

Photometric monitoring of Wolf-Rayet (W-R) binaries revealed that many of them display a shallow eclipse when the W-R star passes in front of its O-type companion (e.g. Lamontagne et al. 1996). These so-called atmospheric eclipses arise when part of the light of the O-type companion is absorbed by the wind of the W-R star. In a few cases, the light curve displays an eclipse at both conjunctions and the analysis of this phenomenon can provide important information about the physical parameters of W-R stars and their winds. In this context, the most famous example is V444 Cyg (WN5 + O6 V), which has been extensively investigated by the Moscow group (e.g. Antokhin & Cherepashchuk 2001 and references therein). Cherepashchuk and coworkers (e.g. Cherepashchuk 1975; Antokhin & Cherepashchuk 2001 and references therein) developed a sophisticated method to handle the ill-posed problem of light curve inversion for V444 Cyg. Based on the minimum a priori assumptions about the transparency function, this method not only yields the radii of both components, but also provides information about the structure of the WN5's stellar wind. However, because of a number of fundamental hypotheses that are not necessarily valid for all eclipsing W-R binaries, this method cannot be readily applied to all eclipsing W-R + O systems. For eccentric systems in particular (e.g. WR22, Gosset et al. 1991) some assumptions (such as spherical symmetry of the problem) break down and a different technique must be used.

We initiated our study of the observed light changes of eclipsing Wolf-Rayet binary systems when two of us (J.B. and C.P.) tried to confirm the 25.56 day period found by Hoffmann et al. (1978) for the SMC star HD 5980, the first extragalactic Wolf-Rayet binary then known to display eclipses. This exercise led to the discovery of the correct orbital period of HD 5980, $P = 19.266 \pm 0.003$ days (Breysacher & Perrier 1980), the light curve revealing, in addition, a rather eccentric orbit of $e = 0.47$ assuming $i = 80^\circ$. However, because of the uncertainties in the depth of both minima, caused by an insufficient number of observations, no detailed quantitative analysis of this preliminary light curve could be attempted.

Its relatively long period and large eccentricity ensure that HD 5980 is an interesting object in which to study the structure of a W-R envelope, and a photometric monitoring of this system was initiated to define the shape of its light curve in a more accurate way. More than 700 observations were collected. After realizing that none of the existing “classical tools” was suited to our purpose – the decoding of the light changes of a *partially-eclipsing system* characterized by an *eccentric orbit* and containing one component with an *extended atmosphere* – we started to develop another approach to the solution of light curves.

The technique of light curve analysis applied to V444 Cyg by Smith & Theokas (1980), which is based on Kopal's fundamental work (cf. Kopal 1975, 1979), appeared as an attractive approach to the solution of our problem. This method is based on the interpretation of the observed light changes in the *frequency-domain*, i.e., not the light curve as a function of time, but its *Fourier-like integral transform*.

<footnote n="1">¹ Observatoire de Grenoble, Université Joseph Fourier, 38041 Grenoble Cedex, France</footnote>

<footnote n="2">² Lieu-dit Petit Riston, 40230 Saubion, France</footnote>

<footnote n="3">³ Institut d'Astrophysique et de Géophysique, Université de Liège, Allée du 6 Août, Bât B5c, 4000 Liège (Sart Tilman), Belgium
e-mail: rauw@astro.ulg.ac.be</footnote>

<footnote>\* Also Research Associate FNRS (Belgium).</footnote>
We now describe in detail the method we developed for the study of Wolf-Rayet eclipsing binaries. A preliminary application of this technique to the light curve of HD 5980 prior to its 1994 LBV-like eruption (see e.g., Bateson & Jones 1994, Barbá et al. 1995, Heydari-Malayeri et al. 1997) was presented by Breysacher & Perrier (1991, hereafter BP91). We reanalyse the pre-outburst observational data using an improved version of the software tool. Revised values for the physical parameters of HD 5980 are derived. A synthetic light curve based on the elements thus obtained allows a quality assessment of the new results.

A more detailed study of HD 5980, including the analysis of the light curve obtained after the eruption (Sterken & Breysacher 1997), will be presented in a forthcoming paper.

## 2. Analysis of the light changes in the frequency-domain

In this section, we present the fundamental equations of our method, we then introduce the mathematical functions to model the transparency and limb-darkening functions and consider the specific problem of eccentric orbits.

### 2.1. The basic equations

We refer to the fundamental work of the Manchester group (cf. Kopal 1975, 1979; Smith 1976), and first consider an eclipsing system consisting of two spherical stars revolving around the common centre of gravity in circular orbits, and appearing in projection on the sky as uniformly bright discs. The system is seen at an inclination angle $i$. When star 1 of fractional luminosity $L_1$ and radius $r_1$ is partly eclipsed by star 2 of fractional luminosity $L_2$ and radius $r_2$ (Fig. 1), the apparent brightness $l$ of the system (maximum light between minima taken as unity) is given by

$$l(r_1, r_2, \delta, J) = 1 - \int \int_A J(r) \, d\sigma, \quad (1)$$

where $\delta$ is the apparent separation of the centres of the two discs, $J$ represents the distribution of brightness over the apparent disc of the star undergoing eclipse, and $d\sigma$ stands for the surface element. The distances $r_1, r_2$, and $\delta$ are expressed in units of the orbital separation. The integral in Eq. (1) provides the apparent “loss of light” displayed by the system when an area $A(r_1, r_2, \delta)$ of star 1 is eclipsed (see Fig. 1). The assumption that star 1 is uniformly bright yields

$$J(r) = \frac{L_1}{\pi r_1^2}. \quad (2)$$

Combining Eqs. (1) and (2), we obtain

$$1 - l(r_1, r_2, \delta, J) = \frac{L_1}{\pi r_1^2} \int \int_A d\sigma = \alpha L_1, \quad (3)$$

where $\alpha$ is the ratio of the mutual area of eclipse to the area of the disc of the eclipsed star, and is a function of $r_1, r_2$, and $\delta$ (Kopal 1975). A generalisation of these concepts to the case of spherical stars with arbitrary limb-darkening laws $J(r)$ was presented by Smith (1976).

For an orbital period $P$ and an epoch of conjunction $t_0$, we define the phase angle $\theta$ at a time $t$ to be

$$\theta = \frac{2\pi}{P}(t - t_0). \quad (4)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Geometry of the eclipse of star 1 by star 2. The integral of Eq. (1) is evaluated over the hatched area $A(r_1, r_2, \delta)$. The points inside this area can be specified by the coordinates $(r, s)$ corresponding to the two intersections of the circle with radius $r$ and centre $O_1$ and the circle with radius $s$ centred on $O_2$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Light curve of an eclipse in the $(l, \sin^2m\theta)$ plane. $\theta_{fc}$ corresponds to the phase angle of first contact, whilst $\theta_t$ represents the phase angle corresponding to the beginning of the total eclipse. The shaded area illustrates the moment $A_{2m}$ and $\lambda = 1 - l(\theta = 0)$.</figcaption></figure>

As proposed by Kopal (1975, 1979), we focus our attention on the area subtended by the light curve in the $(l, \sin^2m\theta)$ plane, where $m$ is a positive integer number ($m = 1, 2, 3, ...$), as shown in Fig. 2. The areas $A_{2m}$ between the lines $l = 1, \sin^2m\theta = 0$, and the true light curve are then given by the integrals

$$A_{2m} = \int_{0}^{\theta_{fc}} (1 - l) d(\sin^2m\theta), \quad (5)$$

which are hereafter referred to as the *moments of the eclipse*, of index $m$, where $\theta_{fc}$ denotes the phase angle of the first contact ($\delta(\theta_{fc}) = r_1 + r_2$) of the eclipse.

Combining Eqs. (1) and (5), Kopal (1975, 1979) and Smith (1976) demonstrated that, based on certain assumptions, it is possible to:

* derive analytical expressions of the moments of the eclipse in terms of the physical parameters $(i, r_1, r_2, L_1, L_2)$ of a binary system with a circular orbit consisting of uniformly bright spherical stars;
* invert these relationships to determine the parameters of the system in terms of the moments $A_{2m}$ that can be empirically obtained from the data;

* extend this treatment to the case of partial eclipses of stars with an arbitrary (yet analytical) limb-darkening law.

Going one step further, Smith & Theokas (1980) generalised the above concepts to derive a convenient mathematical solution to the problem of an atmospheric eclipse, i.e., an eclipse of a limb-darkened star by a star surrounded by an extended atmosphere. In such an eclipse, at each position inside the area A specified by the coordinates r and s (see Fig. 1), a fraction of the light emitted by star 1 is absorbed by the atmosphere of star 2. To account for these transparency effects, a transparency function $F(s)$ is introduced into Eq. (1), so that the total amount of light seen by the observer becomes

$$l(r_1, r_2, \delta, J, F) = 1 - \int_A J(r)F(s)ds. \quad (6)$$

Considering that in the case of an atmospheric eclipse, it might be interesting to give more weight to the data close to mid-minimum, Smith & Theokas (1980) also introduced an alternative set of moments $B_{2m}$ defined by

$$B_{2m} = - \int_{0}^{\theta_{fc}} (1-l) d(\cos^{2m}\theta). \quad (7)$$

The use of kernel $d(\cos^{2m}\theta)$ in Eq. (7) places more emphasis on the data points close to mid-eclipse, which leads to smaller errors than in the case of the $A_{2m}$ moments defined by means of the $d(\sin^{2m}\theta)$ kernel (Theokas & Smith 1983, private communication). We define $\epsilon$ to be the mean error in an individual data point; the relative error in the light curve $\epsilon/(1-l(\theta))$ increases for data points near $\theta_{fc}$. While these points are given more weight by the kernel of the $A_{2m}$ moments, the converse situation holds for the $B_{2m}$ moments, where the kernel reaches its peak for a given m in a zone where $1-l(\theta)$ is closer to its maximum value.

By definition, the $A_{2m}$ and $B_{2m}$ moments are related to each other by means of

$$B_{2m} = \sum_{p=1}^{m} \frac{m!}{(m-p)! p!} (-1)^{p+1} A_{2p}. \quad (8)$$

We now concentrate on the $B_{2m}$ moments because they are equally well suited to the analysis of the primary and secondary minima with the transparency and limb-darkening functions adopted in the present study (cf. Sect. 2.2). Since there are a number of typos in the paper of Smith & Theokas (1980), we provide below the mathematical details of the method.

The infinitesimal element of area $d\sigma$ of Eq. (6) is expressed as

$$d\sigma = \frac{1}{2} \frac{\partial^2}{\partial r \partial s} (\pi r^2 \alpha(r, s, \delta)) dr ds. \quad (9)$$

The following general expression for the $B_{2m}$ moments was derived by Smith & Theokas (1980)

$$B_{2m} = \lambda + \int_{0}^{r_1} \int_{0}^{r_2} J(r) F(s) ds dr \\ \times \frac{\partial^2}{\partial r \partial s} \left( \int_{0}^{\theta_{fc}} \cos^{2m} \theta \frac{\partial(\pi r^2 \alpha(r, s, \delta))}{\partial \theta} d\theta \right), \quad (10)$$

where $1 - l(\theta = 0)$ is defined as $\lambda$ (see Fig. 2).

The analytical expressions obtained for the $B_{2m}$'s for m = 1, 2, 3, 4, and 5 are thus

$$B_2 = \lambda - \csc^2 i (P - I_1 R_1 r_2^2 - \psi_1), \quad (11)$$

$$B_4 = \lambda - \csc^4 i (P - 2 I_1 R_1 r_2^2 + I_2 R_1 r_1^2 r_2^2 + I_1 R_2 r_2^4 - \psi_2), \quad (12)$$

$$B_6 = \lambda - \csc^6 i (P - 3 I_1 R_1 r_2^2 + 3 I_2 R_1 r_1^2 r_2^2 + 3 I_1 R_2 r_2^4 \\ - I_3 R_1 r_1^4 r_2^2 - I_1 R_3 r_2^6 - 3 I_2 R_2 r_1^2 r_2^4 - \psi_3), \quad (13)$$

$$B_8 = \lambda - \csc^8 i (P - 4 I_1 R_1 r_2^2 + 6 I_2 R_1 r_1^2 r_2^2 + 6 I_1 R_2 r_2^4 \\ - 4 I_3 R_1 r_1^4 r_2^2 - 4 I_1 R_3 r_2^6 - 12 I_2 R_2 r_1^2 r_2^4 + I_4 R_1 r_1^6 r_2^2 \\ + I_1 R_4 r_2^8 + 6 I_3 R_2 r_1^4 r_2^4 + 6 I_2 R_3 r_1^2 r_2^6 - \psi_4), \quad (14)$$

and

$$B_{10} = \lambda - \csc^{10} i (P - 5 I_1 R_1 r_2^2 + 10 I_2 R_1 r_1^2 r_2^2 \\ + 10 I_1 R_2 r_2^4 - 10 I_3 R_1 r_1^4 r_2^2 - 10 I_1 R_3 r_2^6 \\ - 30 I_2 R_2 r_1^2 r_2^4 + 5 I_4 R_1 r_1^6 r_2^2 + 30 I_3 R_2 r_1^4 r_2^4 \\ + 30 I_2 R_3 r_1^2 r_2^6 + 5 I_1 R_4 r_2^8 - I_5 R_1 r_1^8 r_2^2 \\ - 10 I_4 R_2 r_1^6 r_2^4 - 20 I_3 R_3 r_1^4 r_2^6 - 10 I_2 R_4 r_1^2 r_2^8 \\ - I_1 R_5 r_2^{10} - \psi_5). \quad (15)$$

The coefficients $P$, $I_m$, $R_m$, and $\psi_m$ are defined (cf. Smith & Theokas 1980) by the equations

$$P(r_1, r_2, J, F) = \int_{0}^{\min(r_1, r_2)} J(r)F(r)2\pi r dr, \quad (16)$$

$$I_m(r_1, J) = \int_0^{r_1} \frac{J(r)}{r^{2m-2}} \frac{\partial}{\partial r} (\pi r^{2m}) dr, \quad (17)$$

$$R_m(r_2, F) = \int_0^{r_2} \frac{F(s)}{r_2^{2m}} \frac{\partial}{\partial s} (s^{2m}) ds, \quad (18)$$

and

$$\psi_m = \psi_m(r_1, r_2, i, J, F) = \frac{L_1}{\pi r_1^2} \left( \frac{1 + 2m(r_1 + r_2)^2}{6(r_1 + r_2)^2} \right. \\ \times (\cos^2 i - (r_2 - r_1)^2)^{3/2} - \sqrt{\cos^2 i - (r_2 - r_1)^2} \\ - \frac{m|r_2 - r_1|}{8(r_1 + r_2)^2} (\cos^4 i - (r_2 - r_1)^4) \\ \left. + |r_2 - r_1| \arctan \left( \frac{\sqrt{\cos^2 i - (r_2 - r_1)^2}}{|r_2 - r_1|} \right) \right). \quad (19)$$

## 2.2. The transparency and limb-darkening functions

In our method, the transparency of the W-R wind is described by an analytical function that depends on a limited number of parameters. Since the functional form of the transparency is adopted a priori, our choice will obviously have a direct influence on the parameters derived for the system. Therefore, it is important to clearly specify the assumptions made in our approach. To avoid confusion with the standard symbols employed by Smith and Theokas (1980), from now on, the various radii in our model will be denoted $\rho_i$ ($i = 1, 2, 3$), where $\rho_1$ stands for the radius of the O-star that undergoes the eclipse. The use of
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematical view of the transparency law across the disc of the W-R component as given by Eq. (21). The solid line shows the amount of light absorbed along the line of sight of impact parameter s. The individual contributions due to the opaque core and the semi-transparent extended atmosphere are illustrated.</figcaption></figure>

the subscripts e and a will indicate whether the components are seen in emission or absorption.

The first simplifying hypothesis is that all composite parts of the system are supposed to be spherically symmetrical. This means that the method is not applicable to close binary systems in which the components depart strongly from a spherical form as a result of tidal distortion and where ellipticity and reflection effects are both present.

Because a principle objective of the proposed technique of light curve analysis is the determination of the structure of a W-R envelope, a composite model consisting of an opaque core and a surrounding extended atmosphere was adopted for the W-R component. As a consequence, while for the transparency function F(s) of the eclipsing W-R star, of radius $r_0$, Smith & Theokas (1980) simply adopted

$$F(s) = F_{y_a}(r_0, v) = y_a \left[ 1 - v \left( \frac{s}{r_0} \right)^2 \right] \quad \text{for } s < r_0, \quad (20)$$

where v is the coefficient of transparency, we chose as a first step a transparency law of the form

$$F(s) = F_{1-y_a}(\rho_{3a}, 0) + F_{y_a}(\rho_{2a}, v), \quad (21)$$

where the radius of the opaque core of the W-R star is $\rho_{3a}$ and that of the extended eclipsing envelope is $\rho_{2a}$ ($\ge \rho_{3a}$, see Fig. 3). This transparency law was used in the preliminary analysis of the light curve of HD 5980 by BP91 and is motivated by the fact that it corresponds to a physically more realistic model of the W-R star than that defined in Eq. (20), although the advantage of relative mathematical simplicity is still preserved.

For the brightness distribution $J(r)$ over the W-R disc (important for the eclipse of the W-R component by the O-star), a law very similar to that of the transparency function is adopted

$$J(r) = J(0) \left[ J_{1-y_e}(\rho_{3e}, 0) + J_{y_e}(\rho_{2e}, u_2) \right], \quad (22)$$

where $J(0)$ is the central surface brightness, $u_2$ is the coefficient of limb-darkening, and $J_{y_e}$ is defined as

$$J_{y_e}(r_0, u) = y_e \left[ 1 - u \left( \frac{r}{r_0} \right)^2 \right]. \quad (23)$$

When the W-R star is eclipsed, the radius of the core, assumed to be of uniform brightness, becomes $\rho_{3e}$ and that of the limb-darkened envelope $\rho_{2e}$. We note that the core and envelope radii of the W-R component seen in emission or absorption may differ.

We briefly address the physical meaning of this transparency law. First, in our model, the semi-transparent envelope of the W-R star has a finite extension given by the radii $\rho_{2a}$ and $\rho_{2e}$. However, the stellar winds of W-R stars do not stop abruptly so close to the star, but instead extend to large distances (much larger than the orbital separation in most binaries) to the shock with the interstellar (or circumstellar) medium. How can we then interpret $\rho_{2a}$ and $\rho_{2e}$? The radius $\rho_{2a}$ corresponds to the farthest position in the stellar wind where the residual optical depth along the line of sight produces a variation in the light curve that can be distinguished against the intrinsic photometric variability of the W-R star and the photometric errors. Similarly $\rho_{2e}$ is the outer radius of the W-R envelope that emits a measurable fraction of the light in the considered waveband.

A clear difference between our approach and that of Antokhin & Cherepashchuk (2001) concerns the functional behaviour of $F(s)$: whereas in our model, $F(s)$ is a convex function over the entire range $\rho \in [0, \rho_{3a}]$, Antokhin & Cherepashchuk (2001) use a convexo-concave function, where the concave part corresponds to the stellar wind. As a consequence, $F(s)$ given by Eq. (21) decreases at a slower rate over the wind than the transparency law inferred by Antokhin & Cherepashchuk.

In contrast to Smith and Theokas (1980) who neglect the effect, we take into account the limb-darkening of the OB-type component. Assuming that the formula employed by these authors to represent the brightness distribution across the W-R disc also applies to normal stars, we adopt the following limb-darkening law for the OB star

$$J(r) = \frac{L_O}{\pi \rho_1^2 (1 - u_1 + u_1^2/3)} \left[ 1 - u_1 \left( \frac{r}{\rho_1} \right)^2 \right]^2, \quad (24)$$

where $L_O$ is the luminosity of the OB star, of radius $\rho_1$, and $u_1$ is the coefficient of limb-darkening at the effective wavelength of the photometric filter considered.

Using these laws of transparency and limb-darkening, we then derived the expressions for $P$, $I_m$, $R_m$ (see Appendix A), and $\psi_m$, and hence the final equations for the moments $B_{2m}$, corresponding to the primary and secondary minima.

## 2.3. The orbital eccentricity

The treatment of elliptical orbits in the frequency-domain was also addressed by Kopal (1979). The problem still concerns the determination of the elements of the eclipse from the moments $B_{2m}$ in the present case – derived from the light curve, but accounting for the eccentricity $e$ and the longitude of periastron $\omega$.

In the definition of the $B_{2m}$ moments, the phase-angle $\theta$ is no longer identical to the mean anomaly $M$ but has rather to be replaced by a linear function of the true anomaly $v$

$$\theta = v + \omega - \frac{\pi}{2}. \quad (25)$$
The $d(\cos^{2m} \theta)$ kernel in Eq. (7) thus becomes

$$d(\cos^{2m} \theta) = d[\sin^{2m}(v + \omega)]. \qquad (26)$$

As a consequence, the empirical values of $B_{2m}$ cannot be derived from the observed data until a proper conversion of the phase angle into the true anomalies has been completed. This can be achieved either by a numerical inversion of Kepler's equation or the well-known asymptotic expansion of elliptical motion (e.g. Danjon 1959; Kopal & Al-Naimiy 1978)

$$v = M + \left(2e - \frac{1}{4}e^3\right) \sin M + \left(\frac{5}{4}e^2 - \frac{11}{24}e^4\right) \sin 2M + \frac{13}{2}e^3 \sin 3M + \frac{103}{96}e^4 \sin 4M + \dots \qquad (27)$$

Regardless of the technique used to compute the true anomaly, this conversion evidently requires an a priori knowledge of $e$ as well as $\omega$. For a given value of the inclination $i$, these parameters can be derived by inversion of the equations (see e.g. Kopal & Al-Naimiy 1978)

$$\Delta\Phi = \frac{1}{2} + \frac{e \cos \omega}{\pi} \left\{ 1 + \csc^2 i - \frac{e^2}{2} \left[ \frac{8}{3} \cos^2 \omega - 2 + O(\cot^2 i) \right] \right\}, \qquad (28)$$

and

$$e \sin \omega = \frac{d_2 - d_1}{4 \sin \left(\frac{d_1+d_2}{4}\right)} \left[1 - \frac{\cot^2 i}{\sin^2 \left(\frac{d_1+d_2}{4}\right)}\right]^{-1}, \qquad (29)$$

where $\Delta\Phi$, $d_1$, and $d_2$ are, respectively, the phase displacement of the minima and the durations of the primary and secondary stellar core eclipses. These latter quantities are determined directly from the observed light curve. Since the orbital inclination of an eclipsing binary system is likely to be rather large, the $O(\cot^2 i)$ term in the coefficient of the $e^3$ term of Eq. (28) can be neglected.

The empirical “elliptical” moments of the light curve then provide the elements of the binary exactly as in the “circular” case. One must ensure that the resulting values of the radii have been reduced to a constant unit of length. In our code, we therefore report all distances in relation to the semi-major axis $a$ of the relative orbit.

Our composite model adopted for the W-R star has different radii ($\rho_{2,3a}$ and $\rho_{2,3c}$) depending on whether the W-R component is seen as an eclipsing or an eclipsed disc. Because of this, although each individual half-eclipse provides an independent solution for the elements, the complete determination of the elements requires a combination of solutions obtained for both minima and because of the non-zero eccentricity, both the descending and ascending branches of each.

# 3. Decoding of the light curve

## 3.1. Empirical determination of the moments

We consider the light curve of a W-R binary system of eccentricity $e$ and period $P$, derived for a given photometric bandpass. The data are assumed to consist of a list of entries that provide for each observation the orbital phase $\Phi_i$ and the measured intensity $l_i$. To normalize the brightness scale, a mean intensity value is derived well outside the eclipses, during a phase-interval where the system is assumed to display (constant) maximum light. The $l_i$ value of each data point is then divided by this mean

to normalize the light curve to unity. It has to be emphasized, however, that this does not necessarily imply that $L_1 + L_2 = 1$ for the W-R binary. The luminosity of a third photometrically unresolved component along the line of sight may indeed contribute to the observed brightness as well, thereby leading to a brightness distribution such as $L_1 + L_2 < 1$ in the final solution for the eclipse. In the case of an eccentric orbit, the eclipse-free mean $l_i$ value is preferably taken around apastron to avoid as much as possible any luminosity increase that could occur around periastron as a result of enhanced interaction effects between the components.

The determination of the moments $B_{2m}$ requires a smoothed light curve, which can be obtained, for instance, from a spline fit to the observed points with special attention to the minima. However, this task can become difficult if the descending or ascending branch of either minimum is ill-defined because of an uneven sampling of the observations or intrinsic photometric variability in the W-R star (see e.g. the case of WR 22, Gosset et al. 1991). A clustering of the points, in particular, is a serious handicap for the method. A preliminary processing performed, by filtering the observational data, to help reduce the scatter allows us then to obtain a smooth light curve $l(\Phi)$.

The next step consists of determining the quantities $d_1$, $d_2$, and $\Delta\Phi$ (see above) from the smoothed light curve. For an assumed value of the orbital inclination $i$, the parameters $e$ and $\omega$ are obtained by means of an inversion of Eqs. (28) and (29). With these values of $e$ and $\omega$, the orbital phases $\Phi_i$ are converted into true anomalies.

The moments $B_{2m}$, which take into account the eccentricity effect, are obtained in practice by summation, using the following expression

$$B_{2m} = \sum_{i=1}^{N-1} (\cos^{2m}(\Theta_i) - \cos^{2m}(\Theta_{i+1})) \times \left(1 - \frac{l(\Theta_i) + l(\Theta_{i+1})}{2}\right), \qquad (30)$$

where $\Theta_i$ are the predefined angles at which the smoothed light curve is sampled, $N$ is defined by the constant step $\Delta\Theta = \Theta_{i+1} - \Theta_i$ adopted, and the value of $\Theta_1$ corresponds to the first contact of the eclipse. The $l(\Theta_i)$ values refer to the normalized smoothed light curve, and by definition $l = 1$ for $|\Theta| > |\Theta_1|$.

Since the individual data points are affected by observational errors, the integration of the empirical moments must itself be affected by errors. The uncertainty associated with the moments $B_{2m}$ can be evaluated using the following equation (see Al-Naimiy 1977; Smith & Theokas 1980)

$$\Delta B_{2m} = \frac{1}{\sqrt{n}} \left\{ \frac{1}{n} \sum_{j=1}^{n} \left[ l_j - l(\theta_j) \right]^2 \right\}^{1/2} \times \left( \cos^{2m} \theta_1 - \cos^{2m} \theta_n \right), \qquad (31)$$

where $n$ is the number of observed points over the considered eclipse, and $l_j - l(\theta_j)$ is the difference between the observed point of index $j$ and the smoothed light curve at $\theta_j$. The angles $\theta_1$ and $\theta_n$ refer, respectively, to the first and last observed data point over the relevant part of the light curve.

## 3.2. Solution for the elements

For each half-eclipse, there are five non-linear algebraic equations to be solved simultaneously for the elements $\rho_1$, $\rho_{2a}$ or $\rho_{2c}$,
$\rho_{3a}$ or $\rho_{3e}$, $i$, $u_1$ or $u_2$, $y_a$ or $y_e$, $v$, $L_1$ or $L_2$ the meanings of which are summarized below for convenience:

$L_1$ = luminosity of the OB-type star,

$L_2$ = luminosity of the W-R star,

$i$ = inclination angle of the orbit,

$\rho_1$ = radius of the OB-type star,

$\rho_{2a,e}$ = radius of the W-R envelope seen in absorption or emission,

$\rho_{3a,e}$ = radius of the W-R opaque core seen in absorption or emission,

$u_1$ = limb-darkening coefficient of the OB-type star,

$u_2$ = limb-darkening coefficient of the W-R envelope,

$v$ = transparency coefficient of the W-R envelope (cf. Fig. 3), and

$y_{a,e}$ = contribution of the W-R envelope in absorption or emission (cf. Fig. 3).

This large number of variables can fortunately always be reduced to a smaller number for both the primary and secondary minima, as explained hereafter. The preliminary determination of the orbit inclination $i$ already eliminates for instance one variable.

At the primary minimum, when the OB star is in front, according to our composite model $v \equiv 0$ by definition. Since we must also have that $0 \le u_2 \le 1$, solutions can be searched for a set of discrete values of the parameter $u_2$ in this interval, so that the remaining variables are $\rho_1$, $\rho_{2e}$, $\rho_{3e}$, $y_e$, and $L_2$.

At the secondary minimum, when the W-R star eclipses the OB component, $u_1$ is the limb-darkening coefficient of the OB star. The value of $u_1$ can be adopted following e.g. the tabulated values supplied by Klinglesmith & Sobieski (1970). According to our model, we now have $0 \le v \le 1$, so that again, after selection of a sample of $v$ values, we can proceed in solving the equations for the remaining parameters $\rho_1$, $\rho_{2a}$, $\rho_{3a}$, $y_a$, and $L_1$ only.

The solution of the system of as many as five non-linear equations

$$B_{2m}(\rho_1, \rho_{2a}, \rho_{3a}, i, u_1, y_a, v, L_1, L_2) = B_{2m}(\text{observed})$$

or

$$B_{2m}(\rho_1, \rho_{2e}, \rho_{3e}, i, u_2, y_e, v, L_1, L_2) = B_{2m}(\text{observed})$$

is obtained by minimizing the $\chi^2$ compiled from the residuals of these equations

$$\chi^2 = \sum_{m=1}^{N} \frac{|B_{2m}(\text{computed}) - B_{2m}(\text{observed})|^2}{\Delta B_{2m}^2} . \quad (32)$$

This minimization is achieved by means of Powell's technique (e.g., Press et al. 1992 and references therein).

# 4. Application to HD 5980

## 4.1. HD 5980: a peculiar system

HD 5980 $\equiv$ AB 5 (Azzopardi & Breysacher 1979) is associated with NGC 346, the largest H II region + OB star cluster in the Small Magellanic Cloud. This remarkable W-R binary, which underwent an LBV-type event in August 1994, is presently recognized as a key-object for improving our understanding of massive star evolution. HD 5980 is a rather complex system because it consists of at least three stars: two stars form the eclipsing binary with the 19.266 day period, whilst the third component, an O-star, which is detected by means of a set of absorption lines

and by means of its third light (see also below), could be a member of a highly eccentric 96.5 day period binary (Schweickhardt 2000; Foellmi et al. 2008). Whether or not the third star is physically bound to the eclipsing binary remains currently unclear. Before the LBV eruption, both components of the eclipsing binary already showed emission lines in their spectra and were thus classified as Wolf-Rayet stars (Niemela 1988). However, as shown by the analysis of the spectra taken during and after the LBV event, at least the star that underwent the eruption (hereafter called star A) was not a classical, helium-burning, Wolf-Rayet object, but rather a WNha star, i.e., a rather massive star with substantial amounts of hydrogen present in its outer layers (Foellmi et al. 2008). These WNha stars have wind properties that are intermediate between those of extreme Of stars and classical WN stars.

A summary of the light changes exhibited by HD 5980 was presented by Breysacher (1997). The technique of light curve analysis described above is applied to HD 5980 prior to the outburst. Given that star A, the component in front of its companion (hereafter called star B) during primary eclipse, was a WNha star and since there are no indications of wind effects in the primary eclipse, we assume that this component behaves as an OB-star in the light curve. 705 measurements with the Stroemgren $v$ filter, described in BP91, are taken into consideration. The shape of the resulting light curve does not allow us to use the ill-defined ascending branch of the primary eclipse (star A in front) for the analysis. Therefore, only three half-minima will be considered. Compared to the preliminary analysis carried out by BP91, the software tool presently used has been upgraded, allowing us for instance to assess the quality of the solution by means of a synthetic light curve.

## 4.2. Solutions of the light curve

From the smoothed light curve, the durations of the primary and secondary stellar core eclipses as well as the separation between the core eclipses are measured first to be $d_1 = 0.062 \pm 0.005$, $d_2 = 0.095 \pm 0.005$, and $\Delta\Phi = 0.362$. All durations are expressed as phase intervals (i.e. fractions of the orbital cycle). The corresponding values of the eccentricity $e$ and the longitude of periastron $\omega$ are: $e = 0.314 \pm 0.007$ and $\omega = 132.5^\circ \pm 1.5^\circ$. These values are in fairly good agreement with those derived for these parameters by Breysacher & François (2000) ($e = 0.30 \pm 0.02$, $\omega = 135^\circ \pm 10^\circ$) from a completely different approach based on the analysis of the width variation of the He II $\lambda$ 4686 line using the analytical colliding-wind model devised by Lührs (1997). From radial-velocity studies, Kaufer et al. (2002) and Niemela et al. (1997) also found that $e = 0.297 \pm 0.036$ and $e = 0.28$, respectively. The values of the moments of the light curve of HD 5980 are listed in Table 1.

The inclination angle of the orbit can easily be derived with reasonable accuracy. For each of the three half-minima, a quick analysis is carried out for a number of plausible discrete values of $i$ ($i = 82^\circ, 83^\circ, \dots, 89^\circ$), and the value of $i$ finally adopted is the one for which closest agreement is obtained between the three solutions provided for the radius of star A and the radius of the opaque core of star B. These quantities are fundamental elements of the system and the combination of the solutions of the descending and ascending branches of both minima must indeed provide, at the end of the detailed analysis, a unique value for $\rho_1$ and very similar – if not identical – values for $\rho_{3a}$ and $\rho_{3e}$. The above conditions are fullfilled for $85^\circ \le i \le 87^\circ$, therefore we adopt $i = 86^\circ \pm 1^\circ$.
<table><caption>**Table 1.** Values of the moments of the eclipses of HD 5980 used throughout this paper.</caption><thead><tr><th rowspan="2"></th><th>Primary eclipse</th><th colspan="2">Secondary eclipse</th></tr><tr><th>Ingress</th><th>Ingress</th><th>Egress</th></tr></thead><tbody><tr><td>B<sub>2</sub></td><td>.00935 ± .00026</td><td>.00832 ± .00041</td><td>.00481 ± .00024</td></tr><tr><td>B<sub>4</sub></td><td>.01801 ± .00048</td><td>.01542 ± .00065</td><td>.00922 ± .00043</td></tr><tr><td>B<sub>6</sub></td><td>.02604 ± .00067</td><td>.02164 ± .00080</td><td>.01326 ± .00059</td></tr><tr><td>B<sub>8</sub></td><td>.03351 ± .00083</td><td>.02721 ± .00088</td><td>.01699 ± .00071</td></tr><tr><td>B<sub>10</sub></td><td>.04046 ± .00097</td><td>.03226 ± .00094</td><td>.02044 ± .00082</td></tr></tbody></table>

These values are derived from the light curve prior to the 1994 LBV event.

The second step in the procedure consists of solving the equations for $i = 86^{\circ}$, each half-eclipse being treated in a completely independent manner. We recall that all radii are reduced to the semi-major axis of the relative orbit.

We first consider the descending branch of the primary eclipse. Solutions for $\rho_1, \rho_{2e}, \rho_{3e}, L_2$, and $y_e$ are searched for discrete values of the parameter $u_2$ (0.1, 0.2, 0.3, ... 1) with rather broad variation ranges allowed to the variable parameters. Each search sequence is based on a series of 1000 trials. A first set of results is obtained to provide convergence rates (i.e., an estimate of the likelihood of the derived solutions) and a $\sigma$ value for each parameter. A second iteration with reduced variation ranges (the original one $\pm \sigma$) for the parameters leads to a second set of solutions with improved convergence rates and lower $\sigma$ values. The process is repeated five times, until a stabilization of the parameter values becomes noticeable. For the last iteration, the grid of $u_2$ values is enlarged significantly and the respective convergence rates of the corresponding solutions are used to determine the best fit $u_2$ value and to estimate the error on this parameter. The obtained results are $\rho_1 = 0.150 \pm 0.004$, $\rho_{2e} = 0.257 \pm 0.018$, $\rho_{3e} = 0.110 \pm 0.005$, $L_2 = 0.300 \pm 0.016$, $u_2 = 0.58 \pm 0.07$, and $y_e = 0.19 \pm 0.03$.

For the secondary eclipse, solutions are searched for $\rho_1, \rho_{2a}, \rho_{3a}, L_1$, and $y_a$. The choice of the parameter $u_1$ for the limb-draking law (Klinglesmith & Sobieski 1970) of star A has little impact on the other parameters. We repeated the fitting procedure for different values of $u_1$ (0.1, 0.3, 0.5, 0.7) and recovered the same solutions within the errorbars. The largest sensitivity was found for $\rho_{2a}$ when $u_1 = 0.7$. In this rather unlikely case, $\rho_{2a}$ exceeds its usual value by $1.5\sigma$. In the following, we thus focus on the results obtained with $u_1 = 0.3$, which seems a reasonable value for star A (Klinglesmith & Sobieski 1970). The descending and ascending branches are analysed separately and solutions are searched for discrete values of the parameter $v$ (0.1, 0.2, 0.3, ... 1). The same iterative procedure as described above for the primary eclipse is applied. The resulting solutions for the secondary descending branch are $\rho_1 = 0.160 \pm 0.006$, $\rho_{2a} = 0.259 \pm 0.018$, $\rho_{3a} = 0.110 \pm 0.005$, $L_1 = 0.385 \pm 0.028$, $y_a = 0.19 \pm 0.03$, and $v = 0.60 \pm 0.20$; and for the secondary ascending branch $\rho_1 = 0.163 \pm 0.008$, $\rho_{2a} = 0.290 \pm 0.019$, $\rho_{3a} = 0.105 \pm 0.003$, $L_1 = 0.412 \pm 0.026$, $y_a = 0.20 \pm 0.03$, and $v = 0.45 \pm 0.15$.

The differences between the solutions provided by the three half-minima are relatively small compared to the errors and average values are derived for the parameters. A distinction between absorption and emission values does not appear to be necessary any longer. The error in the mean for each parameter $p_i$ is computed in a conservative approach using the expression

$$ \sigma = \left( \frac{\sum_{i=1}^{n} (p_i - \langle p \rangle)^2 + \sum_{i=1}^{n} \sigma_i^2}{nN} \right)^{1/2}, \quad (33) $$

<table><caption>**Table 2.** Final “best-fit” values of the model parameters of the stars in the HD 5980 binary system for $i = 86^{\circ}$ and prior to the 1994 LBV outburst.</caption><thead><tr><th></th><th>Primary eclipse<br>Ingress</th><th>Secondary eclipse<br>Ingress & Egress</th></tr></thead><tbody><tr><td>star A</td><td>$\rho_1 = 0.158 \pm 0.005$<br>$L_1 = 0.398 \pm 0.021$</td><td></td></tr><tr><td>star B</td><td>$\rho_3 = 0.108 \pm 0.003$<br>$\rho_2 = 0.269 \pm 0.014$<br>$L_2 = 0.300 \pm 0.016$<br>$y = 0.19 \pm 0.02$<br>$u_2 = 0.58 \pm 0.07$<br>$v = 0.52 \pm 0.14$</td><td></td></tr></tbody></table>

where $n$ is the number of values used in the mean and $N$ the number of independent sets (half-eclipses) of values. For an analysis completed in this way, the values adopted for the parameters of the stellar components in the HD 5980 binary system are given in Table 2.

### 4.3. Discussion

Figure 4 shows how the synthetic light curve derived from the above elements fits the observational data for both the primary and secondary eclipses. While the bottom and the wings are fairly well fitted, the transparency law that we adopted for the extended envelope of star B is probably still too crude to allow a perfect match to the observations of the descending and ascending sides of both minima. As a consequence, the size of this envelope is probably slightly overestimated by our model. An asymmetry in this envelope, inferred by BP91, is difficult to ascertain from the present study. The difference between the values of $\rho_{2a}$ provided by the two branches of the secondary eclipse is indeed marginally significant only given the errors. $L_1 + L_2 \neq 1$ is independent confirmation of an unresolved source of third light along the line of sight that accounts for an additional relative luminosity of $L_3 = 0.302$. In principle, one would expect a superior control of the uniqueness of the solution if the brightness ratios $L_2/L_1$ and $L_3/L_1$ in the optical could be fixed independently of our light curve analysis. However, in the specific case of HD 5980, it is impossible to infer these brightness ratios from spectroscopy. In this system, there have been changes in both the spectroscopic and in photometric properties of the binary components, and the brightness ratios are thus epoch dependent. Spectrophotometric brightness ratios would have to be estimated based on the dilution of the emission lines. However, as shown by recent spectroscopic studies (e.g. Foellmi et al. 2008), these emission lines contain epoch-dependent contributions from components A and B, as well as from the wind-wind interaction region. In our analysis of the light curve, we thus preferred not to constrain the brightness ratios a priori. Our results can however be checked a posteriori against the results obtained in the UV domain. Assuming that an O7 supergiant, which does not partake in the orbit, was responsible for this third light, Koenigsberger et al. (1994) estimated a value of 2.8 for the ratio $(L_1 + L_2)/L_3$, which was considered to be consistent with the 2.0 obtained from the BP91 elements. It is noticeable that the value of 2.3 derived from the present analysis is in even closer agreement with the estimate of Koenigsberger et al. (1994).

Adopting $M_v = -7.5$ for the global absolute visual magnitude of HD 5980 (cf. Koenigsberger et al. 1994), the magnitudes of the binary components are $M_v(\text{starA}) = -6.50$ and $M_v(\text{starB}) = -6.19$. Because of the difficulties in assigning the
<figure><img src="image_5.png" /></figure>

**Fig. 4.** *Top panel:* the observed light curve of HD 5980 in the Stroemgren *v* filter as a function of orbital phase. The *lower left and lower right panels* show the synthetic light curves for the primary and secondary eclipses respectively compared to the actual data. The synthetic light curves were computed using the mean parameters of the solutions found by our programme.

<figure><img src="image_6.png" /></figure>

**Fig. 5.** Distribution of the parameters derived from the inversion of the moments of the light curve in various parameter planes. The *left and middle panels* are derived from the primary eclipse, whilst the *rightmost panel* presents the solutions of the secondary eclipse. The final solutions derived from each eclipse as well as their $1-\sigma$ error bars are overplotted in each panel.

various lines of the spectrum of HD 5980 to a specific component, there have been few attempts to establish a full SB2 orbital solution (Niemela et al. 1997; Foellmi et al. 2008). Therefore, our knowledge of the component masses of the eclipsing binary remains uncomfortably poor, even though we are dealing with an eclipsing system. Foellmi et al. (2008) used a multi-component fit of the N IV λ 4058 and N V λ 4603 lines to derive absolute masses of 58–79 M⊙ for star A and 51–67 M⊙ for star B. These values differ significantly from the estimates of Niemela et al. (1997), who attributed the entire N IV λ 4058 line to star A and

the entire N V λ 4603 emission to star B when inferring minimum masses of $m \sin^3 i = 28$ and $50 M_\odot$, respectively. Adopting the solution of Foellmi et al. (2008), the sum of the masses of the components of the eclipsing binary equals 109–146 $M_\odot$. This corresponds to a semi-major axis of 0.67–0.74 AU (144–159 $R_\odot$) for the binary orbit, and to stellar radii of 22.7–25.1 $R_\odot$ for star A, 15.6–17.2 $R_\odot$ for the core of star B, and 38.7–42.8 $R_\odot$ for its envelope.

Alternatively, one could estimate absolute masses from the visual magnitudes evaluated above by means of mass-luminosity
relations derived e.g., from massive star evolutionary models. However, this approach is hampered by a lack of precise knowledge of the spectral types of the stars (and hence their bolometric corrections) and by the fact that the components of HD 5980 can probably not be considered as “normal” early-type stars.

We considered the correlations between the various free parameters by plotting them in pairs (see Fig. 5). For this purpose, we used the results of 10 000 trials for each eclipse. In most cases, we do not observe obvious correlations; the solutions are scattered over a limited part of the parameter plane (see e.g., the radius of the core of star B, $\rho_3$, versus that of star A, $\rho_1$, derived from the primary eclipse). However, the solutions show a clear trend, if we plot the luminosity of star A ($L_1$) as a function of the radius of the envelope of star B ($\rho_{2a}$, evaluated from the secondary eclipse), and fall even along a clearly defined locus if we plot the luminosity of star B ($L_2$) as a function of the radius of star A ($\rho_1$, derived from the primary eclipse). These trends can be understood at least qualitatively. In fact, for the primary eclipse (the opaque core of star A occulting star B), an increase in the radius of star A implies that a larger fraction of the light of star B is removed. To account for the observed depth of the light curve, the total luminosity of star B must decrease. On the other hand, during secondary eclipse (star B in front of star A), an increase in the radius of the semi-transparent envelope of star B will produce deeper and broader wings of the secondary eclipse. To account for the observed eclipse shape, the model must react in terms of an increase in the surface brightness (and hence the luminosity) of star A.

In principle, the photometric variability of HD 5980 as detected through medium-band filter observations could be affected by the Doppler shift of strong emission lines that fall in the wavelength range covered by the photometric filters. To quantify this effect, we simulated an observation of the star through the ESO Stroemgren filters. For this purpose, we used the spectrum of HD 5980 taken from the spectrophotometric catalogue of Morris et al. (1993) kindly provided to us by Dr. J.-M. Vreux. Our calculations indicate that less than 1% of the flux in the v band comes from emission lines. Therefore, we do not expect the Doppler shift to have any significant effect on the photometric variability discussed here. The situation would be quite different, if we were dealing with data taken in the b band. There, about 30% of the flux is produced by emission lines (especially the strong He II $\lambda$ 4686 line). In this case, the Doppler motion as well as variations in the line intensity caused by the wind-wind interaction (e.g., Breysacher et al. 1982; Breysacher & François 2000) could lead to significant variations in the observed b magnitude.

Finally, we note that the data points of the light curve suggest an increase in the brightness of the HD 5980 system at orbital phases after $\phi \sim 0.04$ (see Fig. 4). Whilst this increase could be caused by light reflections from the wind interaction region (which should have its concavity roughly turned towards the observer) near periastron, we caution that this apparent trend is actually inferred from observations from a single campaign (cycle 39 in Fig. 6) and could therefore be related to intrinsic variability of the WR star rather than represent a genuine phase-locked effect.

## 5. Conclusions

We have presented a method that allows us to treat the problem of atmospheric eclipses in the light curves of moderately wide, eccentric Wolf-Rayet + O binary systems in a semi-analytical way. We have then applied this method to the light curve of the

<figure><img src="image_7.png" /><figcaption>Fig. 6. Observed v magnitudes of HD 5980 in the phase interval 0.04–0.16. The various symbols indicate data from different orbital cycles. Periastron passage occurs at phase 0.061.</figcaption></figure>

peculiar system HD 5980 prior to its 1994 outburst. Despite the non-uniform sampling of the light curve, the fact that we are analysing data from different orbital cycles (hence affected differently by the intrinsic variability of the WR star) and the limitations of our assumptions on the properties of the WR envelope, our method yields consistent results when applied to the primary or secondary eclipse. We have been able to constrain some of the physical parameters of this system, although the lack of a consistent SB2 spectroscopic orbital solution prevents us from obtaining fully model independent parameters. As a next step, we will try to generalize our method to the analysis of the eclipses of a system harbouring two stars, both with extended atmospheres. This should allow us to analyse the light curve of HD 5980 observed after the LBV eruption.

## Appendix A: Analytical expressions for the moments of the eclipses

### A.1. Primary eclipse (OB star in front)

Here the relevant radii are: $\rho_1$, the radius of the OB-star; $\rho_{2e}$ and $\rho_{3e}$, the radii of the W-R core and envelope in emission. We note that we have assumed that $\rho_{2e} \le \rho_1$. Adopting a surface brightness law for the W-R star given by Eqs. (22) and (23) and $F(s) = 1$ (for $s \le \rho_1$), we find that its total luminosity is given by

$$L_{W-R} = J(0) \pi \rho_{2e}^2 \left[ (1 - y_e) \left( \frac{\rho_{3e}}{\rho_{2e}} \right)^2 + y_e \left( 1 - \frac{u_2}{2} \right) \right], \quad (A.1)$$

$$I_m = L_{W-R} \frac{(1 - y_e) \left( \frac{\rho_{3e}}{\rho_{2e}} \right)^{2m} + y_e \left( 1 - \frac{m u_2}{m + 1} \right)}{(1 - y_e) \left( \rho_{2e} \right)^2 + y_e \left( 1 - \frac{u_2}{2} \right)}, \quad (A.2)$$

and

$$P = L_{W-R} \frac{(1 - y_e) \left( \frac{\rho_{3e}}{\rho_{2e}} \right)^2 + y_e g \left[ 1 - \frac{u_2}{2} g \right]}{(1 - y_e) \left( \rho_{2e} \right)^2 + y_e \left( 1 - \frac{u_2}{2} \right)}, \quad (A.3)$$

where $g(\rho_1, \rho_{2e}) = \min(1, (\rho_{3e}/\rho_{2e})^2)$. The expression of P reduces to

$$P = L_{W-R} \quad \text{if } \rho_{2e} < \rho_1, \quad (A.4)$$
and finally,

$$R_m = 1. \tag{A.5}$$

## A.2. Secondary eclipse (W-R star in front)

Here $\rho_{2a}$ and $\rho_{3a}$ are the radii of the W-R core and envelope in absorption, respectively. Again, we have assumed that $\rho_{2a} \le \rho_1$. Adopting a transparency law for the W-R star given by Eqs. (20) and (21), we find that

$$L_O = J(0) \pi \rho_1^2 \left(1 - u_1 + \frac{u_1^2}{3}\right), \tag{A.6}$$

$$I_m = \frac{L_O}{1 - u_1 - \frac{u_1^2}{3}} \left(1 - \frac{2 u_1 m}{m + 1} + \frac{u_1^2 m}{m + 2}\right), \tag{A.7}$$

and

$$P = \frac{L_O}{1 - u_1 + \frac{u_1^2}{3}} \left\{ \begin{aligned} & \left(1 - y_a\right) \left( \frac{\rho_{3a}^2}{\rho_1^2} - u_1 \frac{\rho_{3a}^4}{\rho_1^4} + \frac{u_1^2}{3} \frac{\rho_{3a}^6}{\rho_1^6} \right) \\ & + y_a h \left[ \left(1 - u_1 h + \frac{u_1^2 h^2}{3}\right) - v g \left(\frac{1}{2} - \frac{2 u_1 h}{3}\right) \right. \\ & \left. \qquad + \frac{u_1^2 h^2}{4} \right] \end{aligned} \right\}, \tag{A.8}$$

where $g(\rho_1, \rho_{2a}) = \min(1, (\rho_1/\rho_{2a})^2)$ and $h(\rho_1, \rho_{2a}) = \min(1, (\rho_{2a}/\rho_1)^2)$.

Finally, we obtain

$$R_m = (1 - y_a) \left(\frac{\rho_{3a}}{\rho_{2a}}\right)^{2m} + y_a \left(1 - \frac{mv}{m+1}\right). \tag{A.9}$$

Acknowledgements. We are greatly indebted to Drs. S. A. H. Smith and A. C. Theokas for clarification about their method, and Profs. G. Koenigsberger and A. Maeder for fruitful discussions. G.R. is supported by the FRS/FNRS (Belgium)

and through the XMM/INTEGRAL PRODEX contract (Belgian Federal Science Policy Office) as well as by the Communauté Française de Belgique – Action de recherche concertée – Académie Wallonie – Europe.

## References

Al-Naimiy, H. M. K. 1977, Ap&SS, 46, 261

Antokhin, I. I., & Cherepashchuk, A. M. 2001, Astron. Rep., 45, 371

Azzopardi, M., & Breysacher, J. 1979, A&A, 75, 120

Barbá, R. H., Niemela, V. S., Baume, G., et al. 1995, ApJ, 446, L23

Bateson, F. M., & Jones, A. F. 1994, Pub. RAS New Zealand, 19, 50

Breysacher, J. 1997, in Luminous Blue Variables: Massive Stars in Transition, ed. A. Nota, & H. J. G. L. M. Lamers, ASP Conf. Ser. 120, 227

Breysacher, J., & François, P. 2000, A&A, 361, 231

Breysacher, J., & Perrier, C. 1980, A&A, 90, 207

Breysacher, J., & Perrier, C. 1991, in Wolf-Rayet Stars and Interrelations with other Massive Stars in Galaxies, Proc. IAU Symp. 143, ed. K. A. van der Hucht, & B. Hidayat (Dordrecht: Kluwer), 229

Breysacher, J., Moffat, A. F. J., & Niemela, V. S. 1982, ApJ, 257, 116

Cherepashchuk, A. M. 1975, SvA, 19, 47

Danjon, A. 1959, Astronomie Générale, ed. J. & R. Sennac (Paris)

Foellmi, C., Koenigsberger, G., Georgiev, L., et al. 2008, Rev. Mex. Astron. Astrof., 44, 3

Gosset, E., Remy, M., Manfroid, J., et al. 1991, IBVS, 3571

Heydari-Malayeri, M., Rauw, G., Esslinger, O., et al. 1997, A&A, 322, 554

Hoffmann, M., Stift, M., & Moffat, A. F. J. 1978, PASP, 90, 101

Kaufer, A., Schmid, H.-M., Schweickhardt, J., et al. 2002, in Interacting Winds from Massive Stars, ed. A. F. J. Moffat, & N. St-Louis, ASP Conf. Ser., 260, 489

Klingsmith, D. A., & Sobieski, S. 1970, AJ, 75, 175

Koenigsberger, G., Moffat, A. F. J., St-Louis, N., et al. 1994, ApJ, 436, 301

Kopal, Z. 1975, Ap&SS, 34, 431

Kopal, Z. 1979, Language of the Stars (Dordrecht and Boston: Reidel Publ. Co.)

Kopal, Z., & Al-Naimiy, H. M. K 1978, Ap&SS, 57, 479

Lamontagne, R., Moffat, A. F. J., Drissen, L., Robert, C., & Matthews, J. M. 1996, AJ, 112, 2227

Lührs, S. 1997, PASP, 109, 504

Morris, P. W., Brownsberger, K. R., Conti, P. S., Massey, P., & Vacca, W. D. 1993, ApJ, 412, 324

Niemela, V. S. 1988, in Progress and Opportunities in Southern Hemisphere Optical Astronomy, ed. V. M. Blanco, & M. M. Phillips, ASP Conf. Ser., 1, 381

Niemela, V. S., Barbá, R. H., Morrell, N. I., et al. 1997, in Luminous Blue Variables: Massive Stars in Transition, ed. A. Nota, & H. J. G. L. M. Lamers, ASP Conf. Ser., 120, 222

Press, W. H., Teukolsky, S. A., Vetterling, W. T., et al. 1992, Numerical Recipes in Fortran 77 (Cambridge University Press)

Schweickhardt, J. 2000, Ph.D. Thesis, Landessternwarte Heidelberg/Königsstuhl

Smith, S. A. H. 1976, Ap&SS, 40, 315

Smith, S. A. H., & Theokas, A. C. 1980, Ap&SS, 70, 103

Sterken, C., & Breysacher, J. 1997, A&A, 328, 269