© The Authors 2025

1. Introduction

Over the past two decades, binarity and higher multiplicity has been established as a common property of massive stars (e.g. Offner et al. 2023, and references therein). Whilst this multiplicity leads to more complex evolutionary paths (e.g. Chen et al. 2024), it nonetheless offers new opportunities for the observational determination of the physical properties of massive stars. The most obvious example is the study of eclipsing spectroscopic binaries, which provides the least model-dependent way to determine the masses and radii of stars, as it simply relies on Newton’s law of gravity and geometric models based on first principles (e.g. Mahy et al. 2020).

Furthermore, tidal interactions in binary systems can help probe the internal structure of stars. Among the most important consequences of the tidal interactions in short-period eccentric stellar binary systems is the Newtonian apsidal motion (e.g. Harmanec et al. 2014; Schmitt et al. 2016; Zasche et al. 2020; Rosu et al. 2020a; Claret et al. 2021). Apsidal motion consists of a slow secular change of the argument of periastron, ω, triggered by the non-spherical gravitational field of tidally distorted stars. The rate of this apsidal motion, $\dot{\omega }$, depends on the internal structure constant of order two (k₂; e.g. Shakura 1985). This latter quantity can be determined from theoretical stellar structure models via the resolution of the Radau differential equation (e.g. Rauw et al. 2016; Rosu et al. 2020b). Hence, measuring $\dot{\omega }$ provides an efficient way to constrain the internal structure of stars (Claret & Giménez 2010; Rosu et al. 2020b).

Tidal interactions also have other observational consequences. In close eccentric binary systems, the tides are dynamical, as they vary over the orbital cycle. This can trigger tidally excited oscillations (TEOs), provided there exists a resonance between the frequencies of the eigenmodes of the components and integer multiples of the orbital frequency (see e.g. Willems 2003, 2007; Kołaczek-Szymański & Różański 2023; Fellay & Dupret 2025). TEOs thus manifest themselves in photometric light curves as pulsational modulations of the brightness at harmonics of the orbital frequency. Moreover, the phase-dependent tidal deformation of the stars along with the phase-modulated light reflection and mutual heating effects lead to additional photometric variations such as the so-called heartbeat variation near periastron (Pablo et al. 2017; Engel et al. 2020; Kołaczek-Szymański et al. 2021).

Short-period (P_orb ≲ 20 d) eccentric binaries are thus expected to display a variety of proximity effects in their light curves and spectra. In this context, we revisit the short-period eccentric O-star binary BD+60° 497 (Rauw & De Becker 2004; Hillwig et al. 2006; Rauw & Nazé 2016a) in the young open cluster IC 1805. This non-eclipsing binary system hosts an O6.5 V((f)) primary and a late-O secondary on a 3.96 d orbit of moderate eccentricity (e ∼ 0.16; Rauw & Nazé 2016a). These properties make BD+60° 497 a promising target to search for apsidal motion and other proximity effects.

2. Observations and data reduction

2.1. Spectroscopy

Spectroscopic observations of BD+60° 497 were collected during five observing campaigns at the Observatoire de Haute Provence (OHP) in France using the 1.52 m telescope equipped with the Aurélie spectrograph (Gillet et al. 1994). These observations covered the spectral range between 4450 and 4890 Å at a resolving power of ∼10 000. For the data collected in September 2002, October 2003, and August 2018, the detector was an EEV CCD with 2048 × 1024 pixels of size 13.5 μm squared. In August 2019 and September 2022, the detector was replaced by an Andor CCD camera with 2048 × 512 pixels of the same size. The data from the 2002 and 2003 campaigns were previously described by Rauw & De Becker (2004) and used by Rauw & Nazé (2016a). Integration times of individual exposures ranged between 30 min and 45 min, depending on the weather conditions. The OHP spectra were reduced using version 21FEBpl 1.0 of the MIDAS software developed at the European Southern Observatory.

In 2014 and 2018, we collected a series of observations with the HEROS echelle spectrograph on the robotic 1.2 m TIGRE telescope (Schmitt et al. 2014; González-Pérez et al. 2022) installed at La Luz Observatory near Guanajuato (Mexico). TIGRE/HEROS spectra have a resolving power of 20 000 over the spectral range from 3760−8700 Å with a small gap around 5600 Å. Due to an instrumental issue, the data in the blue channel of the 2018 campaign could not be used. The data reduction was performed with the dedicated HEROS pipeline (Mittag et al. 2011; Schmitt et al. 2014). Telluric absorptions in the spectral regions around the He Iλ 5876, Hα and He Iλ 7065 lines were corrected by means of the telluric routine within IRAF and using the atlas of telluric lines of Hinkle et al. (2000).

The journal of the observations used in this work for spectral disentangling and radial velocity (RV) analysis is provided in Table A.1. Our spectra were continuum normalised using MIDAS, adopting the best-fit spline functions adjusted to a set of continuum windows. Preliminary RVs of both stars were determined by adjusting double Gaussian profiles to the He Iλλ 4471 and 4713 lines. These estimated RVs were then used as input values to perform spectral disentangling. For this purpose, we used our implementation of the shift-and-add spectral disentangling method described by González & Levato (2006). This method takes advantage of the Doppler shifts at different orbital phases to iteratively reconstruct the individual spectra of the components of the binary system and simultaneously improves our estimates of the individual RVs at each epoch of observation (e.g. Rauw & Nazé 2016a). At each iteration, to establish the RVs of a given component at a specific date, the Doppler-shifted reconstructed spectrum of the companion is first subtracted from the spectrum observed on that date. The residuals are then cross-correlated with a synthetic TLUSTY spectrum (Lanz & Hubeny 2003), and the updated RVs are evaluated by adjusting a parabola to the peak of the cross-correlation function (Verschueren & David 1999). For the primary and secondary stars of BD+60° 497, we used TLUSTY spectra respectively with effective temperatures T_eff = 37.5 kK and T_eff = 32.5 kK. Both template spectra assumed a surface gravity log g = 4.0 and were broadened to a projected rotational velocity of v sin i = 150 km s⁻¹.

In the blue part of the spectrum, we disentangled the spectral ranges from 4450 to 4600 Å, from 4600 to 4780 Å, and from 4780 to 4890 Å. Each of these domains contains at least one strong spectral line. We iterated our disentangling routine until the maximum absolute value of the RV correction between two consecutive iterations became less than 1.0 km s⁻¹. For most wavelength ranges convergence was achieved within fewer than 100 iterations.

2.2. Photometry

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) observed BD+60° 497 during four sectors (18 in November 2019, 58 in November 2022, and 85 and 86 from the end of October 2024 until mid December 2024). Broadband photometry over the 6000 Å to 1 μm bandpass was collected at a 2 min cadence. We retrieved the light curves, processed with the TESS pipeline (Jenkins et al. 2016), from the Mikulski Archive for Space Telescopes (MAST) portal¹. In our analysis, we used the pre-search data conditioned (PDC) photometry obtained after removing trends correlated with systematic spacecraft or instrumental effects and correcting for crowding (if any). We filtered the data to keep only those with a quality flag of zero and converted the PDC fluxes into magnitudes.

Because the pixels of the TESS CCD detectors encompass an area of 21″ × 21″ on the sky and since the aperture photometry is extracted over several pixels, we checked the Gaia data release 3 catalogue (DR3, Gaia Collaboration 2023) for possible contaminators within a 1′ radius of our target. The catalogue lists 71 objects, all of which are significantly fainter than our target. The brightest neighbouring source is 6.3 mag fainter than BD+60° 497. Hence, the TESS photometry does not suffer from significant contamination by neighbouring sources.

3. Spectroscopy

3.1. Reconstructed spectra

Spectral disentangling was applied independently to the full set of Aurélie data collected in September 2002 and October 2003 (16 spectra) on the one hand, and all Aurélie data collected after August 2018 (28 spectra) on the other hand. The reconstructed spectra are shown in Fig. 1. We found significant differences in the reconstructed spectra of the 2002−2003 and 2018−2022 epochs which are not due to differences in spectral resolution as all data were taken with the same instrumentation. These differences concern most of all the He IIλλ 4542 and 4686 lines. In 2018−2022, these lines were considerably stronger in the secondary spectrum than in 2002−2003. For completeness, we note that the disentangled blue spectra of the 2014 TIGRE/HEROS campaign are in excellent agreement with those of the 2002−2003 Aurélie campaign.

Fig. 1. Disentangled Aurélie spectra of BD+60° 497. Reconstructed spectra from the 2002−2003 epoch are shown by the blue continuous line for the primary, and the red continuous line for the secondary (shifted vertically by 0.2 units). The black dashed lines represent the result of the disentangling for the 2018−2022 data. The spectra are displayed normalised to the global continuum of the system. The broad double-minimum features around 4502, 4726, 4763, and 4780 Å are artefacts due to the disentangling of the stationary diffuse interstellar bands at those wavelengths.

Such differences might stem from mutual heating effects and a peculiar sampling of the orbital cycle. Indeed, Rauw & De Becker (2004) noted that the visibility of the secondary spectral lines changed with orbital phase. The 2002−2003 data contain mostly data taken at quadrature orbital phases, whilst 35% of the data from 2018−2022 were taken near conjunction phases. To check whether this difference in sampling biases our reconstruction, we repeated the spectral disentangling restricting ourselves to spectra with a RV separation between the two components |RV_p − RV_s|≥200 km s⁻¹. In this way, we retained 14 spectra (out of 16) for the 2002−2003 epoch and 18 spectra (out of 28) for the 2018−2022 epoch. Despite the different samplings, the reconstructed spectra in the second approach were essentially unchanged compared to our results with the full datasets of each epoch.

