© ESO, 2013

1. Introduction

Studing massive stars is important because they play a key role in galaxy evolution: their strong winds interact with the ambient interstellar medium and can trigger the formation of new stars, they produce heavy chemical elements, and they are an important source of UV radiation. Recent research indicates that a large part of massive stars (~50%) form a binary or a multiple system (Mahy et al. 2009; Sana & Evans 2011). The spectra of binary systems are an invaluable source of information for determining the physical properties (such as masses, temperatures, and radii) of stars. It has been shown, however, that the observational analysis of massive binaries is complicated by effects that are linked to interactions in these systems that affect the spectra and the spectral classification (Sana et al. 2005; Linder et al. 2007; Linder 2008, and references therein). Unfortunately, the majority of spectral modelling codes are designed for isolated spherical stars and cannot reproduce the particular effects of the binarity. Therefore, it is important to improve our models of spectral computation to accurately represent the spectra of these stars. In this context, we have developed a method that takes into account some of the effects that are produced by gravitational interactions in the system.

The first models of binary systems in which the stars are not spherical have been proposed by Russell & Merrill (1952) to reproduce the light curves of binaries. Then, Kopal (1959) introduced the Roche potential approach that was first used by Lucy (1968) and Wilson & Devinney (1971) followed by numerous other works. In a first step, we designed a similar code for circular massive binary systems (Palate & Rauw 2012, hereafter Paper I) in which the distorted shapes of the stellar surfaces are caused by gravitational interaction. The emerging spectra are then computed accounting for the shape of the stars and their local properties. In this paper we extend the method to incorporate radiation pressure effects on the stellar surface and the case of eccentric binaries. This new version of our CoMBiSpeC (code of massive binary spectral computation) model now allows spectral computation of (almost) any massive binary system.

2. Improvements of the method

Our method for modelling the spectra of circular binary systems, presented in Paper I, consists of computing the actual stellar surface following the Roche potential approach and accounting for the local gravity and temperature on each surface element. Non-local thermodynamic equilibrium (NLTE) spectral model grids (TLUSTY, Lanz & Hubeny 2003, 2007) are used to compute the integrated spectrum of the star at each orbital phase (see also, e.g., Linnell & Hubeny 1994 and Linnell et al. 2012). Because our spectra are synthetic, they are free of observational noise. In this section we describe the two improvements that have been implemented to the method. The first is the inclusion of radiation pressure effects; i.e., the irradiation of each star by the companion and the effect of the radiation on each star’s own surface. The second is the use of an explicit calculation of the stellar surface, from first principles, which here is applied to the case of eccentric binaries. In eccentric systems, the orbital separation changes as a function of phase, and accordingly, the deformations of the surface and the orientation of the tidal bulge are also variable.

2.1. Radiation pressure effects

There exists quite an extensive body of literature (see e.g. Dermine et al. 2009; Howarth 1997a; Schuerman 1972; Drechsel et al. 1995 and Phillips & Podsiadlowski 2002) on the impact of radiation pressure on the Roche potential in early-type binaries. According to Dermine et al. (2009), the effects of radiation pressure in a binary system can be separated into three parts: radiation pressure on each star’s own surface, the effect on the companion, and radiation pressure on the matter outside the binary orbit. Here we focus on the two first effects, which are called internal and external radiation pressure (Drechsel et al. 1995). The impact of the radiation pressure on the matter outside the binary is, of course, very important for the surrounding medium but is irrelevant for calculating the shape of the stars. Howarth (1997a) pointed out that the effect of the internal radiation pressure can be treated as a simple scaling of the Roche potential of the stars. Another approach, proposed first by Schuerman (1972), consists of scaling the mass of the stars rather than the entire potential. The former approach is more appropriate because the emitted radiation bolometric flux and hence the radiation pressure scale with the local gradient of the potential (which also includes the attraction by the other star and the centrifugal forces). Therefore, according to Howarth (1997a), the stellar topology is not affected by internal radiation pressure and the latter simply scales the local surfaces gravity. The Newtonian gravity is defined by ${g}\mathrm{=}\mathrm{-}{\nabla }\mathit{\phi ,}$(1)where φ is a scalar gravity potential per unit mass.

The von Zeipel (1924) theorem can be written¹${\mathit{T}}^{\mathrm{4}}\mathrm{\propto }\parallel {g}\parallel \mathit{.}$(2)The magnitude of the radiative acceleration is given by (Howarth 1997a) ${\mathit{a}}_{\mathrm{rad}}\mathrm{=}\mathrm{-}\frac{\mathit{\kappa \pi }}{\mathit{c}}\mathrm{ℱ}\mathit{,}$(3)where κ is the flux mean opacity per unit mass and ℱ is the astrophysical bolometric flux (ℱ ∝ T⁴).

If we define Γ as the ratio of the radiative to the gravitational acceleration, a_rad = Γg. According to Howarth (1997a), this ratio is constant for a given star and independent of the position on the stellar surface since both accelerations vary with T⁴. Therefore, we have that the effective potential is the potential without radiation pressure minus the radiative acceleration ${g}\mathrm{eff}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Gamma }\mathrm{)}{g}\mathit{.}$(4)From Eq. (1), φ_eff = (1 − Γ)φ and the ratio Γ can be written as (following the Castor et al. 1975 theory for evaluating κ) $\begin{array}{c}\\ \mathrm{\Gamma }& \mathrm{=}& \frac{\mathit{\kappa \pi }}{\mathit{c}}\mathrm{ℱ}\mathit{/}\mathit{g}\\ & \mathrm{=}& \frac{{\mathit{\sigma }}_{\mathrm{Th}}}{{\mathit{m}}_{\mathrm{H}}\mathit{c}}\mathit{\sigma }{\mathit{T}}^{\mathrm{4}}\mathit{/}\mathit{g,}\end{array}$(5)with $\frac{{\mathit{\sigma }}_{\mathrm{Th}}}{{\mathit{m}}_{\mathrm{H}}}\mathrm{\approx }\mathrm{0.036}$ m² kg^-1, σ_Th the Thomson-scattering cross section and σ the Stefan-Boltzmann constant.

