© ESO, 2012

1. Introduction

Massive stars are rare but very important objects. Their strong stellar winds interact with the ambient interstellar medium and can trigger the formation of new stars. They also play a key role in the production of various chemical elements. However, many open questions remain about these stars, regarding their formation, evolution, and physical properties. One of the most accurate ways to determine the physical properties (such as masses, temperatures, radii ...) of massive stars is to study eclipsing binary systems. However, to get yet more information, we need to improve our models and take into account effects related to the peculiarity of massive binary systems. Several authors (Sana et al. 2005; Linder et al. 2007; Linder 2008; and references therein) have presented evidence of the impact of these peculiar effects on both the spectral classification and the spectra themselves. These effects introduce complications in assigning the proper spectral types and luminosity classes. For example, main sequence stars are seen as giant stars (e.g. CPD-41°7742, Sana et al. 2005). As the spectral type of a star is one of the most important pieces of information that we can obtain from observation, it is crucial to have reliable criteria of classification for both isolated stars and binary systems. One of the most important phenomena that impact the spectral classification is the Struve-Sahade effect that consists in a variation in the relative intensity of some spectral lines of the secondary with respect to the primary star lines with orbital phase (Struve 1937). Another important effect is the surface temperature distribution that could affect the radial velocity determination. On a heavily deformed star, the temperature at the stellar surface is indeed non-uniform and the lines are not all formed at the same place. This leads to different radial velocity determinations as a function of the line that we use.

This work presents a novel way of simulating the spectra of binary systems at different phases of the orbital cycle. In a first approach, we model the case of binary systems with circular orbits containing main sequence stars. First, we use the Roche potential to compute the shape of the stars. We then apply the von Zeipel (1924) theorem and reflection effect to obtain the local temperature at the stellar surface. Finally, we use this geometrical information to compute synthetic spectra from a grid of model atmosphere spectra.

Section 2 describes the assumptions, the modelling of the geometry of the star, and the modelling of the spectra. We also present the important issues of the computation and the way in which we have treated them. Section 3 gives our first models and presents detailed analyses of the spectral classification, the Struve-Sahade effect, the light curves, and radial velocity curves. We provide a summary of our results and the future perspectives in Sect. 4.

2. Description of our method

Circular binary systems can be modelled in the following quite simple way that we can divide into two parts, one part involving the surface, gravity and temperature calculation and a second the spectra calculation. We now develop these two steps.

2.1. Geometric modelling

First of all, we need to calculate the surface of the stars of the binary system. Throughout this paper, we make the assumption of a circular orbit with both stars co-rotating with the binary motion, and acting like point-like masses. Under these assumptions, the stellar surface is an equipotential of the classical Roche potential (see e.g., Kopal 1959) (1)where x is the axis running from the centre of the primary towards the secondary, z the axis from the centre of the primary and perpendicular to the orbital plane, and y is perpendicular to x and z. The x, y, and z coordinates are dimensionless, m₁ and m₂ are respectively the masses of the primary and the secondary, a is the separation between the stars, and G is the constant of gravity.

This can be written in adimensional form using spherical coordinates centred on the primary’s centre of mass (2)where , , x = rcosϕsinθ, y = rsinϕsinθ, and z = rcosθ. We represent the stellar surface with a discretised grid of 240 × 60 points (in θ and ϕ respectively). For a given value of Ω, r(Ω,θ,ϕ) is obtained by an iterative solution of this equation using a simple Newton-Raphson technique. The knowledge of r(Ω,θ,ϕ) over the discretised stellar surface allows us to calculate the gravity. The local acceleration of gravity is given by the gradient of the Roche potential

(3)where  as above, and .

If we assume that we know the temperature at the pole, we can calculate the local temperature of each surface point following the von Zeipel (1924) theorem (4) where g = 1 in the case of massive stars.

We have explicitly made the assumption of a co-rotating system, so that the stars always present the same face to each other. This leads to a local increase in temperature owing to the reflection effect between the two stars. We follow the approach of Wilson (1990) to treat this effect. This approach consists of an iterative process in which we assume a radiative equilibrium at each point of the surface and calculate a coefficient of reflection. Once we reach convergence, we update the temperature following (5)where ℛ is the coefficient of reflection.

2.2. Spectral modelling