Figure 2 compares observed spectra (i.e. prior to any disentangling) collected near quadrature phases during the 2002−2003 and 2018−2022 epochs. One can clearly see a difference in relative strength of the lines of the two stars between those epochs. This is especially apparent for the quadrature phase when the primary star is receding. At those phases, the secondary absorptions are nearly of equal strength to those of the primary in the 2018−2022 data, whereas they were considerably weaker in 2002−2003. An alternative to the mutual heating explanation could be the blend with the spectral signature of a third component. However, Mason et al. (2009) did not detect any companion at an angular separation down to 0.03″ and with a V-band magnitude difference ΔV ≤ 3 mag in their speckle interferometric survey. Similarly, Maíz Apellániz (2010) did not detect a companion using lucky imaging with a separation down to 0.1″ and ΔV ≤ 8 mag. Gaia DR3 (Gaia Collaboration 2023) quotes a Renormalised Unit Weight Error (RUWE) value of 1.19 for BD+60° 497. This parameter expresses the quality of the adjustment of the astrometric data. If a single source model provides a good fit, RUWE is expected to be around 1.0. In our case, the RUWE value does not provide evidence of a non-single source.

Fig. 2. Comparison between the strongest lines in the Aurélie spectra as observed in 2002−2003 (blue) and 2018−2022 (red). The lines are, from left to right He Iλ 4471, He IIλ 4542, He IIλ 4686, and Hβ. The top panel illustrates the spectra observed close to the quadrature with the primary moving towards the observer, whilst the bottom panel illustrates the comparison for phases near the opposite quadrature. Spectra in the top panel were taken on HJD 2452928.641 (blue) and HJD 2459853.5509 (red), that is respectively at phases 0.99 (ω = 150°, position angle 236°) and 0.78 (ω = 261°, position angle 254°) according to our ‘RV dataset 2’ ephemerides in the forthcoming Table 1. Spectra in the bottom panel were obtained on HJD 2452527.555 (blue) and HJD 2458719.6140 (red), that is respectively phases 0.70 (ω = 144°, position angle 111°) and 0.42 (ω = 243°, position angle 130°) according to the same ephemerides.

We attempted to disentangle the whole set of Aurélie spectra accounting for the possible presence of a third star. However, as could be expected from the absence of a clear signature of a third star, these attempts did not result in a meaningful reconstruction of the spectra.

Rauw & De Becker (2004) classified BD+60° 497 as O6.5 V((f)) + O8.5-9.5 V((f)) based on observations taken around opposite quadrature phases. Sota et al. (2011) classified BD+60° 497 as O6.5 V((f)) + O8/B0 V. Rauw & Nazé (2016a) used disentangled spectra to infer an O6.5 V((f)) + O8.5 V classification. Comparison of our disentangled spectra for the 2002−2003 epoch with the standard classification spectra of Sota et al. (2011) indicates a primary spectrum consistent with an O6.5-7 V classification, except for the unusually strong Si IV 4631, C IIIλλ 4647−50−51, O IIλλ 4662−72 and O IIλλ 4699−4705 absorptions which would be more consistent with a slightly later (O7.5) spectral type. On the other hand, the disentangled spectra from the 2018−2022 epoch yield a primary spectrum with weaker H I, He I, and He II absorption lines. The He I/He II ratio remains consistent with a spectral type around O7. These results are confirmed when considering the Conti & Alschuler (1971) and Mathys (1988, 1989) criteria based on the equivalent width (EW) ratios of the He Iλ 4471 and He IIλ 4542 lines. At both epochs, we obtain ratios at the boundary between O6.5 and O7. Comparing with the grid of synthetic TLUSTY O-star spectra (Lanz & Hubeny 2003), we obtain a good match for the primary with an effective temperature of 37.5 kK and a fractional optical luminosity of 65% for the 2002−2003 epoch and of 50% for the 2018−2022 epoch.

The differences in the spectral appearance between the two epochs are more pronounced for the disentangled secondary spectrum. During the 2002−2003 epoch, the secondary spectrum was consistent with spectral type O8.5 V, but the strength of H I, He I and He II lines strongly increased in the 2018−2022 epoch. The intensity ratio of the helium lines of that epoch points towards an O7 V spectral type. Quite interestingly, the metal lines in the disentangled spectra of both stars do not undergo such strong changes. Comparing again with the grid of TLUSTY spectra, we obtain a good match for the secondary in 2002−2003 with T_eff ≃ 32.5 kK and a fractional optical luminosity of 35%. To match the 2018−2022 spectrum we need instead a model with T_eff ≃ 36.25 kK and contributing about 50% of the optical flux.

3.2. Rotational velocity and optical brightness

We applied the Fourier transform method (Simón-Díaz & Herrero 2007) to the Si IIIλ 4552 line in the disentangled primary and secondary spectra to determine the projected rotational velocity, v sin i. The best-fit rotational broadening functions correspond to v sin i = 130 ± 20 km s⁻¹ for the primary and 110 ± 20 km s⁻¹ for the secondary. The quoted uncertainties correspond to the dispersions of the measurements between the 2002−2003 and 2018−2022 epochs; the v sin i values of both stars were found to be higher in 2018−2022.

The V magnitude and B − V colour of BD+60° 497 were given as 8.773 and 0.566 by Sung & Lee (1995). This is in good agreement with the mean values of the measurements compiled by Reed (2003): V = 8.80 ± 0.03 and B − V = 0.565 ± 0.011. Adopting (B − V)₀ = −0.27 (Martins & Plez 2006), we infer a colour excess of 0.835 ± 0.011. Sung et al. (2017) derived a total-to-selective extinction ratio towards IC 1805 of R_V = 3.05 ± 0.06. This yields A_V = 2.547 ± 0.060. Maíz Apellániz & Barbá (2018) inferred a somewhat higher R_V = 3.182 ± 0.056 and thus a larger A_V = 2.692 ± 0.018, which would imply a 13% increase of the luminosity of the system compared to the Sung et al. (2017) extinction ratio. Considering nine O-type stars in IC 1805, the mean geometric distance inferred from the study of Bailer-Jones et al. (2021) based on Gaia EDR3 is (2021 ± 186) pc. For the specific case of BD+60° 497, Bailer-Jones et al. (2021) quoted a somewhat larger geometric distance $({2440}_{-114}^{+128})$ pc. Whilst the astrometric parameters of BD+60° 497 position the system outside the main locus of O-type stars in IC 1805, it still lies inside the locus of OB stars (considering B-type stars up to spectral type B3, see Appendix C) which are presumably members of IC 1805. Adopting the mean IC 1805 distance, and the extinction ratio of Sung et al. (2017), we derive an absolute magnitude of the binary M_V = −5.30 ± 0.21 that corresponds nearly exactly to the combination of the brightness of an O6.5-7 V and an O8 V star as listed by Martins et al. (2005). If we use instead the larger distance estimated for BD+60° 497, we obtain M_V = −5.71 ± 0.13, that corresponds to a 45% higher luminosity.

Comparing the EWs of the He Iλ 4471 and He IIλλ 4542, 4686 lines measured on our disentangled spectra with the EWs computed from TLUSTY synthetic spectra for stars of same spectral type, we can estimate the fractional V-band luminosities of both stars as 0.60 ± 0.10 for the primary and 0.40 ± 0.17 for the secondary. These results imply that any putative third light contribution must be of low level (less than 20% of the total light).

3.3. Orbital solution and apsidal motion

Literature RVs of our target were compiled by Mermilliod (1984)². The original data came from Hayford (1932) and Underhill (1967). With the exception of one measurement by Underhill (1967), these data did not resolve the SB2 nature of the system. This could bias the RVs. Combining these literature RVs with the measurements provided by Hillwig et al. (2006) and our own RVs obtained via spectral disentangling of the Aurélie data, for dates when the RV separation between the two stars was larger than 200 km s⁻¹, results in a set of 62 primary RV points spread over nearly a century. The vast majority of these points (57 out of 62) were acquired within the last 20 years.

To establish an orbital solution for the primary that accounts for the effects of apsidal motion, we used the method of Rauw et al. (2016) and Rosu et al. (2020a), where the primary RV data are adjusted by the mathematical expression

$$\begin{array}{c}\hfill {\mathrm{RV}}_p(t)=\gamma +K_p[cos(\varphi (t)+\omega (t))+ecos(\omega (t))],\end{array}$$(1)

where γ is the systemic velocity, K_p is the semi-amplitude of the primary RV curve, ϕ(t) is the true anomaly (computed via Kepler’s equation), e is the eccentricity, and ω(t) is the longitude of periastron, which explicitly accounts for apsidal motion:

$$\begin{array}{c}\hfill \omega (t)={\omega }_0+\dot{\omega }(t-t_0),\end{array}$$(2)

with ω₀ = ω(t₀) and t₀ as a reference time of periastron passage. The rate of apsidal motion, $\dot{\omega }$, hence becomes an explicit model parameter, along with K_p, e, the anomalistic period P_anom (that is the period between two consecutive periastron passages), ω₀ and t₀.

Our method consists then in scanning a six-dimensional parameter space and searching for the combination of parameters that best adjusts the observed RV_p(t) values. Because the secondary RVs are subject to larger uncertainties, we restricted the apsidal motion analysis to the primary RVs. After some preliminary explorations, we built a grid of models with K_p varying from 137.5 to 160.5 km s⁻¹, P_anom between 3.9587 d and 3.9614 d, t₀ between HJD 2456681.8 and HJD 2456683.8, e between 0.0 and 0.30, ω₀ between 130 and 280°, and $\dot{\omega }$ between 0 and 15° yr⁻¹. The steps were respectively 0.5 km s⁻¹, 6 × 10⁻⁵ d, 0.1 d, 0.01, 2°, and 0.15° yr⁻¹.