Because Γ is independent of the position on the stellar surface, we can compute it for T = T_pole and g = g_pole$\mathrm{\Gamma }\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{Th}}}{{\mathit{m}}_{\mathrm{H}}\mathit{c}}\mathit{\sigma }{\mathit{T}}_{\mathrm{pole}}^{\mathrm{4}}\frac{\mathrm{1}}{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}{g}\mathrm{pole}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}\mathrm{·}$(6)The advantage of this method is that the von Zeipel theorem remains valid. Another advantage is the simplicity of the method which does not require many computational resources.

The external radiation pressure is more difficult to treat. We have based our method on the approach of Drechsel et al. (1995) and Phillips & Podsiadlowski (2002). This method consists of scaling the mass of the companion in the Roche potential. The scale parameter $\mathit{\delta }\mathrm{=}\frac{{\mathit{F}}_{\mathrm{rad}}}{{\mathit{F}}_{\mathrm{grav}}}$ is computed iteratively for each surface point. The external radiation pressure can be seen as a force that decreases the attraction of the companion. Therefore, the companion seems less massive, and consequently the mass has to be scaled. The new “Roche” potential can be written $\begin{array}{c}\\ \mathrm{\Omega }& \mathrm{=}& \frac{\mathrm{1}}{\mathit{r}}\mathrm{+}\frac{\mathit{q}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\delta }\mathrm{(}\mathit{r,\theta ,\varphi }\mathrm{\left)}\mathrm{\right)}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{r}\cos \mathit{\varphi }\sin \mathit{\theta }\mathrm{+}\mathrm{1}}}\mathrm{+}\frac{\mathit{q}\mathrm{+}\mathrm{1}}{\mathrm{2}}\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }\\ & & \mathrm{-}\mathit{qr}\cos \mathit{\varphi }\sin \mathit{\theta }\end{array}$(7)where $\mathit{r}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$, $\mathit{q}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}$, x = rcosϕsinθ, y = rsinϕsinθ, and z = rcosθ. Here, θ and ϕ are the colatitude and longitude angle in the spherical coordinates centred on the star. θ = 0 and ϕ = 0 correspond respectively to the north pole of the star and the direction towards its companion.

The radiation pressure can be written ${\mathit{P}}_{\mathrm{rad}}\mathrm{=}\frac{\mathrm{1}}{\mathit{c}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{\mathit{\omega }}{\mathit{I}}_{\mathit{\nu }}\hspace{0.17em}{\cos }^{\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{1}}\hspace{0.17em}\mathrm{d}\mathit{\omega }\mathrm{d}\mathit{\nu ,}$(8)where Ψ₁ is the angle between a local surface normal on the irradiated star and the direction to a surface element on the other star. I_ν is the specific intensity of the external radiation field, and ν is the frequency.

The radiation pressure gradient per unit mass is equal to (in plane-parallel approximation) $\begin{array}{c}\\ {\mathit{F}}_{\mathrm{rad}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathrm{d}{\mathit{P}}_{\mathrm{rad}}}{\mathrm{d}\mathit{r}}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathit{c}}{\mathrm{\int }}_{\mathit{\omega }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{\kappa }}_{\mathit{\nu }}\frac{\mathrm{d}{\mathit{I}}_{\mathit{\nu }}}{\mathrm{d}{\mathit{\tau }}_{\mathit{\nu }}}{\cos }^{\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{1}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{d}\mathit{\omega }\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathit{c}}{\mathrm{\int }}_{\mathit{\omega }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{I}}_{\mathit{\nu }}\hspace{0.17em}\cos {\mathrm{\Psi }}_{\mathrm{1}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{d}\mathit{\omega }\mathit{.}\end{array}$(9)Following Castor et al. (1975), we assume that the dominant source of opacity comes from electron scattering, thus, we have ${\mathit{\kappa }}_{\mathit{\nu }}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{Th}}}{{\mathit{m}}_{\mathrm{H}}}$.