The second part of the algorithm computes the spectrum of the binary by summing the incremental contributions of each surface point. We use a grid of OB star synthetic spectra computed with TLUSTY (Hubeny & Lanz 1992) assuming solar metallicity to create a database of spectra (Lanz & Hubeny 20032007). Each spectrum is defined by two parameters: gravity and temperature. As we know these parameters for each point at the stellar surfaces, we can compute the local spectrum. The computation consists of a linear interpolation between the four nearest spectra in the grid. The appropriate Doppler shift is then applied to the spectrum accounting for the orbital and rotational velocity (the latter under the assumption of co-rotation) of the surface element.

We multiply the spectrum by the area of the element projected along the line of sight towards the observer. The last corrective factor applied to the local contribution to the spectrum is the limb-darkening. The limb-darkening coefficient is based on the tabulation of Al-Naimiy (1978) for a linear limb-darkening law. Finally, we sum the contribution to the total spectrum. It has to be stressed that we assume that there is no (or little) cross-talk between the different surface elements as far as the formation of the spectrum is concerned.

We emphasize that we compute the spectra at different phases, only one “hemisphere” of each star being visible for an observer at any given time. However, as the stellar surface is not a sphere, the surface gravity and effective temperature of the surface elements visible at any given phase change with orbital phase. This leads to variations in the mean visible temperature of the stars and thus variations in the line properties. Furthermore, when the primary shows its hottest part, the secondary shows its coolest. The contrast between the spectra of the two stars therefore increases. In the following, phase zero corresponds to the eclipse of the primary by the secondary. Over the first half of the orbital cycle (phase =  [0,0.5 [), the primary star has a negative radial velocity. It has to be emphasized that our simulations allow us to compute the individual spectra of each component at each phase. This information is not available for real binary systems, where disentangling can provide average individual spectra, but usually does not allow us to derive individual spectra for each phase. In this first version, we take into account neither the radiation pressure of the star that can possibly modify the shape of the atmosphere with respect to the Roche potential of very luminous stars, nor the possibility of eccentric systems. These improvements will be done in a subsequent paper.

2.3. Computational issues

For the contact binaries, some problems occur at the contact point (L1 point). If we follow the von Zeipel (1924) theorem, the temperature normally drops to zero because the gravity is null at the L1 point. However, this result is not very physical. In addition, numerical instabilities occur in the Newton-Raphson scheme for stars filling their Roche lobe at positions near L1, instabilities that affect the determination of the local radius and surface gravity. To avoid this problem, we do not follow the equipotential strictly passing through the L1 point but an equipotential passing near this point. Another problem is the lack of some synthetic spectra in the model grid. This can occur, for instance, when a high temperature is coupled with a low gravity. The algorithm treats these problems by interpolating only between the two nearest spectra in the grid and notifies the user. However, this kind of situations occurs only for very few points located near the L1 point, which have a low temperature and a very small contribution to the total stellar spectrum. These points generally have little effect on the global result. Another surprising problem is caused by reflection. In the highly deformed binary, a coupling of some points near L1 leads to a high increase in temperature. At these positions, the convergence of the reflection computation is not guaranteed. In the approach of Wilson (1990), the convergence is evaluated for each point and the iterative procedure is stopped once all the points have converged. We modify this global convergence criterion to ensure that the end of the iteration process is reached. This criterion is that at least 99% of the points converged. In some cases, despite this revised criterion, very few points have an abnormal temperature unrelated to that of the surrounding points. These points are located near the L1 point and no spectra in the grid are available for these combinations of gravity and temperature.

3. First models

The number of well-known massive O-type binary systems with a circular orbit is rather small. Thus, we selected five well-studied massive binaries to both test our algorithm and try to model the effects of binarity on the spectrum. The five systems are: HD 159176, HD 165052, HD 100213, CPD-41°7742, and AO Cas. These systems were studied by Linder et al. (2007, Paper I) for the first three, Linder (2008, Paper III) for the fifth and Sana et al. (2005, Paper II) for the fourth, respectively. For the first three systems, Paper I reports a Struve-Sahade effect, which does not appear to affect the last two systems. These systems are main-sequence systems of massive OB-stars except for AO Cas in which the secondary is a giant star. Table 1 lists the adopted parameters for each star. The parameters come from the observational analysis carried out in Papers I, II, and III. However, some systems (HD 159176, HD 165052) are not eclipsing and some parameters are unavailable. To overcome this problem, we use the calibration of Martins et al. (2005) to assign “typical” values of radius and mass. Hence, our models do not strictly correspond to the five observed systems. To avoid any confusion, the models of HD 159176, HD 165052, HD 100213, CPD-41°7742, and AO Cas are named respectively models 1, 2, 3, 4, and 5. The first step of the computation is the surface, gravity, and temperature evaluation (see Tables 2 and 3). With these results, we can compare the simulated local gravity to the classical value for a single star given by .

The most compact systems, model 1, model 2, model 4, and the primary star of model 5, are affected by small variations in their radius, log (g), and temperature. The agreement between the single star gravity and the gradient of the Roche potential is quite good. However, for model 3 and the secondary of model 5, the differences are more important. For model 3, the difference between the maximum and minimum radius is about 2.5   R_⊙ for the primary and 2.2   R_⊙ for the secondary, log (g) varies from 4.1 to 2.7 for both stars (Fig. 1 shows log (g) variations for model 3). These variations in the surface gravity lead to a change in temperature of about 8000 K for the primary star and 5000 K for the secondary star (Fig. 2 shows temperature variations for model 3). The difference between the single star log (g) and the binary model log (g) is 1.4. However, this difference appears only for the most deformed part of the star. At the poles, the difference between the single star model and binary model is negligible. For the secondary star of model 5, the difference between the maximum and minimum radius is 3   R_⊙, and log (g) varies from 3.5 to 2.6. These variations in the surface gravity lead to a change in temperature of about 8000 K. The difference between the single star log (g) and the binary system value is 1.24. This difference, again, concerns only the most deformed part of the star. Near the poles, the difference between the two models is negligible.

We also apply a disentangling method to our combined spectra with the goal to compare its output with the true spectra of the binary components computed in our simulations. For this purpose, we use the disentangling algorithm implemented by Linder et al. (2007) and based on the González & Levato (2006) approach. We then compare the results to the individual spectra computed with our method. The disentangled spectra are in very good agreement with the computed ones.

We now consider two problems encountered in observational studies of massive binaries: the spectral classification and the Struve-Sahade effect.

Fig. 1 Example of variations in log(g) for our model 3.

Fig. 2 Example of variations in temperature for our model 3.

Fig. 3 Example of normalized simulated spectra computed for our model 1. From bottom to top: phases 0 to 0.9 by step of 0.1, the spectra are shifted vertically by 0.4 continuum units for clarity.

3.1. Spectral classification

First, we apply the quantitative Conti (1971)–Mathys (19881989) 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. The spectral types that we obtain for our simulated spectra with this criterion are in good agreement with those obtained by the observational studies of the corresponding binary systems. However, the luminosity classes are very different. The criterion applied to simulated spectra identify supergiants, whilst observations find main-sequence stars (Table 4 gives the classification of the stars with the Conti-Mathys criterion). The problem seems to come from the He I λ 4143 line. This line is a singlet transition of He I. Najarro et al. (2006) reported pumping effects between He I and Fe IV lines that are not properly taken into account in model atmosphere codes and lead to incorrect He I line strengths for singlet transitions. To classify our synthetic stars, this line is therefore unsuitable and we use instead the Walborn & Fitzpatrick (1990) atlas. This atlas is more qualitative than the Conti-Mathys criterion but allows us to use a larger wavelength domain to classify the stars. The Conti-Mathys criterion was applied to the combined spectrum of the system, the individual spectra of each component, and the disentangled spectra. The disentangling is performed on 20 simulated spectra sampling the full orbital cycle. The classification based on the Walborn & Fitzpatrick (1990) atlas was performed for individual synthetic spectra. For observational data, we had to disentangle the spectra before being able to use the atlas. The atlas of Walborn & Fitzpatrick (1990) draws attention to many lines such as He I λ 4471, He I λ 4387, He II λ 4542, Si IV λ 4088, and He II λ 4200. By observing these lines (and their relative depth) and comparing our spectra to the references of the atlas, we reclassify the stars (see Table 5) achieving now a closer agreement with the classification of the real binary systems.

3.2. Struve-Sahade effect

The Struve-Sahade effect is usually defined as the apparent strengthening of the secondary spectrum when the star is approaching the observer and its weakening as it moves away (Linder et al. 2007). This effect was first reported by Struve (1937). He proposed that the effect is caused by streams of gas moving towards and obscuring the secondary. Gies et al. (1997) suggested instead that the phenomenon might be due to heating by the back-scattering of photons by the stellar-wind interaction zone. Stickland (1997) proposed that the Struve-Sahade effect could be linked to a combination of several mechanisms. Gayley (2002) and Gayley et al. (2007) argued that this effect could be linked to flows at the surface of the stars. For a more extensive discussion, we refer the reader to Linder et al. (2007). In our sample of binaries, three are known to present this effect: HD 159176, HD 165052, HD 100213 that respectively inspired our model 1, model 2, and model 3.

As a first step, we investigate the variation of some lines in the spectra at 20 different phases to look for variations in the profile. We inspect in particular the extremal phases 0.25 and 0.75, where differences in the relative strength of the lines can most easily be distinguished as the two components are the most separated at these phases. The result of this qualitative examination is given in Table 6. We must pay attention to possible blends that can modify the strength of a line leading to an apparent Struve-Sahade effect (see below).

Fig. 4 Variations in the He I λ 4026 line in model 2. From bottom to top, in dotted blue: phases from 0.1 to 0.4 in steps of 0.1. From top to bottom, in green: phases from 0.6 to 0.9 in steps of 0.1.

First of all, we detect quite some of variations in the lines. However, many of these variations can be explained by blends with weak nearby lines rather than true variations in the line intensities. We stress that, a priori, we do not expect our simulations to display a genuine Struve-Sahade effect, because the simulations are axisymmetric about the binary axis and the largest difference, which is that caused by non-uniform surface gravity and temperature, is thus expected when comparing phases near 0.0 and 0.5, rather than 0.25 and 0.75, which are the phases most concerned with the Struve-Sahade effect (see for example Fig. 4). We also performed a quantitative study on several lines to evaluate the real variations. To check the influence of a possible blend, we measure the equivalent width in the simulated combined spectra of the binary at different phases as well as in the spectra of the individual components of the binary at the same phases. The measurements of the line EWs are performed with the MIDAS software developed by ESO. The EWs of lines in the spectra of the individual stars are determined directly by simple integration with the integrate/line routine, whilst for the binary system, we use the deblend/line command as we would do for actual observations of binary system. The latter routine fits two Gaussian line profiles to the blend of the primary and secondary lines. We stress again that for real observations, we do not have access to the individual spectra of the primary and secondary at the various orbital phases. Our simulations allow us to compare the “observed” EWs (those measured by deblend/line) to the “real” ones (those measured in the individual spectra with integrate/line) and to distinguish real variations from apparent variations due to the biases of the deblending procedure. We now consider the simulations for the various systems listed in Table 1.

3.2.1. Model 1

We first study the individual spectra of each component. The lines listed in Table 6 display small variations of 1–3% over the orbital cycle. These variations are phase-locked. However, two lines vary by more than 20%: He I λ 4921 and C III λ 5696 (predicted in photospheric emission¹). Unfortunately, the He I λ 4921 line is a singlet transition of He I and is thus subject to the effect described by Najarro et al. (2006), who showed that the singlet helium transitions are, in atmosphere codes such as TLUSTY, very sensitive to the treatment of the spectrum of Fe IV. The variation in the second line is more difficult to explain. There effectively exists a small blend in the blue component, but the blended line is very weak and it seems improbable that this small blend can explain the variations of 25 − 30%. The variations in the combined spectra are more difficult to see. Furthermore, the C III line is the only one that displays large variations in the combined spectra (see Fig. 5). We stress here that this line is rather easy to measure in our simulated data: the lines of the primary and secondary are very well-separated, and the continuum is well-defined. The results are then more reliable than for other lines. The variations are phase-locked and the lines are stronger when the stars display their hotter face. However, the predicted strength of the C III lines are rather weak and would make the measurement of real data (with noise) quite difficult.

Fig. 5 EWs (in Å) of C III λ 5696 of model 1. Upper panel: EWs of primary and secondary measured in individual spectra. Lower panel: EWs of primary and secondary measured in combined spectra.

Fig. 6 EW ratio of primary to secondary for He I λ 4026 in model 2.

The primary/secondary EW ratios (all EW ratios in this paper have been measured on combined spectra with the deblend/line routine) also present variations such as these shown on Fig. 6 except for three lines. This pattern of variations appears in other models and we call it pattern 1. The He I λ 4713 line displays an inverse V variation pattern with a peak centred on phase 0.5. He I λ 4921 undergoes irregular variations with large differences between the first and the second half of the orbital cycle. Finally, C III λ 5696 shows a V-shape variation, which is expected in respect of previous results (see Fig. 5).

3.2.2. Model 2

First we study the individual spectra of each component. We can see variations in the intensities of all the lines. However, in most cases, these variations are very small of the order of 1−2% and will probably escape detection in real observations. The variations seem phase-locked, as we can see in Fig. 7, and reflect the temperature over the visible part of the stellar surface.

The analyses of the combined spectrum are affected by the same difficulties as for model 1, and are sometimes in opposition with what we observe in the individual simulated spectra: the EW of the primary line is always smaller than that of the secondary, while the combined spectrum shows the primary line to be deeper than the secondary at phases from 0.0 to 0.5 (e.g. He I λ 4026, He I λ 4471, see Fig. 7).

The He I λ 4921 line displays a more significant variation of the order of 15% (see however the remark about this line in the discussion of model 1). Moreover, in the combined spectrum, this line is blended with an O II line at λ 4924.5 so it is very difficult to determine whether the variation comes from the blend or from another effect.

Finally, we also study the ratios of the EWs of primary to secondary lines, which we find display pattern 1 variations (see Fig. 6).

Figure 8 illustrates the He I λ 4026 line at two opposite phases. Our individual spectra clearly reveal this line to be asymmetric with a steeper red wing for both stars. When fitting the combined binary spectra with two symmetric Gaussian profiles, the red component will be systematically assigned a lower flux, because the blue component will be apparently broader than the red one. Although the resulting fit can be of excellent quality, it will provide systematically larger EWs for the blue component and smaller EWs for the red. This situation therefore leads to an artificial Struve-Sahade behaviour.

Fig. 7 Same as Fig. 5 for the He I λ 4026 line in model 2.

Fig. 8 He I λ 4026 line at phases 0.25 and 0.75 in model 2. Upper panel: Primary spectra. Middle panel: secondary spectra. Lower panel: combined spectra.

3.2.3. Model 3

The analysis of model 3 is more complex than for the two other systems. The results show different behaviours. The EW ratios of the primary to secondary measured in the combined spectra display different variation patterns (see Fig. 9). The He I λ 4471 line presents pattern 1-like variations. For the first part of the orbital cycle, He II λ 4542, He I λ 4713, and He I λ 5016 lines show variations that seem phase-locked. However, the second part of the orbital cycle is either constant or irregular. Finally, the He I λ 5876 line shows phase-locked variations, which are symmetric with respect to phase 0.5. The variations have a W-shape pattern (see Fig. 9).

Fig. 9 The EW ratio of primary to secondary for several lines (from left to right, top to bottom: He I λ 4471, He II λ 4542, He I λ 4713, He I λ 5876) in model 3.

The EWs of the He I λ 4471 line for the individual spectra of the primary show a maximum at phase 0.5 when the star displays its coolest face (the front part is never visible because of the eclipse). In contrast, the secondary shows a minimum at this phase when only the hottest parts of the star are visible (the coolest are occulted by the primary). The variations are of the order of 10%. The analysis of the combined spectra is more complex. For the first half of the orbital cycle, the EWs of the primary star have a similar behaviour to the measurements in individual spectra. For the second part of the orbital cycle and for the secondary star, the variations do not appear to be phase-locked. The EWs of the He II λ 4542 line for the individual spectra display a V-shape pattern for the primary and an inverse V-shape pattern for the secondary. The variations are of the order of 20% and phase-locked. The maximum (resp. minimum) of the curves occurs when the stars present their hottest (resp. coolest) face. The combined spectra again undergoes irregular variations of the same order of magnitude as found in the study of the individual spectra. For the He I λ 4713 and He I λ 5016 lines, the proximity of other lines renders the analysis more complex because of blends. The variations seem irregular, except for the He I λ 5016 line in the individual spectra of the primary star which exhibits an inverse V-shape pattern. Finally, the EWs of the He I λ 5876 line for the individual spectra of the primary are constant. The EWs for the individual spectra of the secondary are more interesting undergoing variations of about 10% and an inverse V-shape pattern with a peak at phases near 0.5. The EWs of the combined spectra of the primary present a maximum peaked at phase 0.5 and are quite constant over the rest of the orbital cycle. The secondary shows a minimum at phase 0.5 and small variations at other phases. As we see, this system is more complex than the other two. The variations are more significant but often symmetric with respect to phase 0.5. They also seem to be phase-locked and more significant than for the other stars (for the individual spectra at least).

3.3. Broad-band light curves

The analysis of the broad-band light curves is quite interesting because we can see variations in the intensity over the orbital cycle. In the case of an eclipsing binary, the light curve is, of course, dominated by the eclipses (Fig. 11). However, even for non-eclipsing systems such as model 1, variations occur as we can see in Fig. 10. These feature are the so-called ellipsoidal variations, which are well-known. Moreover, an interesting result is the asymmetry of the curve. In the first half of the orbital cycle, the primary star is fainter than during the second half. In contrast, the secondary star is brighter in the first half and fainter in the second. Adding both light curves together results in a roughly symmetrical light curve for this particular system.

Fig. 10 Synthetic light curves of model 1 in the range λλ 3500, 7100. Upper panel: light curves of the primary and secondary stars (normalized to the maximum of emission). Lower panel: combined light curve of the binary system (normalized to the maximum of emission).

Fig. 11 Same as Fig.10, but for model 3.

Over the first half of the orbital cycle, the primary star has a negative radial velocity, and its spectral energy distribution is Doppler shifted to the blue, whilst that of the secondary is red-shifted. The behaviour of the light curves can indeed be explained by the displacement of the emission peak of the continuum. The difference between the flux emitted over the wavelength range λλ 4000, 6000 Å by an O9 star at rest and by an O9 star with a radial velocity of 100 km s^-1 is about 1%. Therefore, the shift in the maximal peak of emission towards the blue (resp. red) in the first (resp. second) half of the orbital cycle for the primary (resp. secondary) and conversely for the second half of the orbital cycle induces the asymmetries in the broad-band light curve.

3.4. Radial velocity curve

Finally, we analyse the radial velocity (RV) curve. As we know the velocity of each visible point at a given phase, we can compute the mean value of these points and thus obtain a mean radial velocity. The curves are symmetric for all the five systems and the semi-amplitudes (K) are in good agreement with the observed ones. This shows that the lines that were used in the analyses of the observations represent the true RV of the star rather well. Table 7 compares the observed and simulated semi-amplitudes of the radial velocity curves.

A non-uniform surface temperature distribution in HD 100213 was reported by Linder et al. (2007), thus we investigated in more detail our model 3 to check whether such a distribution exists in our simulation. We studied the temperature distribution by computing the amplitude of the radial velocity curves of He I and He II lines. As we can see in Table 8, the agreement between the model and the observation is not very good. The observations suggest that the He I lines are formed over the rear part of the stars (because of the larger K values). This region then has to correspond to the lowest temperatures. However, in our model, the lowest temperatures are localized in both the front and rear parts of the stars. Hence, our model shows an intermediate amplitude resulting from the weighted average of the front and rear formation regions.

For the He II lines, the observations suggest that they are formed in the intermediate and lower radial-velocity amplitude region. Thus, the lines are formed over the front and middle parts of the star. However, our model shows that the lines are only formed in the middle part of the star and not its front. We can thus conclude that there is probably an additional mechanism that enhances the reflection process and heats the front part of the stars. This mechanism could be linked to the stellar wind interaction that might warm the stellar surface.

In Sect. 2.3, we have mentioned the problem of the increase in the temperature near the L1 point in the reflection treatment. We stress that, a priori, this increase cannot be responsible for the heating of the front parts of the stars. The problem instead concerns only the convergence and affects very few points. For model 3, after 50 iterations, the model has not yet converged and the temperature of the problematic points (3 points among 14   162 for this model) were at more than 100 kK. In comparison, the convergence for the non-problematic points is reached within fewer than five iterations.

4. Summary and future perspectives

We have presented a mathematical model, based on first principles, that allows us to compute the physical properties on the surface of stars in circular-orbit massive binary systems containing main-sequence O stars. We have assumed that the stars are corotating with their orbital motion, allowing us to use the Roche potential formalism to infer the shape of the stars and the local acceleration of gravity. We have included the effects of gravity darkening and reflection prior to inferring the surface temperature distribution of the stars. These results are then combined with the TLUSTY model atmosphere code to compute synthetic spectra of each binary component and the combined binary as a function of orbital phase. Many of the spectral lines are found to display a phase-locked profile and/or intensity variations in our simulated spectra. These variations are quite small for well-detached systems, but their amplitude strongly increases for heavily deformed stars. This variability mostly reflects the non-uniform temperature distribution over the stellar surface, which is seen under different orientations as a function of orbital phase. Given the assumptions that we have made in our model, we expect these variations a priori to be symmetric with respect to the conjunction phases. Our simulations demonstrate the impact of line blending on the measurement of lines in massive binary systems. Indeed, for a large number of lines, we have found that the true variations in the EWs (which are symmetric with respect to the conjunction phases) differ from those that are measured on the combined binary spectra using conventional deblending methods and are often found to be asymmetric with respect to the conjunction phases. The most dramatic effect is seen when the intrinsic line profiles of the individual stars are asymmetric (e.g. owing to blends with weaker lines in the wings of a strong line). In this case, deblending the binary spectra by assuming some type of symmetric line profile (Gaussian or Lorentzian) introduces a systematic effect that mimics the so-called Struve-Sahade effect. This result hence leads to an alternative explanation of the Struve-Sahade effect: in some binary systems and at least for some lines, this effect could simply be an artifact caused by the fitting of blended asymmetric lines with symmetric profiles. Therefore, at least in some cases, the Struve-Sahade effect does not stem from genuine physical processes, but rather reflects a bias in the measurement of the line profiles.

Whilst deblending appears problematic, our simulations demonstrate that spectral disentangling, assuming a reasonable sampling of the orbital cycle, yields good results, which provide a good representation of the average spectrum of a binary component and can hence be used for spectral classification and model atmosphere fitting. Comparing the RV amplitudes measured for different lines in the spectrum of the contact binary system HD 100213 with the predicted RV amplitudes for this system, we found that our model apparently fails to reproduce the surface temperature distribution of this system. Some additional heating of the parts of the stars that face each other would be needed to account for the observed distribution.

In subsequent future steps, we intend to generalize our code in different ways. We plan to extend our approach to more evolved stars that have extended atmospheres. This would requires both a revision of the Roche potential formalism to account for the effects of radiation pressure (see e.g. Howarth 1997) and the usage of a model atmosphere code that includes the effects of a stellar wind. Another important step would be to adapt our code to asynchronous and/or eccentric binary systems. In these systems, the Roche potential formalism is no longer valid and must be replaced with a proper handling of tidal effects (see e.g. Moreno et al. 2011). Finally, we plan to account for the presence of a wind interaction zone between the stars. The latter could contribute to the heating of the stellar surface in two different ways, either by the backscattering of the photospheric photons or the irradiation of X-ray photons emitted by the shock-heated plasma in the wind interaction zone. These improvements will be the subject of forthcoming papers.

Acknowledgments

We 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.

References

All Tables

All Figures

	Fig. 1 Example of variations in log(g) for our model 3.
In the text

	Fig. 2 Example of variations in temperature for our model 3.
In the text

	Fig. 3 Example of normalized simulated spectra computed for our model 1. From bottom to top: phases 0 to 0.9 by step of 0.1, the spectra are shifted vertically by 0.4 continuum units for clarity.
In the text

	Fig. 4 Variations in the He I λ 4026 line in model 2. From bottom to top, in dotted blue: phases from 0.1 to 0.4 in steps of 0.1. From top to bottom, in green: phases from 0.6 to 0.9 in steps of 0.1.
In the text

	Fig. 5 EWs (in Å) of C III λ 5696 of model 1. Upper panel: EWs of primary and secondary measured in individual spectra. Lower panel: EWs of primary and secondary measured in combined spectra.
In the text

	Fig. 6 EW ratio of primary to secondary for He I λ 4026 in model 2.
In the text

	Fig. 7 Same as Fig. 5 for the He I λ 4026 line in model 2.
In the text

	Fig. 8 He I λ 4026 line at phases 0.25 and 0.75 in model 2. Upper panel: Primary spectra. Middle panel: secondary spectra. Lower panel: combined spectra.
In the text

	Fig. 9 The EW ratio of primary to secondary for several lines (from left to right, top to bottom: He I λ 4471, He II λ 4542, He I λ 4713, He I λ 5876) in model 3.
In the text

	Fig. 10 Synthetic light curves of model 1 in the range λλ 3500, 7100. Upper panel: light curves of the primary and secondary stars (normalized to the maximum of emission). Lower panel: combined light curve of the binary system (normalized to the maximum of emission).
In the text

	Fig. 11 Same as Fig.10, but for model 3.
In the text