We applied this method both to the full set of 62 primary RVs (dataset 1), keeping in mind the caveat concerning the oldest measurements, and restricting ourselves to the 57 RVs from Hillwig et al. (2006) and our Aurélie RVs (dataset 2, see Fig. 3). In both cases, we obtained the best adjustment for K_p near 153 km s⁻¹. The best-fit parameters of the two datasets overlap, with the fit of dataset 2 being of better quality (see Table 1). Our determination of the rate of apsidal motion for dataset 2 yields a period of apsidal motion $U={59}_{-10}^{+16}$ yr. The eccentricity is very much consistent with the values previously inferred by Hillwig et al. (2006) and Rauw & Nazé (2016a).

Fig. 3. Comparison between the best-fit model for the primary RVs (black solid curve) computed with relation 1 for dataset 2. Cyan, blue, and violet filled symbols stand respectively for primary RVs from Hillwig et al. (2006) and our disentangling of the Aurélie and HEROS spectra. The orange, red, and maroon open symbols show the secondary RVs (not used in the fit) respectively from Hillwig et al. (2006), our Aurélie and HEROS data, whilst the dashed line illustrates the secondary RV solution inferred from the primary RV curve for a mass ratio q = mp/ms = 1.28 ± 0.03.

Finally, we used the RVs of the primary and secondary stars to adjust the mass ratio q = m_p/m_s = 1.28 ± 0.03. Our value of q is in excellent agreement with those of Rauw & De Becker (2004) and Hillwig et al. (2006), but smaller than the 1.56 value of Rauw & Nazé (2016a). With our new value of q, we infer minimum masses of m_p sin³i = (9.26 ± 0.90) M_⊙ and m_s sin³i = (7.23 ± 0.68) M_⊙. Comparing these minimum masses with the ‘typical’ masses taken from Martins et al. (2005) for stars in the range of spectral types inferred above, we conclude that the orbital inclination is likely close to (44.8 ± 2.5)°.

Fig. 4. Evolution of the total $\dot{\omega }$ (Newtonian plus general relativity) as a function of age. The red curve displays this evolution for solar-metallicity stars of 26.2 M⊙ and 20.4 M⊙ in a binary with e = 0.15 and Panom = 3.95984 d according to the models of Claret (2019). The solid blue line shows our best estimate of $\dot{\omega }$ for BD+60° 497, and the hatched area provides the associated uncertainty. The brown and orange curves correspond to systems respectively with stars of masses 30.3 M⊙ + 23.7 M⊙ and 21.9 M⊙ + 17.1 M⊙.

We can compare the inferred rate of apsidal motion to theoretically expected values. For this purpose, we assume the stellar rotation axes to be aligned with the axis perpendicular to the orbital plane which is seen under a 45° inclination angle. The contribution of each star to the Newtonian component of $\dot{\omega }$ depends on several parameters: the orbital eccentricity, ${\left(\frac{R_{\ast }}{a}\right)}^5$ that is the ratio of the stellar radius to the semi-major axis at the power five, and the internal structure constants k₂ (see Eqs. (5)–(7) of Rauw et al. 2016). Whilst the eccentricity is well constrained, the lack of photometric eclipses prevents us from using well-determined values for the stellar radii. The ${\left(\frac{R_{\ast }}{a}\right)}^5$ factor thus constitutes by far the biggest source of uncertainties in our calculation. Still, to get an order of magnitude estimate, we can rely on radii and internal structure constants from theoretical stellar evolution models. For this purpose, we used the tables of Claret (2019) for solar metallicity³. If we consider an age of 3.5 Myr, that is the main-sequence turn-off age of the IC 1805 cluster determined by Sung et al. (2017), the theoretical radii are R_p = 8.85 R_⊙ and R_s = 6.95 R_⊙. We hence estimated a value of ${\dot{\omega }}_N$ (corrected for the effect of rotation according to Claret 2024) near 4.90° yr⁻¹. To this, we need to add the general relativity contribution which can be evaluated from Eq. (10) in Rauw et al. (2016). This contribution amounts to 0.27° yr⁻¹. The total theoretical rate of apsidal motion for an age of 3.5 Myr thus amounts to 5.17° yr⁻¹.

Alternatively, we can compare the observational value of $\dot{\omega }$ with the theoretically expected evolution of this quantity as a function of age. For this purpose, we again used the model predictions of Claret (2019) along with the correction to k₂ for the effects of rotation following the formalism of Claret (2024). For the latter correction, we assumed that the rotational and orbital axes are aligned and are both seen under an inclination of 45°. The results are shown by the red line in Fig. 4. Comparing the theoretical values of $\dot{\omega }$ with our determination results in an age estimate of $4.{13}_{-1.37}^{+0.42}$ Myr. Sung et al. (2017) evaluated an age spread of ≤1.5 Myr about the main-sequence turn-off age of 3.5 Myr for the massive star population in IC 1805. Our evaluation is thus in reasonable agreement with the cluster age.

Propagating the uncertainties on m_p sin³i and i yields $m_p=26.2_{-4.4}^{+4.1}$ M_⊙. We repeated the above procedure adopting the extreme values of m_p and assuming an m_s consistent with the observational q value. The results are shown by the orange and brown curves in Fig. 4. The more massive the stars, the earlier the system reaches the observational value of $\dot{\omega }$. For the most extreme masses considered here, ages between 0.25 Myr and 6.25 Myr would agree within the uncertainties with our observational estimate of $\dot{\omega }$.

3.4. CoMBiSpeC simulations

An important question is whether the apparent changes in the spectral properties between the 2002−2003 and 2018−2022 epochs reflect another subtle manifestation of the apsidal motion. According to our best orbital solution, the longitude of periastron changed from around 145° in the earlier epoch to about 250° in the later epoch. Hence, in 2002−2003, periastron occured at times when the secondary’s hemisphere that faces the primary (and thus undergoes the strongest heating and reflection effects) was turned away from us. This constrasts with the more recent epoch where this hemisphere was facing the observer at periastron passage. Hence apsidal motion adds to the aspect changes of the system between the two epochs.

Fig. 5. Surface temperature distribution for the stars of the BD+60° 497 binary system as computed with CoMBiSpeC at periastron (top) and apastron (bottom). The primary and secondary input radii were taken to be respectively 9.8 R⊙ and 7.3 R⊙, and their effective temperatures are respectively 37.5 kK and 35 kK.

To quantify the mutual heating effects, we used the CoMBiSpeC code (Palate & Rauw 2012; Palate et al. 2013) to simulate the temperature distribution at the surface of both stars at different orbital phases. CoMBiSpeC discretises the surface of each star into a finite grid of elements. For each cell, the local acceleration of gravity and surface temperature are computed based on the binary gravitational potential accounting for radiation pressure, gravity darkening, and reflection effects (see Palate & Rauw 2012; Palate et al. 2013, for details). For our simulations, we adopted the following stellar parameters: T_eff, p = 37.5 kK, T_eff, s = 35 kK, M_p = 26.2 M_⊙, M_s = 20.4 M_⊙, R_p = 9.8 R_⊙ and R_s = 7.3 R_⊙. The effective temperatures are taken from Martins et al. (2005) accounting for our spectral classification, and the radii are inferred from the models of Claret (2019) for an age of 4.13 Myr that best matches the observationally determined rate of apsidal motion. The synthetic surface temperature distributions at periastron and apastron are shown in Fig. 5.

Owing to gravity darkening and mutual heating, both stars display a non-uniform temperature distribution at their surface. For the secondary the difference in temperature between the hottest and the coolest areas reaches values between 1.4 kK at apastron and 2.4 kK at periastron in our simulation. In comparison, for the primary these differences range between 0.8 kK and 2.1 kK. The temperatures along the primary’s equator reach their lowest values at periastron because, at this orbital phase, the star fills up a larger fraction of its Roche lobe. This results in a lower effective surface gravity⁴ and thus a lower effective temperature in the equatorial region. Over the hemisphere that faces the secondary, this effect is only partially compensated by the stronger irradiation from the secondary star at periastron. At this orbital phase, the hemispheres of both stars that face each other have nearly identical temperatures.

Because of the irradiation by the hotter primary, the hottest side of the secondary is always the hemisphere facing the primary. The CoMBiSpeC simulations predict a maximum temperature on this hemisphere that remains nearly constant as a function of orbital phase (near 36.1 kK). The difference between periastron and apastron only amounts to 200 K. At first sight, this result might seem counterintuitive as one would expect enhanced heating at periastron. However, as outlined above, the enhanced cooling effect of gravity darkening around periastron, lowers the relevant primary temperature. The lower primary temperature compensates the reduction in separation between the stars that would otherwise boost the mutual heating effect. This also explains why the lowest secondary surface temperatures are predicted at periastron: as a result of gravity darkening, the rear side of the secondary has its temperature varying between 33.8 kK at periastron and 34.5 kK at apastron.

We also computed the surface temperature distributions for the orbital configurations corresponding to the spectra illustrated in Fig. 2. The results are illustrated in Fig. E.1. The quadrature phases in 2002−2003 correspond to configurations where the lower limit of the secondary’s temperature distribution is lower by 300 to 500 K than at the quadrature phases of the 2018−2022 epoch. This is due to the fact that the quadrature phases in 2002−2003 happened at smaller orbital separations. Whilst this effect goes into the right direction, it seems questionable whether such a modest difference in temperature could account for the differences in the properties of the disentangled spectra. The simulations predict also a ∼500 K lower apparent temperature for the primary star at quadrature phases in 2002−2003 compared to the same phases in 2018−2022.