The value of δ is therefore computed with an iterative method. We start from the classical unmodified Roche potential and compute the radiation received by each surface element from its companion. The radiation is computed by summing the contribution of all points that are visible from the primary (resp. secondary) on the secondary (resp. primary). The method applied to compute the received bolometric flux by a point is derived from Eq. (9) and is similar to the reflection effect treatment of Wilson (1990). Indeed, we have I_ν = I_ν,0 × (1 − u + ucosΨ₂), where Ψ₂ is the angle between a local surface normal on the emitting star and the direction towards an irradiated surface element on the other star, and u is the linear limb-darkening coefficient based on the tabulation of Claret & Bloemen (2011). Moreover, we can write $\mathrm{d}\mathit{\omega }\mathrm{=}\frac{\mathrm{d}{\mathit{S}}_{\mathrm{2}}\cos {\mathrm{\Psi }}_{\mathrm{2}}}{{\mathit{s}}^{\mathrm{2}}}$, where s is the distance between the irradiated and emitting surface. Finally, if the radiation comes only from the companion, we have ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{I}}_{\mathit{\nu ,}\mathrm{0}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{=}\frac{\mathit{\sigma }{\mathit{T}}^{\mathrm{4}}}{\mathit{\pi }}$. The local δ factor is thus defined by $\mathit{\delta }\mathrm{=}\frac{{\mathit{F}}_{\mathrm{rad}}\mathrm{\ast }{\mathit{d}}^{\mathrm{2}}}{\mathit{GM}}$, where d is the distance between the point and the centre of the emitting star. This factor can be written (for the radiation emitted by the secondary and received by the primary) as ${\mathit{\delta }}_{\mathrm{2}}\mathrm{=}\frac{\mathit{\sigma }{\mathit{\sigma }}_{\mathrm{th}}{\mathit{d}}^{\mathrm{2}}}{\mathit{\pi c}{\mathit{m}}_{\mathrm{H}}\mathit{G}{\mathit{M}}_{\mathrm{2}}}\begin{array}{}\sum \\ \cos {\mathrm{\Psi }}_{\mathrm{1}}\mathit{>}\mathrm{0}\cos {\mathrm{\Psi }}_{\mathrm{2}}\mathit{>}\mathrm{0}\end{array}{\mathit{T}}_{\mathrm{2}}^{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{u}\mathrm{+}\mathit{u}\cos {\mathrm{\Psi }}_{\mathrm{2}}\mathrm{)}\cos {\mathrm{\Psi }}_{\mathrm{1}}\cos {\mathrm{\Psi }}_{\mathrm{2}}\frac{\mathrm{d}{\mathit{S}}_{\mathrm{2}}}{{\mathit{s}}^{\mathrm{2}}}\mathrm{·}$(10)Then, with the δ factor, we recompute the stellar surface with the modified Roche potential. Again, we evaluate the δ factor and surface until we reach convergence.

2.2. Eccentricity

The extension to eccentric systems uses the TIDES² code (Moreno et al. 1999, 2005, 2011). This code computes the time-dependant shape of the stellar surface for eccentric and/or asynchronous systems. It also provides velocity corrections for rigid body rotation. The method consists of modelling a deformable and perturbed surface layer that lies upon a rigid body that is in uniform rotation by solving the equations of motion of a grid of surface elements. These equations take into account centrifugal and Coriolis forces, gas pressure, and viscous effects in addition to those of the gravitational potential of both stars. It is important to note that in asynchronously rotating or eccentric systems, the linear approximation of the Roche potential is no longer valid and the presence of viscosity in the medium may lead to non-linear effects³. The surface shape and local velocities are obtained as a function of orbital phase from the solution of the equations of motion. Full details of the model are given in Moreno & Koenigsberger (1999), Toledano et al. (2007), and Moreno et al. (2005, 2011).

With the derived values of radius and velocity perturbations, we compute the surface gravity and temperature. Assuming the perturbations to be small, the gravity is computed by using a classical Roche gradient.

The temperature distribution is computed following the von Zeipel theorem ${\mathit{T}}_{\mathrm{local}}\mathrm{=}{\mathit{T}}_{\mathrm{pole}}{\left(\frac{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\underline{\mathrm{\nabla }}{\mathrm{\Omega }}_{\mathrm{local}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\underline{\mathrm{\nabla }}{\mathrm{\Omega }}_{\mathrm{pole}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}\right)}^{\mathrm{0.25}\mathit{p}}\mathit{,}$(11)where the gravity-darkening parameter p = 1 in the case of massive stars.

We also include reflection effects that heat the surface of the star through irradiation of each component by the other. The effect was treated by following the approach of Wilson (1990). We neglected any cross-talk between the surface elements (e.g. due to horizontal advection or radiative exchanges between neighbouring surface elements of different temperature). This is probably an important approximation because this effect might smooth out some of the temperature variations, especially near periastron (for eccentric systems) where the stars move faster and where the gravitational interactions are stronger. These interactions and velocity imply that energy dissipation due to viscosity effects is more significant near periastron and thus the cross-talk is expected to be higher.

The spectral computation was performed with the CoMBiSpeC code, which was described in Paper I. For this paper, some minor changes have been applied to take into account particularities linked to the eccentricity, such as the variation of the separation between the two stars.

3. Circular-orbit models

We present here models of four binary systems with circular orbits. Three of them were previously studied in Paper I, i.e., models 1, 2, and 3. These models are based on the HD 159176, HD 165052, and HD 100213 systems analysed by Linder et al. (2007). The fourth system, model 4 in this paper, is based on Sk-67°105 studied by Bonanos (2009). All four of these systems are in synchronous rotation and, thus, are in equilibrium configuration so that the Roche potential approach is valid. For all models, we compared the shape, temperature, and gravity variations caused by radiation pressure. Table 1 gives a summary of the parameters (from Linder et al. 2007 and Bonanos 2009) used for the computation. We have chosen the polar temperature (rather than a surface mean temperature) as the reference effective temperature in Table 1. This is because T_pole is used to compute the local temperature via the von Zeipel theorem and remains essentially unaffected by radiation pressure effects. Finally, for model 4, we investigated the Struve-Sahade effect for some lines.

3.1. Models 1 and 2

We have analysed the spectra⁴ of models 1 and 2 at 20 orbital phases and found that the impact of radiation pressure is rather low. The comparison was made by subtracting the normalized spectra with radiation pressure to the spectra without radiation pressure. On each of these “residual” spectra, we searched for the maximum, mean, and median values and the standard deviation. Then we searched for the maximum of the maxima (over the orbital phase) and computed a mean of the medians, the standard deviations, and means. We also looked for the spectral bins that present deviations higher than 1% (over a total of 360 000 bins). Table 2 gives a summary of these values for the combined spectra of the studied models. For example, the highest differences are 0.0113 and 0.0191 for models 1 and 2. These maxima occur in both cases at phase 0.5, i.e. at the conjunction with the primary in front of the secondary. These differences are significant but, the mean⁵ difference is only  ~7−9 × 10^-4, which is very small. The general appearance of the spectra is preserved and the conclusions reached in Paper I remain valid. We can see in Fig. 1 that the spectra of model 1 with and without radiation pressure are nearly indistinguishable (top panels).

The largest and mean radii of the stars are smaller if radiation pressure effects are included. The radiation pressure also decreases the mean log (g) by  ~0.05 despite the smaller radii. Finally, the mean temperatures increase by 100 K. This increase of the temperature and decrease of gravity cause a weak increase of a few per cent of the total flux emitted by the stars. Figure 1, lower panel, displays the log  of the ratio of the fluxed spectra with and without radiation pressure and underlines the small differences in the continuum. Table 3 gives a summary of the radii and temperatures of circular models.

Fig. 1 Top: comparison between the normalized spectra computed without radiation pressure (in black) and with radiation pressure (in grey) for model 1. To better distinguish spectra, the black line is wider than the grey line. Orbital phases are, from bottom to top: phases 0.1, 0.25 and 0.4. The spectra are shifted vertically by 0.25 continuum units for clarity. Middle: corresponding residuals (spectrum with radiation pressure minus spectrum without radiation pressure). The residual plots are shifted by 0.015 for clarity. Bottom: logarithm of the ratio of the fluxed spectra with radiation pressure divided by the one without. These plot are shifted by 0.006 for clarity.

3.2. Model 3

This model leads to an (over)contact binary if radiation pressure is not included, but when it is included, the system becomes detached⁶. Surprisingly, however, the impact of radiation pressure on the spectrum is weak. This is because the parts of the star that contribute most to the spectrum come from the rear and side parts of the star which are less modified by the radiation pressure because δ is small (see Fig. 2). Moreover, the shape of the star is not strongly affected by the radiation of the companion. This is a priori unexpected but the reason for it is that the most deformed parts of the stars mainly see the coldest part of the companion, which leads to a rather low value of the maximum δ. We stress here that in systems where the two stars are in (over)contact (or nearly so) and are of similar size and temperature, the reflection effect is not sufficient to counterbalance the gravity darkening. Thus, the coldest parts of the stars are those near the L1 point.

In Paper I, we pointed out that we failed to reproduce the surface temperature distribution deduced by Linder et al. (2007) for HD 100213. These authors showed in their analysis through radial velocity measurements that the He II and He I lines were not formed in the same region of the stellar surfaces. The small changes of the stellar surface induced by the radiation pressure increase the importance of the reflection process and are sufficient to now explain the surface temperature distribution observed by Linder et al. The He II lines have a smaller radial velocity amplitude and thus are apparently stronger in the hemisphere facing the companion, whilst the He I lines have larger velocity amplitudes, indicating that they are stronger on the opposite hemisphere. This suggests that the hemisphere facing the companion is hotter than the opposite one. When radiation pressure is added, the stars are less deformed. This implies a weak temperature increase of the facing hemispheres caused by the decreasing gravity darkening and a higher effectiveness of the reflection process. This weak increase is therefore sufficient to explain the observations.

Fig. 2 Value of the radiation pressure parameter, δ, over the stellar surface for model 3. The highest values are located near the L1 point.

Espinosa Lara & Rieutord (2011) and Maeder (1999) have suggested that von Zeipel’s (1924) classical theorem requires modifications. These theoretical studies present more complex models that could be, at first approximation, equivalent to the classical von Zeipel theorem with a gravity-darkening parameter (GDP) p smaller than 1.00. Therefore, we chose this very deformed system in which the von Zeipel theorem implies strong temperature variations to test the influence of p on the spectra. We computed the spectrum for values of the GDP of 0.25, 0.50, 0.75, 1.25 and 1.50, all other parameters having the same value as model 3 in Table 1. We found that for a moderate change of p (of  +/− 0.25) the spectra are not notably affected. We also found that the variations remain within the noise level of real data for variations of up to  +/− 0.50. Therefore, since we made some other important assumptions (no cross-talk, assumptions underlying the Roche potential), it is not useful to adopt another more complex formalism for gravity darkening because it does not notably affect our simulations.

3.3. Model 4

This system was chosen as an example of one in which the radiation pressure considerably modifies the shape of the stars (see Fig. 3), underlining the potential significance of this effect. The system parameters (see Table 1) are based on the observational study of Sk-67°105 that was carried out by Bonanos (2009).

Fig. 3 Surface modification by the radiation pressure for model 4. Top: star surface without radiation pressure effect. Bottom: star surface with radiation pressure effect.

3.3.1. Spectral classification

Observational studies have shown that the primary and secondary are an O7-8V and an O8-8.5 III-V star. We first determined the spectral type of our synthetic spectra with the Conti (Conti & Alschuler 1971) – Mathys (1988, 1989) criterion (hereafter Conti-Mathys criterion). This criterion is based on the ratio of the equivalent widths (EWs) of He I λ 4471 to He II λ 4542 for the spectral-type determination and on the ratio of the EWs of Si IV λ 4089 to He I λ 4143 for the luminosity-class determination. We found that our model spectra that include radiation pressure correspond to an O6.5-7I type for the primary spectrum and an O7.5-O8I type for the secondary. The determination of the spectral type gives the same result whether the measurement is made on individual spectra or on combined spectra⁷. In our model, the stars seem hotter, and the problem encountered in Paper I concerning the luminosity class is again present because we overestimate the luminosity. In a second approach, we used the Walborn & Fitzpatrick (1990) atlas to compare individual spectra to the reference atlas and thus re-classify the stars. With the atlas we found for the primary and secondary an O(6.5)-7 V and an O8-8.5 (III)-V classification. The classification using the atlas is of course more qualitative, therefore we give a range and the best correspondence is in brackets. Within the uncertainties, the agreement with the observations is good.

3.3.2. Struve-Sahade effect

The Struve-Sahade effect (hereafter S-S effect) was originally defined as the apparent strengthening of the secondary spectrum when the star is approaching the observer and its weakening as it moves away (Howarth et al. 1997b; Linder et al. 2007). A more general definition is the apparent variation of the line strengths of either of the binary components as a function of the orbital phase (Howarth et al. 1997b). Bonanos (2009) reported a quite considerable S-S effect for several spectral lines of Sk-67°105. We studied eight of these lines in detail: He I λλ 4026, 4143, 4471, 4713, 5016, Si IV λ 4089, He II λλ 4542, and 5411.

First, we analysed the EWs on individual spectra of the primary and the secondary⁸. The measurements of the line EWs were performed with the MIDAS software developed by ESO. The EWs of lines in the spectra of the individual stars were determined directly by simple integration with the integrate/line MIDAS routine, whilst for the binary system, we used the deblend/line command as we would do for actual observations of a real binary spectrum. The latter routine fits two Gaussian line profiles to the blend of the primary and secondary lines.

The EWs of the different lines measured on the individual spectra display phase-locked variations. For the primary, the EWs of the He I lines reach a maximum at phase 0.5 (when the primary eclipses the secondary) because the coldest part of the star is the rear part. Two minima are visible at phases 0.05 − 0.1 and 0.9–0.95 because the observer starts to see the front part of the primary star. This part of the star is the hottest because of the very effective reflection in this system. The minimum is not observed at phase 0 because of the roughly annular eclipse that hides a substantial part of the stellar surface. The variation amplitude is line-dependent with  ~5% for He I λ 4026, 60% for He I λ 4143, and 20% for the other He I lines. As expected for heating effects, the He II lines exhibit a reverse variation compared to the He I lines (see Fig. 4). The Si IV line displays a similar variation as the He I lines.

For the secondary, the EWs display the same variation as the primary, but shifted by 0.5 in phase. We can notice the zero value at phase 0.5 when the secondary is totally eclipsed. The variation amplitude of the He I λλ 4026, 4713, 5016, and Si IV lines is small and less than  ~10%. The amplitude is greater than  ~20% for the other lines (we did not account for the zero value).

The measurements of the combined spectra using a two-Gaussian fit to the blended lines agree relatively well with the variations measured on individual spectra for several lines: He I λλ 4026, 4143, 4713, 5016, He II λλ 4542, and 5411. The variation amplitude is sometimes overestimated, though. The EWs of He I λ 4471 are asymmetric before and after phase 0.5. The variation does not seem regular but is similar to the variations measured on individual spectra, i.e. a decrease of the EW of the secondary and an increase for the primary during the first half of the orbital cycle. Finally, the Si IV line displays irregular variations.

Bonanos (2009) reported a strong S-S effect in Sk-67°105 for several lines that we studied in the model. The variations observed are consistent with the definition of Linder et al. (2007), i.e. in the first half of the orbital cycle the secondary (resp. primary) line is deeper than the primary (resp. secondary) line and the situation is reversed in the second half of the orbital cycle. However, even though some lines in our synthetic spectra display strong EW variations, we did not observe this type of phase dependence.

Fig. 4 Example of EW variations of the He II λ 4542 line. The EWs have been normalized to the highest value. Black: primary. Grey: secondary. Left: variations measured on individual spectra (EWp,max = 0.53 Å, EWs,max = 0.47 Å). Right: variations measured on combined spectra by fitting two Gaussians to the blended lines (EWp,max = 0.52 Å, EWs,max = 0.16 Å).

3.3.3. Radial velocity and synthetic broad-band light curve

The last two characteristics that we studied are the light curve in the wavelength range 3800–7100 Å and the radial velocity curve. The two eclipses are clearly visible in the light curve. The depth of the eclipses agrees well with the observations. We stress that light curves are “by-products” of the spectral computation and that CoMBiSpeC is not primarily designed for light curve computation. The semi-amplitude of the radial velocity also agrees well with the observations. The radial velocity was computed by taking the mean value of the velocity of the visible points of the stellar surface at a given phase weighted by the surface projected along the line of sight. The values are given for the primary and the secondary (observational value in brackets): 138.3 (137) km s^-1 and 332.9 (326) km s^-1. Finally, the synthetic light curve and radial velocity curve agree well with those derived by Bonanos (2009).

4. Eccentric models

We present here five models of eccentric binary systems. These models are inspired by the following massive binaries: HD 93205, HD 93403, HD 101131, HD 152218, and HD 152248, and are named models E1, E2, E3, E4, and E5. The stars were studied by: Antokhina et al. (2000), Rauw et al. (2000), Gies et al. (2002), Sana et al. (2008), Mayer et al. (2008), and Sana et al. (2001). The parameters used for the computation are given in Table 4. We studied the spectral classification, the radial velocity curves, the broad-band light curve, and finally the S-S effect. We computed the surface of these stars with the TIDES code (Table 4 bottom gives the specific parameters used in the TIDES code) and used a modified version of the CoMBiSpeC code to compute the gravity and temperature. These quantities are not uniformly distributed across the stellar surface and, in addition, change over the orbital cycle.

Fig. 5 Variations of the mean, highest and lowest visible radius (left: the primary radii, right: the secondary radii). Top: variations for model E1. Bottom: variations for model E3.

Fig. 6 Variations of the mean, highest and lowest visible radius, temperature, and log (g) for model E4. Left: values for the primary star. Right: values for the secondary star.

A sample of the phase-dependent behaviour of the mean radius, gravity, and temperature as well as their maximum and minimum values are shown in Figs. 5–7. The maximum radius is largest near periastron (just before and/or just after). A model based on an instantaneous Roche lobe, such as the Nightfall⁹ program, presents the strongest variation strictly at the periastron passage. However, in the TIDES code, viscous effects that lead to a delay are taken into account. The most compact systems, model E1, the secondary of model E2, and model E3, remain nearly spherical during the entire orbital cycle (see Fig. 5) because the mean radii are nearly constant and the highest relative difference between the maximum and minimum radii is less than 0.5%. This is not surprising because models E2, and E3 have the longest periods, hence the widest separation, and model E1 contains relatively compact main-sequence stars.

Fig. 7 Variations of the mean, highest and lowest visible radius and temperature for model E5. Left: values for the primary star. Right: values for the secondary star.

4.1. Spectral classification

First, we applied the quantitative Conti-Mathys criterion (Conti & Alschuler 1971, Mathys 1988, 1989). The spectral types that we found for our simulated spectra either agree well with the observations or are hotter than the observations. The hotter stars present a weak He I λ 4143. This makes measuring the spectra of this system difficult and uncertain. The luminosity classes are, as in Paper I, overestimated for models E1, E3, and E5. Model E2 and the primary of model E1 are too hot for the Conti-Mathys luminosity criterion. The luminosity classes agree well for model E4. The problem of simulating singlet transition of He I reported by Najarro et al. (2006) and discussed in Paper I is probably responsible for the general luminosity classification problem. To refine our classification, we again used the Walborn & Fitzpatrick (1990) atlas. This atlas was used with the individual spectra, which for real observational data can only be assessed through spectral disentangling. The atlas is more qualitative than the Conti-Mathys criterion but relies on a large number of lines. We can now achieve a better agreement with the classification of the real binary systems (see Table 5).

4.2. Struve-Sahade effect

In our sample of binaries, two are known to present this effect: HD 93403 and HD 152248. They inspired our models E2 and E5. In our analysis, we detected variations of the line strength in models E1, E2, E4, and E5. However, we observe an S-S effect in the sense of Linder et al. (2007) only for some lines in model E2.

For all models, we investigated eight “well-chosen” lines in the spectra at 20 phases (the phase zero corresponds to the periastron passage). The lines were chosen following three criteria: first, we selected the lines for which an S-S effect is reported in the literature. The second criterion is based on a visual detection of variations in the line profile in the synthetic spectra, and finally, the third criterion is based on the fact that some lines present the S-S effect more often than others. As in Paper I, we draw attention to possible blends with nearby lines that can modify the strength and lead to misinterpretations.

We followed the same procedure as in Paper I and measured the EWs in the simulated combined spectra of the binary at different phases as well as in the simulated spectra of the individual components of the binary at the same phases. This approach allows us to compare the line strengths that are deblended from the combined spectrum with the actual individual spectra, which is not possible with the real observational data. We found variations in the EWs in many of the individual spectra as well as in the combined spectra. However, the agreement between the results of the deblend/line routine and the actual individual spectra is mostly poor. We present in Figs. 9 and 10 examples of typical variations observed in our analysis. The (inverse) U-pattern (Fig. 9 right, grey or Fig. 10 left, black) is often observed and is directly related to the orbital motion. Constant EWs or irregular variations are also often observed. A constant EW is mostly observed in long-period systems. We also encountered other types of variation patterns such as those presented in Fig. 10 (right, grey) during our analysis. The amplitude of the EW variations depends on the line and the system. The variations observed on individual spectra can be explained by the physical variations of the temperature and gravity at the stellar surface during the orbital cycle (we give an example of the interpretation of the variation in Sect. 4.2.1). The explanation of the S-S effect that we suggested in Paper I seems to be reinforced by these new analyses. The S-S effect seems to be due to the combination of the spectra in which the lines have a non-Gaussian/non-Lorentzian and even asymmetric profile. Under these conditions, the deblend/line routine introduces systematic errors that mimic the S-S effect. This is illustrated in Fig. 8, which shows that the deblending routine erroneously selects the bump that is created by the superposition of the two lines as the secondary star’s absorption line. The secondary’s line appears as a weak contribution on the red wing of the primary star’s absorption. However, a fit with two Gaussian profiles will consider the bump as the core of the secondary line and thereby overestimate its strength.

Fig. 8 Combination of the O III λ5592 line of the primary and the secondary for model E2 at phase 0.25. The deblend/line routine will not fit the line of the secondary properly. Top: the resulting lines of the two Gaussian fits. Middle: individual lines of the primary and secondary computed with our model and normalized to the continuum of the entire system. Bottom: spectrum of the system and corresponding fit of MIDAS.

Fig. 9 Example of EW variations for the He I λ 4471 line for model E2 and for the O III λ 5592 line for model E4. The EWs have been normalized to the highest value. Black: primary. Grey: secondary. Top: variations measured on individual spectra (EWp,max,HeI = 0.17 Å, EWp,max,OIII = 0.18 Å, EWs,max,HeI = 0.34 Å, EWs,max,OIII = 0.11 Å). Bottom: variations measured on combined spectra (EWp,max,HeI = 0.17 Å, EWp,max,OIII = 0.14 Å, EWs,max,HeI = 0.10 Å, EWs,max,OIII = 0.08 Å).

Fig. 10 Example of EW variations for the He I λ 4471 and He II λ 5412 lines for model E1. The EWs have been normalized to the highest value. Black: primary. Grey: secondary. Top: variations measured on individual spectra (EWp,max,HeI = 0.04 Å, EWp,max,HeII = 1.07 Å, EWs,max,HeI = 0.59 Å, EWs,max,HeII = 0.98 Å). Bottom: variations measured on combined spectra (EWp,max,HeI = 0.05 Å, EWp,max,HeII = 0.80 Å, EWs,max,HeI = 0.21 Å, EWs,max,HeII = 0.29 Å).

4.2.1. Model E1

In this model, we mainly observed variations illustrated in Fig. 10. These variations can be explained by the orbital motion of the stars and the variation of the visible part of the stars. At periastron passage, we see the rear side of the secondary. Because the stars are close, they are quite deformed and therefore the visible part of the secondary is cool, which leads to a reinforcement of the He I lines. For the primary the hemisphere facing the companion is visible at this phase. However, because the primary is deformed, and owing to the large temperature difference between the two stars, the reflection process does not heat the primary significantly. Thus it is also the coldest part of the primary that is visible at this phase. When the binary separation increases, the stars become nearly spherical. The temperature of the primary becomes nearly constant and the EW of its lines do not change. However, for the secondary, we begin to see the front part of the star. Whilst the reflection is not very effective for heating the primary, it is very effective for the secondary and accordingly, when we see the front part of the latter, the temperature increases. The increase is stronger because the separation between the stars decreases during the second half of the orbital cycle.

The EWs measured on combined spectra (see Fig. 10, bottom) are not consistent with those measured on the individual spectra (except for the He I λ 4143 line of the secondary). As mentioned above, this inconsistency is due to the deblending routine, which does not properly separate the lines of the primary and secondary stars. For example, the blue-shifted component may have a larger EW than when it is red-shifted. This inconsistency is reduced when the lines are measured at orbital phases at which they are clearly separated.

4.2.2. Model E5

For this last model, we investigated the He I λλ 4026, 4143, 4471, 4713, 5016, Si IV λ 4089, He II λ 4200, 4542, and O III λ 5592 lines. The EWs of the He I and Si IV lines in the primary display a U-pattern with a maximum at phase 0.0–0.05 and a minimum at phase 0.45–0.55. The EWs of the secondary display a symmetric variation with respect to the primary, however, the amplitude of the variation is smaller. The O III line displays the inverse behaviour compared to the previous lines. Finally, the He II lines in the primary display two maxima at phases 0.3 and 0.7 with an intermediate minimum at phase 0.5 that give us an M-pattern variation. The secondary presents a low-amplitude U-pattern variation.

Fig. 11 Example of EW variations for He I λ 5016 and He II λ 4200 for model E5. The EWs have been normalized to the highest value. Top: variations measured on individual spectra (EWp,max,HeI = 0.20 Å, EWp,max,HeII = 0.48 Å, EWs,max,HeI = 0.21 Å, EWs,max,HeII = 0.48 Å. Bottom: variations measured on combined spectra (EWp,max,HeI = 0.12 Å, EWp,max,HeII = 0.27 Å, EWs,max,HeI = 0.11 Å, EWs,max,HeII = 0.24 Å).

The EWs measured on the combined spectra are nearly constant for the He I λλ 4026, 4471, 5016, He II λ 4200, 4542, and O III λ 5592 lines for the primary and secondary stars. The Si IV line displays irregular variations but a mean value of the primary higher than the secondary for phases before 0.5 and a mean value of the secondary higher for the second part of the orbital cycle. The He I λ 4713 line displays constant EWs for phases before and after 0.5 but the value is different before and after phase 0.5. Before phase 0.5 the primary has a lower EW than the secondary and after phase 0.5 the primary has a larger EW than the secondary. Finally, the variations observed for the He I λ 4143 line are consistent with the variations measured on individual spectra. All variations measured in our model underestimate the observed S-S effect.

In conclusion, our models display, in general, some variations of EWs during the orbital cycle. The variations on the spectra of individual components are caused by the orbital motion and the modulation of the visible part of the star. However, the variations that we measure on combined spectra and the visual variation of the relative strength of the line of the primary and secondary stars are not strong enough to completely explain the S-S effect as observed by Sana et al. (2001).

4.3. Synthetic broad-band light curves and radial velocity curve

The synthetic light curves were computed in the wavelength range 3800–7100 Å. We compared our light curves to those obtained with Nightfall. Nightfall provides synthetic light curves based on an instantaneous Roche potential and thus does not take any viscous stress into account. We point out here the similarity and difference of the two models. For models E1 and E4, our light curves present qualitatively the same characteristics as Nightfall. The TIDES code leads to more compact objects than the instantaneous Roche potential. Therefore, because the stars are more compact, all other parameters being equal, the gravity darkening is less strong and the reflection effect is more effective and consequently the stars are hotter. Near periastron, the difference between the two models is stronger because our model predicts a weak variation of 0.02 mag and Nightfall predicts no variation. The very low inclination (60°) explains this small change of magnitude. Because the stars are less deformed (than in the Nightfall model), the reflection effect is more effective near periastron and accordingly, it also increases the temperature of the stars which enhances the difference between the two models. Some other differences arise because of the time delay induced by viscous stress. For model E1, the agreement between the two models and the light curve observed for HD 93205 by Antokhina et al. (2000) is poor, however. Models E2 and E3 give very different results with our code and Nightfall. These two systems have a longer period and are more compact in our model than in Nightfall. The discrepancies probably come from the different radii of the stars. Finally, the light curve of model E5 displays the two eclipses in both models but in our model, the eclipses are deeper than evaluated with Nightfall (see Fig. 12). Again the difference can be explained by the smaller radius of the star in our model. Some other small differences can be easily explained by the viscous effects. The agreement between our model and the observations of Mayer et al. (2008) is rather good, though we slightly overestimate the depth of the eclipse of the primary.

Fig. 12 Synthetic light curves of model E5 in the range λλ 3500, 7100. The light curves have been computed with the CoMBiSpeC + TIDES models (black) and with the Nightfall model (grey).

The mean brightness ratios in the wavelength range 3800 − 7100 Å are similar to the literature values. Finally, we determined the amplitude of the radial velocity curves and found a good agreement with the observations (see Table 6) for all simulated systems.

5. Summary and perspectives

We have presented improvements of our mathematical model that allow us to compute the physical properties on the surface of massive stars in binary systems that contain main-sequence O(B) stars. The first improvement is the inclusion of the radiation pressure effect on the shape of the stars. The second improvement is the use of the TIDES code to compute the shape of the stars in eccentric and/or asynchronous systems. In both cases, we took into account various effects such as gravity darkening, reflection, and limb-darkening which allowed us to compute the temperature distribution at the stellar surface. Then we used the TLUSTY OSTAR2002 and BSTAR2006 grids (Lanz & Hubeny 2003, 2007) to compute the spectra of each star of the system as a function of orbital phase. We showed that in a number of cases, the radiation pressure does not have a strong impact on the shape of the stars and therefore on the spectra for the models studied in Paper I. This implies that the conclusions of Paper I for the S-S effect remain valid. Our results showed that the radiation pressure has a weak impact on the shape of highly deformed stars of (over)contact binary systems like model 3. However, in model 3, the small changes in the shape of the stars have resolved our previous problem of the surface temperature distribution that we failed to reproduce with our first version of the algorithm. We also studied the impact of the gravity-darkening parameter variation on the spectra of model 3. Our results indicate that, given the assumptions made for the computation, there is no need to use a more complex treatment of the gravity darkening than the classical von Zeipel theorem. The model 4 clearly showed that radiation pressure can have a strong impact on the shape of the stars. Many of the spectral lines in this model display phase-locked profile and/or strength variations. However, these variations are not sufficient to reproduce the S-S effect observed by Bonanos (2009).

The second part of this paper investigated the eccentric systems. We studied many lines and detected phase-locked profiles and/or strength variations in many cases in the individual spectra of both components. These variations are caused by the change of orientations of the stars as a function of orbital phase. They thus reflect the non-uniform temperature distribution across the stellar surface. The variations measured on the combined spectra of the systems, however, often disagree with the measurements made on individual spectra. Our results also showed that the deblending routines that fit two Gaussian or Lorentzian profiles often fail to measure the line EWs properly and lead to incorrect interpretations. The variations measured in models E2 and E5 underestimate the S-S effect observed in HD 93403 and HD 152248, which in turn might be caused by our not predict intrinsic asymmetric line profiles in these cases. In Paper I, this asymmetry was our explanation for the S-S effect: the sum of asymmetric intrinsic line profiles of individual stars are not properly deblended by Gaussian/Lorentzian profile fits. The absence of this asymmetry in eccentric systems is, however, not easy to explain and could be linked to weaker interactions between the stars in eccentric binaries than in close circular ones. Therefore, we cannot generalize the interpretation of the S-S effect of Paper I to the eccentric systems studied here. However, this highlights that analyses of binary system spectra are very difficult and some observational effects could be generated by analysing techniques and not by physical processes.

For clearly detached systems, spectral disentangling (Hadrava 1995; González & Levato 2006; Simon & Sturm 1994) helps to overcome the difficulties of the deblending routines. However, for systems with a strong temperature gradient at their surface (e.g. HD 100213, Linder et al. 2007), the technique fails because different lines have different radial velocity amplitudes. In eclipsing binaries, disentangling cannot be used for the phases near the eclipses. This technique provides mean spectra of both components, and even if we can measure variations with respect to this mean, we cannot compute the spectra at each orbital phase.

A future step could be to introduce cross-talk, which is surely substantial in eccentric and asynchronous systems. Finally, we plan to include a wind interaction zone between the stars, which could contribute to the heating of the stellar surface in two ways, either by backscattering of the photospheric photons, or by irradiation of X-ray photons emitted by the shock-heated plasma in the wind interaction zone. This requires, however, more sophisticated atmosphere models than what we have employed so far.

Acknowledgments

M.P. and G.R. acknowledge support through the XMM/INTEGRAL PRODEX contract (Belspo), from the Fonds de Recherche Scientifique (FRS/FNRS), as well as by the Communauté Française de Belgique – Action de recherche concertée – Académie Wallonie – Europe. G.K. and E.M. acknowledge support from UNAM/PAPIIT 107711.

References

All Tables

All Figures

Fig. 1 Top: comparison between the normalized spectra computed without radiation pressure (in black) and with radiation pressure (in grey) for model 1. To better distinguish spectra, the black line is wider than the grey line. Orbital phases are, from bottom to top: phases 0.1, 0.25 and 0.4. The spectra are shifted vertically by 0.25 continuum units for clarity. Middle: corresponding residuals (spectrum with radiation pressure minus spectrum without radiation pressure). The residual plots are shifted by 0.015 for clarity. Bottom: logarithm of the ratio of the fluxed spectra with radiation pressure divided by the one without. These plot are shifted by 0.006 for clarity.

In the text

	Fig. 2 Value of the radiation pressure parameter, δ, over the stellar surface for model 3. The highest values are located near the L1 point.
In the text

	Fig. 3 Surface modification by the radiation pressure for model 4. Top: star surface without radiation pressure effect. Bottom: star surface with radiation pressure effect.
In the text

	Fig. 4 Example of EW variations of the He II λ 4542 line. The EWs have been normalized to the highest value. Black: primary. Grey: secondary. Left: variations measured on individual spectra (EWp,max = 0.53 Å, EWs,max = 0.47 Å). Right: variations measured on combined spectra by fitting two Gaussians to the blended lines (EWp,max = 0.52 Å, EWs,max = 0.16 Å).
In the text

	Fig. 5 Variations of the mean, highest and lowest visible radius (left: the primary radii, right: the secondary radii). Top: variations for model E1. Bottom: variations for model E3.
In the text

	Fig. 6 Variations of the mean, highest and lowest visible radius, temperature, and log (g) for model E4. Left: values for the primary star. Right: values for the secondary star.
In the text

	Fig. 7 Variations of the mean, highest and lowest visible radius and temperature for model E5. Left: values for the primary star. Right: values for the secondary star.
In the text

Fig. 8 Combination of the O III λ5592 line of the primary and the secondary for model E2 at phase 0.25. The deblend/line routine will not fit the line of the secondary properly. Top: the resulting lines of the two Gaussian fits. Middle: individual lines of the primary and secondary computed with our model and normalized to the continuum of the entire system. Bottom: spectrum of the system and corresponding fit of MIDAS.

In the text

Fig. 9 Example of EW variations for the He I λ 4471 line for model E2 and for the O III λ 5592 line for model E4. The EWs have been normalized to the highest value. Black: primary. Grey: secondary. Top: variations measured on individual spectra (EWp,max,HeI = 0.17 Å, EWp,max,OIII = 0.18 Å, EWs,max,HeI = 0.34 Å, EWs,max,OIII = 0.11 Å). Bottom: variations measured on combined spectra (EWp,max,HeI = 0.17 Å, EWp,max,OIII = 0.14 Å, EWs,max,HeI = 0.10 Å, EWs,max,OIII = 0.08 Å).

In the text

Fig. 10 Example of EW variations for the He I λ 4471 and He II λ 5412 lines for model E1. The EWs have been normalized to the highest value. Black: primary. Grey: secondary. Top: variations measured on individual spectra (EWp,max,HeI = 0.04 Å, EWp,max,HeII = 1.07 Å, EWs,max,HeI = 0.59 Å, EWs,max,HeII = 0.98 Å). Bottom: variations measured on combined spectra (EWp,max,HeI = 0.05 Å, EWp,max,HeII = 0.80 Å, EWs,max,HeI = 0.21 Å, EWs,max,HeII = 0.29 Å).

In the text

Fig. 11 Example of EW variations for He I λ 5016 and He II λ 4200 for model E5. The EWs have been normalized to the highest value. Top: variations measured on individual spectra (EWp,max,HeI = 0.20 Å, EWp,max,HeII = 0.48 Å, EWs,max,HeI = 0.21 Å, EWs,max,HeII = 0.48 Å. Bottom: variations measured on combined spectra (EWp,max,HeI = 0.12 Å, EWp,max,HeII = 0.27 Å, EWs,max,HeI = 0.11 Å, EWs,max,HeII = 0.24 Å).

In the text

	Fig. 12 Synthetic light curves of model E5 in the range λλ 3500, 7100. The light curves have been computed with the CoMBiSpeC + TIDES models (black) and with the Nightfall model (grey).
In the text