Fig. 6. TESS light curves of BD+60° 497. The panels illustrate, from top to bottom, the PDC photometry from Sectors 18, 58, 85, and 86.

Fig. 7. Phase-folded TESS light curves of BD+60° 497 along with our best eBEER + TEO model (blue curve, see Sect. 5.2). Phase 0.0 corresponds to periastron passage (t0 = 2456682.7, Panom = 3.95984 d, see Table 1). Conjunction with the primary star in front occurred at phases 0.030, 0.992 and 0.967 respectively in Sectors 18, 58, and 85 and 86 (red short-dashed vertical lines). Conjunction with the secondary star in front occurred at phases 0.554, 0.486 and 0.441 respectively in Sectors 18, 58, and 85 and 86 (red long-dashed vertical lines). The eBEER model was computed using the RV dataset 2 orbital solution from Table 1. We assumed an orbital inclination i = 45° and radii of 9.8 and 7.3 R⊙ for the primary and secondary stars. The albedo coefficient of both stars was 0.3. The best-fit TEOs have A2 = 4.4 mmag, ψ2 = 2.39 rad, A3 = 2.4 mmag and ψ3 = 3.39 rad.

4. Photometry

As expected from the low orbital inclination (near 45°) inferred from our RV analysis, the TESS light curves (Fig. 6) do not exhibit any strong photometric eclipses that could help constrain the rate of apsidal motion. Instead, the light curves display only a low-level modulation with a single shallow minimum occurring once per orbital cycle. Between Sectors 18 and 86, there is a progressive change in the shape of the light curve. To highlight those differences, we built a phase-folded light curve for each TESS sector. We binned the data into 200 intervals spanning 0.005 in orbital phase. Phase 0.0 was taken as the time of periastron according to the ephemerides derived from our RV study accounting for apsidal motion (RV dataset 2, Table 1). For each phase bin, the error on the average magnitude was evaluated as $\sqrt{{\sigma }_{\mathrm{disp}}^2+{\sigma }_{\mathrm{err}}^2}$, where σ_disp stands for the dispersion of the individual data points falling into that phase bin and σ_err is the mean photometric error of the data points in that phase bin. The corresponding light curves are illustrated in Fig. 7.

Fig. 8. Orbital solution and light curve of BD+60° 497 in the autumn 2022. Top panel: RV data along with our best-fit orbital solution accounting for apsidal motion (see Table 1). Primary RVs are shown in blue, secondary RVs in red. Only the primary RVs at phases with a RV separation larger than 200 km s−1 were included in the RV fit. RVs at phases when the lines are heavily blended are shown here but were not used in the fit. Bottom panel: TESS Sector 58 light curve folded with the same ephemerides. Phase 0.0 is periastron. Conjunctions occurred at phases 0.992 (primary in front) and 0.486 (secondary in front).

Although there exists some scatter in the phase-folded data points, the light curves from the different sectors are not identical as already noted above. In the Sector 18 data, the minimum is not very well defined and consists in a broad feature made of two consecutive minima of unequal depths. During Sector 58, a single minimum is observed, which appears better defined (sharper and deeper) although it clearly remains asymmetric with a slow brightness decrease followed by a faster recovery. Finally, the combined light curve of Sectors 85 and 86 again yields a sharper and deeper minimum followed immediately by a sharp maximum, which is more prominent in Sector 86 but hardly visible in Sector 85.

Fig. 9. Fourier amplitude spectra of the TESS light curves of BD+60° 497. For each sector or combination of sectors, the upper panel shows the Fourier spectrum of the photometric time series, whereas the lower panel displays the spectral window function associated with the sampling. The dashed vertical lines indicate the orbital frequency corresponding to Panom and its n = 2 and n = 3 harmonics.

The folded light curves indicate that the shallow minima do not correspond to grazing eclipses which would be expected to occur at the conjunction phase closest to periastron passage. At the epochs of the TESS observations, ω was in the range between 255° (Sector 18), 274° (Sector 58), and 286° (Sector 86). Hence, the conjunction phases closest to periastron were expected around phases 0.030 (Sector 18), 0.992 (Sector 58) and 0.967 (Sectors 85 and 86). Instead, Fig. 7 shows the shallow minima to occur around phase 0.6−0.7.

To check that there is no significant phase shift between the RV and photometric data we considered the RV data collected in August 2019 and September 2022. These data were taken within about 2 months of the photometric observations of Sectors 18 and 58. Over such a time interval, the uncertainties on P_anom would result in a phase shift of at most 0.001. Figure 8 illustrates the RV data and our best-fit orbital solution accounting for apsidal motion along with the TESS curve from Sector 58 folded with the same ephemerides. This figure clearly shows that the minimum in the light curve is offset from conjuction phases and occurs at a time when the primary moves towards the observer and the secondary moves away from the observer.

To further quantify the behaviour of the light curves, we performed a Fourier analysis. We determined the frequency content of the data of each individual sector using the discrete Fourier power spectrum method developped by Heck et al. (1985) and amended by Gosset et al. (2001). This algorithm explicitly accounts for the irregularities (gaps) in the sampling of the data. The Fourier spectrum, computed up to the Nyquist frequency of 360 d⁻¹, is essentially empty above 2 d⁻¹. Hereafter, we thus focus on its low frequency content. Whilst there is clearly some low-level intrinsic variability, the strongest signals are related to the orbital cycle.

The results are illustrated in Fig. 9 and summarised in Table 2. The dominant peak near 0.26 d⁻¹ is seen in all three Fourier spectra. The visibility of its harmonics changes with epoch. At first, in Sector 18, only the third harmonic frequency is detected. In Sector 58, the second harmonics also shows up, but the associated peak appears significantly broader than expected (full width at half maximum nearly twice the natural width), suggesting that it might correspond to a blend with another frequency. Finally, the Fourier analysis of the combined dataset from Sectors 85 and 86 reveals second and third harmonics with strongly increased amplitudes compared to previous sectors. The width of the second harmonics is now consistent with the natural width. The last column of Table 2 shows the amplitude-weighted mean of ν₁ and of 1/n of the n^th harmonic frequency. For Sectors 18 and 58, the contribution of the harmonics to this mean was restricted to n = 3, whilst n = 2 and n = 3 were considered for the combined dataset from Sectors 85 and 86. The uncertainty on the resulting fundamental frequency is taken to be Δ ν_nat/10.

Formally speaking, the fundamental frequencies found for each sector are close to – but different from – the anomalistic orbital frequency (0.2525 d⁻¹) or the sidereal orbital frequency (0.2526 d⁻¹). To check whether this could be due to genuine period variations, we performed dedicated simulations with the binary package of the Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), which is a one-dimensional stellar structure and evolution code. Using this tool, we verified if the interplay between tides and wind mass loss of both components may result in significant changes in the orbital period of BD+60° 497. The details of these computations can be found in Appendix F. Depending on the individual initial rotation rates of the components, the system may experience either positive or negative ${\dot{P}}_{\mathrm{orb}}$. However, these experiments gave $\frac{|{\dot{P}}_{\mathrm{orb}}|}{P_{\mathrm{orb}}}$ of order 10⁻⁸ to 10⁻⁷ yr⁻¹. These theoretical values thus imply $|{\dot{P}}_{\mathrm{orb}}|\le 4\times {10}^{-7}$ d yr⁻¹, several orders too small to explain the differences between fundamental photometric frequencies and orbital frequency as a result of genuine orbital period changes.

5. Discussion

There are a number of effects that can affect the photometry and spectroscopy of a system such as BD+60° 497. In a non-eclipsing eccentric binary, photometric modulations can occur because of ellipsoidal variations, mutual reflection effects and Doppler beaming (e.g. Engel et al. 2020, and references therein). In addition, TEOs with frequencies close to resonance with harmonics of the orbital frequency can also play a non-negligible role in the photometric variability (Kołaczek-Szymański & Różański 2023). Because the binary components are O-type stars, they should have strong stellar winds that interact in between the stars. This wind-wind collision can impact the light curve (e.g. Antokhin & Antokhina 2024; Kołaczek-Szymański et al. 2024). Additional effects could arise in a hierarchical triple system or from asynchronous rotation of a spotted rotator. Hereafter, we discuss the pros and cons of these various scenarios.

5.1. Proximity effects

Binary codes such as Nightfall (Wichmann 2011) and PHOEBE (Prša 2018), or the semi-analytical eccentric BEaming Ellipsoidal and Reflection formalism (eBEER, Engel et al. 2020), all predict a double-wave light curve for an eccentric non-eclipsing binary. This shape is due to the ellipsoidal deformations that result in a minimum at each conjunction phase, whilst beaming and reflection result in a single wave modulation (e.g. Nazé et al. 2023). The exact morphology of the synthetic light curve depends on the relative importance of reflection and ellipsoidal effects, which in turn depends on the assumed reflection coefficient and the adopted limb-darkening and gravity darkening coefficients (Engel et al. 2020). Among the TESS light curves of BD+60° 497 (Fig. 7), the data from Sectors 85 and 86 come closest to this description. Allowing for phase shifts and letting the codes adjust the secondary temperature, the filling factors of both stars and the value of ω, we found indeed some Nightfall and PHOEBE models that could reasonably reproduce the shape of the Sectors 85 and 86 light curve (see Fig. D.1). This model predicts two photometric minima at conjunction phases: the deepest and broadest one after apastron and a narrower and shallower minimum near periastron. However, there are several major caveats here. First of all, matching the amplitude of the observed minimum requires significantly smaller stars than expected from their spectral type. Setting instead the radii to their expected values results in a minimum twice as deep as observed. Second, matching the phase of the deepest minimum either requires a shift in phase of −0.12, which is at odds with the RV variations as shown in Sect. 4, or a value of ω around 255° instead of 286°, expected at that epoch from our apsidal motion study. Figure 10 illustrates the expected behaviour of the light curve as ω evolves from 0 to 330°. The phase of the deepest minimum and the relative height of the two maxima of the light curve vary with ω. The maxima have equal strengths for ω = 90° and ω = 270°. At ω = 255°, the phase of the predicted minimum better matches the observations and the strongest maximum is expected at phase 0.9 (see Fig. D.1). However, such a large error on the longitude of periastron would require pushing the uncertainties on ω₀ and $\dot{\omega }$ to their most extreme values. What is even worse is that the observed epoch-dependence of the light curve is at odds with the dependence on ω shown in Fig. 10, as can be seen on Fig. D.1. In Sectors 18 (resp. 58), the expected value of ω would be about 23° (resp. 9°) lower than in Sectors 85 and 86. For such values of the longitude of periastron, we expect a strong peak near phase 0.90−0.92, which clearly is not observed.

Fig. 10. Variation in the theoretical light curve computed with the Nightfall code as a function of ω. The parameters of the model were taken as i = 45°, fillp = fills = 0.6, Teff, p = 37.5 kK, Teff, s = 32.5 kK, e = 0.15.

Therefore, proximity effects alone cannot explain the shape of the observed photometric light curve. Given the properties of the system, such effects ought to exist with an amplitude larger than that of the observed variations. Hence, an important question is what could actually inhibit such effects or hide them. In the next subsections, we consider various possibilities.

5.2. TEOs

The picture could be blurred if BD+60° 497 exhibits TEOs in addition to ellipsoidal variations, reflections and beaming. Our Fourier analysis of the light curves did not unveil harmonics beyond n = 2 and n = 3, except possibly for Sectors 85 and 86, where the n = 4 harmonic is likely present at 1.0168 d⁻¹ (amplitude 0.0011 mag), although it stands out only very weakly against the noise level. The n = 2 and n = 3 harmonics fall into a frequency regime where l = 2 stellar eigenmodes should have a significant overlap integral with the forcing frequency. Hence, they could indeed excite significant TEOs that would come on top of the proximity effects.

To test this scenario, we simulated a grid of synthetic light curves including the proximity effects computed with the eBEER (Engel et al. 2020) model as well as TEOs with n = 2 and n = 3. For the proximity effects we assumed the ω values computed from our orbital solution including apsidal motion. The only free parameter associated with the eBEER model was the reflection coefficient taken equal for both stars. Each TEO brings in another two free parameters which are its amplitude A_n and phase constant ψ_n. We then attempted to adjust the light curves from all three epochs with a single set of amplitudes and phase constants using an expression of the kind

$$\begin{array}{cc}\hfill m_{\mathrm{TESS}}(t)=& \mathrm{eBEER}(\omega (t),\psi (t))+A_{2}cos(4\pi \psi (t)+{\psi }_2)\hfill \\ \hfill & +A_{3}cos(6\pi \psi (t)+{\psi }_3),\hfill \end{array}$$(3)

where ψ(t) is the orbital phase measured from periastron. We varied the albedo coefficient from 0.1 to 1.0 in steps of 0.1. This parameter appears to be anti-correlated with the amplitude A₃ of the n = 3 TEO: the lower the albedo coefficient, the higher A₃. The same rule applies to A₂ although here the relationship is less well defined. The smallest residuals were obtained for a surprisingly low albedo coefficient of 0.3. The corresponding curves are shown in Fig. 7. The agreements are clearly not perfect, but provide reasonable fits for Sector 58 and Sectors 85 and 86. However, for Sector 18 the model predicts variations between phases 0.95 and 1.35 that clearly are not observed.

We thus conclude that TEOs with n = 2 and n = 3 could alleviate some of the difficulties encountered when fitting the light curves. However, fitting all three epochs of TESS observations with the same set of TEO parameters remains challenging.

5.3. Stellar wind effects

The presence of stellar winds in the BD+60° 497 system was diagnosed thanks to low-dispersion IUE spectra that revealed prominent C IVλλ 1548, 1551 P-Cygni profiles. Burki & Llorente de Andrés (1979) estimated a terminal wind velocity v_∞ = 3600 km s⁻¹, whilst Llorente de Andrés et al. (1982) evaluated a mass-loss rate of Ṁ = 6.5 × 10⁻⁷ M_⊙ yr⁻¹. These early estimates neither accounted for the binarity of the system nor for clumps. Given their spectral types, each component of the system is likely to have a stellar wind, but we expect the wind of the hotter and more luminous primary to be the denser and more energetic one. Using the Björklund et al. (2023) recipe we estimated Ṁ_p = 1.2 × 10⁻⁷ M_⊙ yr⁻¹ for the primary, and Ṁ_s = 5.6 × 10⁻⁹ − 1.9 × 10⁻⁸ M_⊙ yr⁻¹ for the secondary depending on its effective temperature. Björklund et al. (2021) found a rather wide range of values between 2.5 and 5.5 for the ratio v_∞/v_esc where v_esc is the effective escape velocity. If we adopt a value of 4 for this ratio, we obtain asymptotic wind velocities of 3500 km s⁻¹ for the primary and 3800 km s⁻¹ for the secondary. These values are in good agreement with the observational result of Burki & Llorente de Andrés (1979).

The presence of stellar winds could impact the photometric light curves as a result of atmospheric eclipses where free electron scattering removes some part of the photons coming from the star that is located behind the wind of its companion. This effect is expected to be strongest at conjunction when the star with the stronger wind is in front. Therefore, given the orbital configuration of BD+60° 497 at the time of the TESS observations, atmospheric eclipses should be maximum near periastron. However, the minima of the light curves are not observed around those phases. Moreover, if we evaluate the typical optical depth due to electron scattering ${\tau }_0={\sigma }_T\frac{\dot{M}}{4\pi m_p{\mu }_{e}a{v}_{\infty }}$ (Eqs. (5) and (6) of Antokhin & Antokhina 2024) for the primary wind, it only amounts to 0.0003, which is clearly not sufficient to explain the observed variations. Hence, atmospheric eclipses are unlikely to play a significant role here.

Another consequence of the winds is the existence of a wind-wind interaction zone. Indeed, the stellar winds of both stars collide in between the stars. As the kinetic energy of the winds is converted into heat at the shock, this can lead to a strong phase-dependent X-ray emission on top of the intrinsic emission of each star (see e.g. Rauw & Nazé 2016b, for a review). In our case, there is no evidence of such an X-ray overluminosity: BD+60° 497 has $log\left(\frac{L_{\mathrm{X}}}{L_{\mathrm{bol}}}\right)=-7.03\pm 0.06$ following nearly exactly the canonical relation for the intrinsic X-ray emission of O-type stars (Rauw & Nazé 2016a). However, the lack of strong X-ray emission does not rule out the existence of a wind-wind interaction. Indeed, given the short orbital separation, the stellar winds will not have the time to accelerate to v_∞ before they collide. The wind-wind interaction zone should be in the radiatively cooling regime emitting mostly at longer wavelengths (EUV, UV). From our estimates of the mass-loss rates, we can estimate a wind momentum ratio $\eta =\frac{\dot{M_p}{v}_{\infty ,p}}{\dot{M_s}{v}_{\infty ,s}}\simeq 5.5-20$. Hence, the primary wind clearly dominates the wind of the secondary, and the collision between the two must be located very close to the surface of the secondary star. This could lead to a back-warming of the hemisphere of the secondary facing the apex of the wind-wind collision zone. Because of the Coriolis deflection, the shock region, and thus the part of the surface affected by the back-warming effect would be skewed by an angle, θ, with respect to the binary axis. Following Parkin & Pittard (2008), we have $tan\theta =\frac{v_{\mathrm{orb}}}{v_{\infty }}$, where v_orb is the orbital velocity of the secondary with respect to the primary and v_∞ is the terminal velocity of the slowest wind. From K_p = 153.5 km s⁻¹, q = 1.28 and i ≃ 45°, we can evaluate v_orb as 495 km s⁻¹. We thus estimated a deflection of about 8°. Such a situation would then lead to a bright spot at the secondary surface which would have its visibility modulated by the orbital motion of the stars. In such a configuration, minimum light would be reached when the bright spot is hidden from our view, that is shortly after conjunction as the secondary starts receding. This is not observed. The skewing angle should shift the minimum away from conjunction, but only by a few hundredths in phase.

A crucial question is whether the power of the wind-wind collision is sufficient to drive a significant heating of the secondary’s surface. The power radiated by the wind interaction zone comes from the conversion of the wind’s kinetic energy into heat. According to Pittard & Stevens (2002), the power radiated by the dominant wind for the case where η ≃ 10 is given by $0.033\left(\frac{1}{2}\dot{M}{v}_{\infty }^2\right)$, and the contribution of the weaker wind is about four times lower. In our case, the total power of the wind interaction zone would amount to about 2 × 10³⁴ erg s⁻¹, that is 5 L_⊙ or 2.6 × 10⁻⁵ L_bol, p. This corresponds to the power emitted in all directions, and thus it provides only an upper limit on the power radiated towards the secondary. This is several orders of magnitude less than the power of order 10⁻² L_bol, p that the secondary receives from the primary via direct illumination. Our CoMBiSpeC simulations indicated that the direct illumination of the secondary by the primary leads to an increase in effective temperature by ∼2 kK. Hence the small additional contribution due to the wind-wind interaction is unlikely to account for the ∼3 kK difference in temperature needed to explain the changes in apparent spectral type found in our spectral analysis.

5.4. Asynchronous rotational modulation

Assuming the secondary star has a non-uniform surface brightness distribution could help solve both the issue of its epoch-dependent spectroscopic signature and the problem of the shape of the light curves. Because this non-uniform brightness distribution would rotate with the star, the changing aspect of the light curve from one epoch to another could then hint at a slightly asynchronous rotation. Bright surface spots with temperature differences of a few kiloKelvin could arise from the surface emergence of magnetic fields generated in a subsurface iron convection zone (Cantiello & Braithwaite 2011). Still, to explain the changing properties of the secondary spectrum in this way, would require spots covering a significant fraction of the stellar surface.

Tidal interactions in binary systems not only tend to circularise the orbit, but also to synchronise stellar rotation and orbital motion (e.g. Hut 1981). In an eccentric binary system, such as BD+60° 497, the orbital angular velocity varies with orbital phase. Strict synchronisation can thus not be reached until the orbit has circularised. The strongest interactions take place around periastron passage, and one expects the rotational velocities to quickly synchronise with the instantaneous orbital angular velocity at periastron (Hut 1981). This corresponds to a period given by $P_{\mathrm{peri}-\mathrm{sync}}=P_{\mathrm{anom}}\frac{{(1-e)}^{3/2}}{\sqrt{1+e}}$. In our case, this period amounts to 2.8937 d. Hut (1981) introduced pseudo-synchronisation as a state where the rotation rate remains constant. The corresponding pseudo-synchronisation period is the result of averaging the tidal effects around periastron and is thus slightly longer than the instantaneous rotation period at periastron (see Eq. (45) of Hut 1981). In our case P_{pseudo − sync} = 1.205 P_{peri − sync} = 3.4875 d.

From our determinations of the projected rotational velocities, v sin i = 130 ± 20 km s⁻¹ for the primary and v sin i = 110 ± 20 km s⁻¹ for the secondary, we can estimate the rotational periods assuming the stellar radii and the inclination of the rotational axes. Regarding the rotational axes, we assume that they are aligned with the axis perpendicular to the orbital plane for which we have estimated i = 45°. For the radii, we again take 9.8 and 7.3 R_⊙ respectively for the primary and secondary stars. In this way, we estimate rotational periods of (2.70 ± 0.41) d for the primary and (2.37 ± 0.43) d for the secondary. These values suggest that BD+60° 497 has achieved near synchronisation with the instantaneous rotation period at periastron. The estimated rotational periods of the components correspond to frequencies of (0.37 ± 0.05) d⁻¹ and (0.42 ± 0.07) d⁻¹. Neither of these frequencies shows up in our Fourier periodograms of the TESS data (Fig. 9) when aliasing is accounted for. The Fourier analysis suggest that any rotational modulation would instead occur with a period not too far away from the orbital period. Hence, it seems rather difficult to explain the peculiarities of BD+60° 497 by pseudo-synchronised rotation of a spotted secondary star.

5.5. The possibility of a triple system

We consider the possibility that BD+60° 497 might be a hierarchical triple system formed by the inner close binary and a third star on a wider orbit. Such a scenario could explain some of the astrometric properties of the system. Indeed, Bailer-Jones et al. (2021) derived a somewhat larger distance from the Gaia EDR3 parallaxes than for other O-type stars of IC 1805, and Hunt & Reffert (2024) did not consider BD+60° 497 to be a member of IC 1805. If this conclusion is confirmed by future Gaia data releases, BD+60° 497 could then be part of the Cas OB6 association but not of the IC 1805 cluster. Taken at face value, the distance inferred by Bailer-Jones et al. (2021) would leave some room in the derived absolute magnitude for a third light contribution by a tertiary star (see Sect. 3.1). Moreover, the proper motion in right ascension of BD+60° 497 also deviates from that of most other O-type stars in IC 1805, though it remains in the range of OB stars of the cluster (see Appendix C). Such deviations could reflect the motion in a triple system.

We investigated whether such a scenario could explain the spectroscopic and photometric oddities found in our analysis. Concerning the epoch dependence of the reconstructed individual spectra, the spectral signatures of a third star would definitely bias the results of the disentangling procedure, which assumes the contribution of only two stars. However, there is no obvious reason why this should preferentially affect the reconstructed secondary spectrum. Given the mass ratio we derived, the contamination should impact the reconstructed spectra of both stars nearly equally.

As to the light curves, the presence of a tertiary star would not inhibit the proximity effects expected for the close binary, but the third light would dilute their signatures in the global photometry of the triple system. If the inclination of the outer orbit is close to 90°, one could also expect eclipses of the third star by the components of the inner binary. However, given the low level of the photometric variations, such events would have to be grazing eclipses. The changing configuration of the binary, due to apsidal motion, could then explain the sector-to-sector changes, but the duration of such grazing eclipses should be rather short, unlike the minima seen in the light curves which actually last at least 0.25 in phase.

In conclusion, whilst we cannot rule out the possibility of a triple system, such a scenario does not seem to easily explain the oddities of the BD+60° 497 system. The forthcoming Gaia data releases should shed new light on the relationship between this system and the IC 1805 cluster, and hopefully give further constraints on the possible presence of a tertiary component.

5.6. Comparison with other massive binaries

Finally, we briefly compare the light curve of BD+60° 497 with those of other eccentric, short-period, non-eclipsing, massive binaries. Shenar et al. (2022) present light curves for O-type binaries in the Tarantula Nebula. Among the non-eclipsing binaries in their sample, VFTS 475 (O9.7 V + B0 V, P_orb = 4.05 d, e = 0.573) and VFTS 887 (O9.7 V + O9.5 V, P_orb = 2.67 d, e = 0.056) are those that most closely ressemble BD+60° 497. Their light curves display peak-to-peak amplitudes of ≤10 mmag (VFTS 887) and ∼15 mmag (VFTS 475), which are close to the peak-to-peak amplitude of the variations seen in BD+60° 497. Unfortunately, Shenar et al. (2022) did not present a modelling of these light curves, preventing a more detailed comparison.

Among the well-studied O-type binaries of our Galaxy, HD 165052 is probably the best case to compare with BD+60° 497. HD 165052 is an O6.5 V((f)) + O7 V((f)) system with an eccentric, short-period, orbit (e = 0.098, P_orb = 2.95585 d; Rosu et al. 2023, and references therein). We analysed the TESS photometry of this non-eclipsing system (see Appendix G). As for BD+60° 497, we found variability at the 10 mmag level, part of which could not be reproduced with a binary model assuming the system parameters found in the spectroscopic study.

It thus seems that the difficulties to reproduce the light curves of short-period, non-eclipsing, eccentric, massive binaries with binary models are not restricted to the sole case of BD+60° 497. This might indicate that current models lack some ingredients in the description of such binaries. Whilst an in-depth investigation of the possible causes is beyond the scope of the present work, a promising candidate could be the effect of external radiation pressure on the shape of the stars (Palate et al. 2013). External radiation pressure, that is due to illumination by a luminous companion, acts against the gravitational pull of the companion over the illuminated part. Compared to models without radiation pressure, the drop-shaped deformation of the stellar surface facing the L₁ point is attenuated when radiation pressure is accounted for (see Fig. 3 of Palate et al. 2013).

6. Summary and conclusions

We have presented the first analysis of the RVs and TESS photometry of the close eccentric O-star binary system BD+60° 497 accounting for the presence of apsidal motion. This study enabled us to determine a total rate of apsidal motion of $\dot{\omega }=(6.{15}_{-1.65}^{+1.05})°$ yr⁻¹. When comparing with the theoretical expectations from the stellar structure models of Claret (2019), this led to an age estimate of the system of $4.{13}_{-1.37}^{+0.42}$ Myr. This estimated value is in very reasonable agreement with the age of the IC 1805 cluster (3.5 Myr) found by Sung et al. (2017).

However, our analysis also revealed a number of oddities that challenge our understanding of this system. First, we found a strange change in the reconstructed secondary spectrum of the epochs 2002−2003 and 2018−2022. When comparing with standard spectral classification criteria, the secondary had its spectral type change from O8.5 in 2002−2003 to O7 in 2018−2022. Also, the relative optical brightness of the secondary would have increased from 35% to 50%. Such a radical change is difficult to explain. It is most probably related to the apparent changes in the visibility of the secondary lines with orbital phase, which were previously observed by Rauw & De Becker (2004). This could hint at a non-uniform surface temperature distribution. We investigated several possible origins of such a feature (mutual heating, back-warming by radiation from a wind interaction zone), but it is currently unclear whether these effects can explain the observations.

The second challenging result of our study concerns the light curves, which display a single broad minimum and changes in the overall shape. Proximity effects, especially the tidal deformation of the stars, are expected to produce a light curve with two minima around conjunction phases. Not only do the observed light curves lack the second minimum, but the single broad minimum appears offset in orbital phase compared to the expectations. We investigated a number of possible explanations (action of TEOs, wind effects, magnetic spots), but none of them provides an obvious solution to the problem, though some (especially TEOs and magnetic spots) produce effects that go in the right direction.

To date, BD+60° 497 remains a challenging object for our understanding of proximity effects in massive binaries. Future spectroscopic observations will be extremely important to monitor the behaviour of the secondary spectrum and will hopefully provide a better understanding of the causes of its unusual behaviour. Such observations will further allow additional RV data to be accumulated, thereby extending the time frame over which apsidal motion can be studied. Another important aspect are the oddities of the light curve, especially since similar effects seem to be present in other non-eclipsing short-period eccentric massive binaries, such as HD 165052. This could imply that current models for proximity effects in such systems are too simplistic.

Acknowledgments

YN acknowledges support from the Fonds National de la Recherche Scientifique (Belgium). This research was supported by the University of Liège under the Special Funds for Research, IPD-STEMA Programme. This work uses data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA’s Science Mission directorate. ADS and CDS were used during this research.

References

Appendix A: Table of new RV data

Table A.1 presents the heliocentric RVs determined via spectral disentangling applied to the Aurélie and HEROS spectra of BD+60° 497. These RVs were used to determine the rate of apsidal motion. Typical errors, estimated from the dispersion of the RVs obtained via disentangling over three different wavelength domains, are 10 km s⁻¹ for the primary and 15 km s⁻¹ for the secondary star.

Appendix B: Uncertainties on the orbital solution

For a given value of K_p, the best-fit orbital solution accounting for apsidal motion was searched for in a five-dimensional parameter space by minimising the χ² between the observed and computed primary RVs. To assess the uncertainties on the various parameters, we computed the χ² values corresponding to a 1 σ and 99% confidence level. The associated iso-χ² surfaces in the five-dimensional space were then projected onto all two-dimensional planes, that is all combinations of two parameters that can be defined. Figure B.1 illustrates the corresponding contours for the K_p = 153.5 km s⁻¹ case and for the RV dataset 2. These contours allow the 1 σ uncertainties quoted in Table 1 to be established. They also highlight correlations between $\dot{\omega }$ and P_anom on the one hand and between ω₀ and t₀ on the other.

Fig. B.1. Corner plot for the five parameters (Panom, e, t0, ω0 and $\dot{\omega }$) for Kp = 153.5 km s−1 and dataset 2. The black dots indicate our best solution (see Table 1), whilst the blue and cyan contours correspond respectively to 1 σ and 99% confidence levels.

Appendix C: Astrometric parameters

Figure C.1 displays the Gaia DR3 astrometric parameters of BD+60° 497 (blue symbols) compared to those of OB stars within a radius of 30′ of the centre of IC 1805. The OB stars were selected using the TGAS catalogue⁵ (Gaia Collaboration 2016) restricted to those objects with spectral types up to B3 in the SIMBAD⁶ database. This sample includes eight O-type stars thought to belong to IC 1805 (Joshi & Sagar 1983, HD 15558, HD 15570, HD 15629, BD+60° 498, BD+60° 499, BD+60° 501, BD+60° 513, and IC 1805 188, shown by black open circles), one O9.5 star (KM Cas, open red circle) whose membership in IC 1805 was questioned by Sung et al. (2017), and 21 B-type stars with spectral types earlier than B3 (shown by green triangles). All objects with large error bars on their astrometric parameters do have a RUWE that strongly exceeds the usual limit of 1.4. For instance, the only O-star system with very big error bars is the long-period binary HD 15558 (Rauw & Nazé 2016a) for which the astrometric solution has a RUWE of 3.06.

Fig. C.1. Gaia astrometric parameters of OB stars within a radius of 30′ around the centre of the IC 1805 open cluster. The blue symbols stand for BD+60° 497. Black open circles indicate O-type star, whilst the green triangles stand for B0-B3 stars. The red circle stands for the O9.5 star KM Cas whose membership in the cluster was questioned.

Based on proper motions from the Gaia DR1 (Lindegren et al. 2016), Sung et al. (2017) concluded that the proper motions of B-type stars in IC 1805 displays a large scatter and follows an elongated distribution in the proper motion plane. With the Gaia DR3 data, the scatter is considerably reduced compared to Fig. 32 of Sung et al. (2017), though we still find a somewhat elongated distribution. BD+60° 497 lies a bit off the main locus of the O-type stars both in proper motion and parallax, but remains within the region occupied by the early B-type stars. A definite answer regarding the membership of BD+60° 497 in the IC 1805 cluster will have to await a full re-analysis of the astrometry of IC 1805 based on the forthcoming Gaia Data Release 4 (due in December 2026).

Appendix D: Nightfall fits of Sectors 85 and 86

Figure D.1 displays the TESS light curves compared to Nightfall models accounting for proximity effects and designed to adjust as well as possible the data from Sectors 85 and 86. The model parameters are q = 1.28, i = 45°, fill_p = fill_s = 0.6, T_eff, p = 37.5 kK, T_eff, s = 32.5 kK, e = 0.15, and ω = 255°. We note here that the ω value was left as a free parameter, and that its value does not agree with our apsidal motion solution. The Roche lobe filling factors correspond to mean stellar radii of 7.3 R_⊙ for the primary and 6.5 R_⊙ for the secondary, which are significantly lower than the radii from the evolutionary models of Claret (2019) – 9.8 R_⊙ for the primary and 7.3 R_⊙ for the secondary. The secondary temperature corresponds to the lower value found by comparing our disentangled spectra with synthetic spectra. This relatively low value is required to keep the amplitude of the predicted photometric minimum near periastron at a low level. In this figure, we also show the light curves for Sectors 18 and 58 with the same parameters, except for ω, which we took to be respectively 232° and 246°. The data from these sectors are clearly not well fitted.

Fig. D.1. Phase-folded TESS light curves of BD+60° 497. The blue curve shows the Nightfall model that offers a reasonable agreement with the data from Sectors 85 and 86 (bottom panel). Phase 0.0 corresponds to periastron passage. The top and middle panel illustrate the same model with the ω values propagated backwards to Sectors 18 and 58 (see text).

Appendix E: Additional CoMBiSpeC simulations

Figure E.1 illustrates the surface temperature distributions computed with CoMBiSpeC for the orbital configurations corresponding to the spectra displayed in Fig. 2.

Fig. E.1. Surface temperature distribution computed with CoMBiSpeC for the stars of the BD+60° 497 binary system at the phases corresponding to the spectra in Fig. 2. The stars are displayed as they would be seen from Earth, assuming an orbital inclination of 45°. The top left (respectively top right) panel stands for the configuration corresponding to the 2002-2003 (respectively 2018-2022) epoch spectrum in the top row of Fig. 2. The bottom left and bottom right panels likewise correspond to the two epochs in the bottom row of Fig. 2.

Appendix F: Simulations in MESA

In our simulations, we used MESA version 24.08.1, compiled with the MESA Software Development Kit (SDK version 24.7.1; Townsend 2024). The orbital parameters of BD+60° 497 and physical parameters of its components were the same as those used in our CoMBiSpeC simulations (Sect. 3.4).

We assumed that both components had identical initial chemical composition with a solar-scaled mixture of elements according to Asplund et al. (2009), and the corresponding opacity tables for this mixture. Our models are characterised by the initial mass fraction of helium, Y = 0.28, and metallicity, Z = 0.02. We assumed both stars rotate rigidly. We tested many combinations of the rotation rates of the components, experimenting with values that were in the range of 50% to 200% of the rotation rate estimated in our paper. MESA uses the mathematical formalism of Heger et al. (2000) and Heger et al. (2005) to apply structural corrections, perform different types of rotationally induced mixing and ‘diffusion’ of angular momentum between adjacent shells. We used the convective pre-mixing scheme (Paxton et al. 2019, their Sect. 5) in combination with the Ledoux criterion to define the boundaries of convective instability. Convective mixing was incorporated into the models via mixing length theory in the formalism provided by Henyey et al. (1965), with the value of the solar-calibrated mixing length parameter α_MLT = 1.82 (Choi et al. 2016). We ignored the effects of semiconvection and thermohaline mixing. The overshooting of the material above the convective core was treated in the exponential diffusion approximation developed by Herwig (2000) with an adjustable parameter, f_ov = 0.015. Mass loss rate due to the radiation-driven stellar wind has been calculated according to the prescription of Björklund et al. (2023). During the integration, we allowed for the circularisation of the orbit and synchronisation of the rotation rate of the components, following the formulae incorporated in the MESA binary package for stars with radiative envelopes (Paxton et al. 2015). For detailed information on all the above processes and their implementation in the MESA code, we refer the reader to the series of MESA instrumental papers.

Appendix G: TESS light curve of HD 165052

HD 165052 is an eccentric (e = 0.098), non-eclipsing, short-period (P_orb = 2.95585 d), O6.5 V((f)) + O7 V((f)) binary system displaying a relatively fast apsidal motion ($\dot{\omega }=11.3°$ yr⁻¹, Rosu et al. 2023). The TESS satellite observed this system in full frame mode, with a 200 s cadence, during two consecutive sectors (91 and 92, April – May 2025). We retrieved the full-frame images (FFI) data from the MAST portal and extracted the light curves with the Python software package Lightkurve. Aperture photometry was extracted on 51 × 51 pixels image cutouts. For the source mask, we adopted a flux threshold of 15 times the median absolute deviation over the median flux. The background was evaluated from those pixels of the image cutouts that were below the median flux. Background subtraction was performed by means of a principal component analysis including five components. Whilst the Gaia DR3 catalogue lists over 300 sources within a 1 arcmin radius, only one of these sources has a G-magnitude brighter than 13. Even this source is 6 magnitudes fainter than HD 165052. We thus conclude that there is no significant contamination of the TESS photometry by neighbouring sources.

The 17169 photometric data points from both sectors were combined and phase-folded using the ephemerides of Rosu et al. (2023). The resulting light curve, binned into 200 phase bins (see the grey points in Fig. G.1), displays a well-defined peak shortly before periastron. The brightness than decreases and reaches a plateau, until it drops towards a broad minimum centred on apastron. If we compute a synthetic light curve with the Nightfall code, adopting the system and stellar parameters from Rosu et al. (2023), we find some important discrepancies. First of all, we find a shift by −0.080 in phase between the synthetic curve and the observed one. Such a shift could be explained by the 5 × 10⁻⁵ d uncertainty on P_orb = 2.95585 d propagated over the 15 years between the time of periastron given by Rosu et al. (2023) and the epoch of the TESS data. We thus applied this shift to the data (Fig. G.1) and started a fitting procedure with Nightfall. Since the predicted amplitude of the variations exceeded the observed value, we let the Roche lobe filling factors, the orbital inclination and ω vary in the fitting procedure. This resulted in a very low orbital inclination (i = 11.4°), which would imply unrealistically large masses of 180 and 162 M_⊙ for the primary and secondary, respectively. The best-fit model still failed to reproduce the plateau-like part of the observed light curve between phases 0.0 and 0.35.

In a second attempt, we froze the inclination to i = 22.7° (Rosu et al. 2023), leaving only the filling factors and ω as free parameters. The resulting best-fit adjustment was of lower quality than in the previous case (see Table G.1). It failed again to reproduce the plateau in the light curve. Moreover, the best-fit Roche-lobe filling factors were smaller than those assumed by Rosu et al. (2023), which, given the (R/a)⁵ dependence of $\dot{\omega }$, would imply a lower than observed rate of apsidal motion.

In summary, we find that the observed photometric variations of HD 165052 cannot easily be explained with the binary model assuming the parameters inferred from the spectroscopic study of Rosu et al. (2023), or leaving some of them free. This situation is very much reminiscent of the problems we encountered for BD+60° 497. Unfortunately, there is no multi-epoch high-precision photometry available for HD 165052. Therefore, we cannot check whether or not the observed light curve displays variations on timescales of years as we found for BD+60° 497.

Fig. G.1. Phase-folded TESS light curve of the HD 165052 binary system as observed during sectors 91 and 92 (grey symbols), and shifted in phase by −0.08 (black symbols). The mean error per phase-binned data point is 1.6 10−4 mag. The red continuous line illustrates the synthetic curve computed with the set of parameters taken from Rosu et al. (2023). The magenta long-dashed curve shows the best-fit synthetic curve with the inclination left as a free parameter, whilst the blue short-dashed curve corresponds to the best-fit for an inclination of i = 22.7°. The model parameters are given in Table G.1.

All Tables

All Figures

Fig. 1. Disentangled Aurélie spectra of BD+60° 497. Reconstructed spectra from the 2002−2003 epoch are shown by the blue continuous line for the primary, and the red continuous line for the secondary (shifted vertically by 0.2 units). The black dashed lines represent the result of the disentangling for the 2018−2022 data. The spectra are displayed normalised to the global continuum of the system. The broad double-minimum features around 4502, 4726, 4763, and 4780 Å are artefacts due to the disentangling of the stationary diffuse interstellar bands at those wavelengths.

In the text

Fig. 2. Comparison between the strongest lines in the Aurélie spectra as observed in 2002−2003 (blue) and 2018−2022 (red). The lines are, from left to right He Iλ 4471, He IIλ 4542, He IIλ 4686, and Hβ. The top panel illustrates the spectra observed close to the quadrature with the primary moving towards the observer, whilst the bottom panel illustrates the comparison for phases near the opposite quadrature. Spectra in the top panel were taken on HJD 2452928.641 (blue) and HJD 2459853.5509 (red), that is respectively at phases 0.99 (ω = 150°, position angle 236°) and 0.78 (ω = 261°, position angle 254°) according to our ‘RV dataset 2’ ephemerides in the forthcoming Table 1. Spectra in the bottom panel were obtained on HJD 2452527.555 (blue) and HJD 2458719.6140 (red), that is respectively phases 0.70 (ω = 144°, position angle 111°) and 0.42 (ω = 243°, position angle 130°) according to the same ephemerides.

In the text

Fig. 3. Comparison between the best-fit model for the primary RVs (black solid curve) computed with relation 1 for dataset 2. Cyan, blue, and violet filled symbols stand respectively for primary RVs from Hillwig et al. (2006) and our disentangling of the Aurélie and HEROS spectra. The orange, red, and maroon open symbols show the secondary RVs (not used in the fit) respectively from Hillwig et al. (2006), our Aurélie and HEROS data, whilst the dashed line illustrates the secondary RV solution inferred from the primary RV curve for a mass ratio q = mp/ms = 1.28 ± 0.03.

In the text

Fig. 4. Evolution of the total $\dot{\omega }$ (Newtonian plus general relativity) as a function of age. The red curve displays this evolution for solar-metallicity stars of 26.2 M⊙ and 20.4 M⊙ in a binary with e = 0.15 and Panom = 3.95984 d according to the models of Claret (2019). The solid blue line shows our best estimate of $\dot{\omega }$ for BD+60° 497, and the hatched area provides the associated uncertainty. The brown and orange curves correspond to systems respectively with stars of masses 30.3 M⊙ + 23.7 M⊙ and 21.9 M⊙ + 17.1 M⊙.

In the text

	Fig. 5. Surface temperature distribution for the stars of the BD+60° 497 binary system as computed with CoMBiSpeC at periastron (top) and apastron (bottom). The primary and secondary input radii were taken to be respectively 9.8 R⊙ and 7.3 R⊙, and their effective temperatures are respectively 37.5 kK and 35 kK.
In the text

	Fig. 6. TESS light curves of BD+60° 497. The panels illustrate, from top to bottom, the PDC photometry from Sectors 18, 58, 85, and 86.
In the text

Fig. 7. Phase-folded TESS light curves of BD+60° 497 along with our best eBEER + TEO model (blue curve, see Sect. 5.2). Phase 0.0 corresponds to periastron passage (t0 = 2456682.7, Panom = 3.95984 d, see Table 1). Conjunction with the primary star in front occurred at phases 0.030, 0.992 and 0.967 respectively in Sectors 18, 58, and 85 and 86 (red short-dashed vertical lines). Conjunction with the secondary star in front occurred at phases 0.554, 0.486 and 0.441 respectively in Sectors 18, 58, and 85 and 86 (red long-dashed vertical lines). The eBEER model was computed using the RV dataset 2 orbital solution from Table 1. We assumed an orbital inclination i = 45° and radii of 9.8 and 7.3 R⊙ for the primary and secondary stars. The albedo coefficient of both stars was 0.3. The best-fit TEOs have A2 = 4.4 mmag, ψ2 = 2.39 rad, A3 = 2.4 mmag and ψ3 = 3.39 rad.

In the text

Fig. 8. Orbital solution and light curve of BD+60° 497 in the autumn 2022. Top panel: RV data along with our best-fit orbital solution accounting for apsidal motion (see Table 1). Primary RVs are shown in blue, secondary RVs in red. Only the primary RVs at phases with a RV separation larger than 200 km s−1 were included in the RV fit. RVs at phases when the lines are heavily blended are shown here but were not used in the fit. Bottom panel: TESS Sector 58 light curve folded with the same ephemerides. Phase 0.0 is periastron. Conjunctions occurred at phases 0.992 (primary in front) and 0.486 (secondary in front).

In the text

	Fig. 9. Fourier amplitude spectra of the TESS light curves of BD+60° 497. For each sector or combination of sectors, the upper panel shows the Fourier spectrum of the photometric time series, whereas the lower panel displays the spectral window function associated with the sampling. The dashed vertical lines indicate the orbital frequency corresponding to Panom and its n = 2 and n = 3 harmonics.
In the text

	Fig. 10. Variation in the theoretical light curve computed with the Nightfall code as a function of ω. The parameters of the model were taken as i = 45°, fillp = fills = 0.6, Teff, p = 37.5 kK, Teff, s = 32.5 kK, e = 0.15.
In the text

	Fig. B.1. Corner plot for the five parameters (Panom, e, t0, ω0 and $\dot{\omega }$) for Kp = 153.5 km s−1 and dataset 2. The black dots indicate our best solution (see Table 1), whilst the blue and cyan contours correspond respectively to 1 σ and 99% confidence levels.
In the text

	Fig. C.1. Gaia astrometric parameters of OB stars within a radius of 30′ around the centre of the IC 1805 open cluster. The blue symbols stand for BD+60° 497. Black open circles indicate O-type star, whilst the green triangles stand for B0-B3 stars. The red circle stands for the O9.5 star KM Cas whose membership in the cluster was questioned.
In the text

	Fig. D.1. Phase-folded TESS light curves of BD+60° 497. The blue curve shows the Nightfall model that offers a reasonable agreement with the data from Sectors 85 and 86 (bottom panel). Phase 0.0 corresponds to periastron passage. The top and middle panel illustrate the same model with the ω values propagated backwards to Sectors 18 and 58 (see text).
In the text

Fig. E.1. Surface temperature distribution computed with CoMBiSpeC for the stars of the BD+60° 497 binary system at the phases corresponding to the spectra in Fig. 2. The stars are displayed as they would be seen from Earth, assuming an orbital inclination of 45°. The top left (respectively top right) panel stands for the configuration corresponding to the 2002-2003 (respectively 2018-2022) epoch spectrum in the top row of Fig. 2. The bottom left and bottom right panels likewise correspond to the two epochs in the bottom row of Fig. 2.

In the text

Fig. G.1. Phase-folded TESS light curve of the HD 165052 binary system as observed during sectors 91 and 92 (grey symbols), and shifted in phase by −0.08 (black symbols). The mean error per phase-binned data point is 1.6 10−4 mag. The red continuous line illustrates the synthetic curve computed with the set of parameters taken from Rosu et al. (2023). The magenta long-dashed curve shows the best-fit synthetic curve with the inclination left as a free parameter, whilst the blue short-dashed curve corresponds to the best-fit for an inclination of i = 22.7°. The model parameters are given in Table G.1.

In the text