Apsidal motion and absolute parameters of GV Nor and V881 Sco based on mid-resolution spectroscopy and TESS photometry

ABSTRACT

The veracity of stellar evolution models and theoretical internal structure constants may be subjected to stringent tests by using physical parameters obtained from spectroscopic and photometric observations of eclipsing binary systems that exhibit apsidal motion. Two eclipsing binary stars from the southern hemisphere with no prior published spectroscopic analyses were selected for this study: GV Nor and V881 Sco. Spectroscopic observations of these two eccentric binary systems were made at the South African Astronomical Observatory (SAAO) in 2018 and 2019, and their mid-resolution spectra were obtained. The radial velocities were measured using the cross-correlation method as well as by disentangling the spectra. The rotational broadening of the disentangled spectra of the components was also determined. The absolute parameters of these two systems were obtained by analyzing the SAAO-based spectroscopic data in conjunction with photometric data from Transiting Exoplanet Survey Satellite. Apsidal motion parameters for these two eccentric binary systems were calculated by analysing eclipse timing variations. Granada and Padova evolution models indicated ages of 340 ± 35 Myr for GV Nor and 760 ± 15 Myr for V881 Sco, in agreement with the observational results. The theoretical internal structure constants estimated from the Granada models were compatible with the observational values within the error limits. However, for both systems, it was found that their circularization and synchronization timescales were shorter than their respective evolutionary ages. Therefore, the fact that the components rotate faster than their synchronized values and still have rather large eccentric orbits (e ∼ 0.11) challenges present theories of circularization and spin–orbit synchronization.

1 INTRODUCTION

1.1 General

The role of detailed analyses of detached eclipsing binary systems (dEBs) in providing answers to a set of vital astrophysical questions has been discussed in many studies (e.g. Southworth 2012, and references therein). Such studies have generally adopted the assumption of circularized binary orbits to simplify the analysis. However, uniquely valuable information may be obtained from the study of eccentric binary systems. Various important questions related to long-term binary system evolution and tidal interactions in eccentric systems appear in the recent study by Kim et al. (2018), which includes evidence for a relation between the orbital period and the apsidal motion period among the systems with known apsidal motion. Furthermore, the observed period-eccentricity diagram implies the existence of an upper limit on the eccentricity that is correlated with the orbital period of a system. Such a correlation may be ascribed to the tidal circularization process, which operates differently in early-type and late-type stars. Corresponding features in the observational P − e diagram have been detected recently (Van Eylen, Winn & Albrecht 2016).

Studies of eccentric binary systems with apsidal motion can shed light on the form of the internal stellar density concentration that should be incorporated in stellar structure and evolution models. Such studies also enable tests of general relativity by measuring changes in the argument of periastron (e.g. De Laurentis et al. 2012; Baroch et al. 2021). These applications require spectroscopic and photometric observations of sufficient precision and quantity. In recent years, the procurement of Transiting Exoplanet Survey Satellite (TESS) photometric data and mid/high-resolution spectroscopic data has allowed the aforementioned two problems to be addressed by studying eccentric binary systems showing apsidal motion (e.g. Baroch et al. 2021; Claret et al. 2021). The advances in our understanding of stellar structure and evolution resulting from the Kepler and TESS missions, coupled with substantial improvements in stellar modelling resources, have raised the importance of detailed observational studies of eccentric binary systems, so that the unique windows into stellar evolution offered by these systems may be exploited. This detailed study of GV Nor and V881 Sco addresses this need.

1.2 GV Nor

GV Nor (= HD 146375, RA 16 18 41.09, DE −53 59 42.35) is an eclipsing system with an orbital period of 2.971866 d, originally discovered as a variable star by Kruytbosch (1932). Other than a rough spectral estimation of B5/8 by Houk & Cowley (1975), there have been only two brief studies of the system, producing an estimate of the period of its apsidal motion of approximately 200 yr by Zasche (2012) and Kim et al. (2018). The substantial increase in the amount of available data calls for a well-defined analysis of the entire system.

1.3 V881 Sco

V881 Sco (= HD 150384, RA 16 42 08.08, DE −33 46 00.52) was discovered as a variable star by Strohmeier & Knigge (1973). Its spectral type was classified as A1/2 IV by Houk (1982). Published studies of V881 Sco are sparse. Zasche (2012) and Kim et al. (2018) examined the apsidal motion in the system by performing O − C analyses of the times of minima. Most of the other literature sources on V881 Sco consist of catalogues.

2 SPECTROSCOPY

Spectroscopic observations of GV Nor and V881 Sco were made with the SpUpNIC [Spectrograph Upgrade: Newly Improved Cassegrain; see Crause et al. (2016), Crause et al. (2019) for details] instrument mounted at the Cassegrain focus of the 1.9-m telescope at the South African Astronomical Observatory (SAAO). We selected a grating of the spectrograph which has 1200 lines/mm, a wavelength coverage of 400–525 nm, a blaze peak at 510 nm, and a resolution of 0.06 nm (corresponding to a resolution of approximately 40 km s⁻¹ in radial velocity (RV) and an approximate resolving power R = 8000 at the wavelength of the Hβ line). We used a slit width of 1.35 arcsec for all spectroscopic observations of these two eclipsing binary systems and the associated RV standards, which were observed contemporaneously. A total of 64 and 94 spectra were obtained for GV Nor and V881 Sco, respectively, across the 2018 and 2019 observing seasons.

HR 3383 (A1V, V_r = 2.80 km s⁻¹), HD 693 (F8V, V_r = 14.81 km s⁻¹), and HR 6031 (A1V, V_r = −5.10 km s⁻¹) were observed contemporaneously as standard stars for RV measurements of the binary systems. The exposure times varied from 600 to 1200 s for each spectrum of GV Nor and V881 Sco, depending on the weather conditions. Arc spectra using a Cu/Ar lamp were taken for wavelength calibration reasons before and after each stellar image. A set of Quartz-Iodine lamp images was also taken every night for flat-field calibrations. Standard IRAF procedures were used for the spectral data reduction and calibrations.

An example of an observed spectrum for each of GV Nor and V881 Sco is shown in Fig. 1. The characteristic lines of late B-type stars are seen in the spectrum of GV Nor, whereas the characteristic lines of early A-type stars are seen in the spectrum of V881 Sco.

2.1 Cross-correlation

The cross-correlation measurements of RVs of the component stars of these two eccentric binaries were performed using the iraf package Fxcor which is based on the method of Tonry & Davis (1979) and Popper & Jeong (1994). The Mg ii (4481) line, which is the most prominent line after the hydrogen Balmer lines (see the upper panel of Fig. 1), was used for RV measurements of GV Nor. In contrast, the metallic lines (Ti ii 4468, Mg ii 4481, Fe ii / Ti ii 4550, and Fe ii 4584) in the spectral region 4450–4600 Å (see the lower panel of Fig. 1) were used for the cross-correlation measurements of V881 Sco. The derivation of RVs from Mg I triplet lines, which are prominent in the spectra of V881 Sco, was performed using Korel. The spectra of HR 6031 were also used as a template for deriving RVs of the binary components. The RVs derived from cross-correlation (and those obtained with Korel) are shown in Fig. 2.

The program Spel (Horn et al. 1994, 1996) was applied to derive orbital parameters of the two eccentric systems, using the measured RVs of the components. During the iterations, the orbital period (P) of each system was fixed at the value listed in Table 6 in Section 3.2. The RV semi-amplitudes (K₁ and K₂) of the components, and the center-of-mass velocity (V₀), eccentricity (e), periastron longitude (ω), and epoch of maximum RV (T_max) of the system were taken as the free/adjusted parameters. Here, T_max is the Heliocentric Julian Day (HJD) time at which the RV curve of the primary component is at its maximum. Preliminary orbital parameters of GV Nor and V881 Sco, together with their errors, are listed in Table 1. Throughout the text, we follow a consistent notation for the primary and secondary components, respectively. The primary component is denoted by a ‘1’, and the secondary by a ‘2’, for all the parameters. However, this notation does not necessarily imply that the primary component is the more massive star, as discussed later in this section (for clarification: the primary is the component that is eclipsed when the phase is 0.0).

2.2 Disentangling of spectra

We ascertained the orbital parameters from the spectroscopy and separated the spectra of the component stars by using the Korel Fourier disentangling code (Hadrava 1995), which incorporates the periastron advance. We used the code to fit the line strengths of both components in each exposure (using the s-factors,¹ see Hadrava 1997, for the mathematical details). In the case of GV Nor and V881 Sco, we did not find any systematic phase-dependent variations for eclipses and proximity effects; the only (slight) variations that we detected were clearly caused by the unevenness of the chosen continua with respect to which the spectra were normalized. The sub-pixel precision of the code (Hadrava 2009) coupled with an oversampling of the input spectra mitigated the low resolution of the observed spectra.

We disentangled three spectral regions around the Balmer lines Hβ and Hγ, and Mg ii 4481 Å, respectively, for GV Nor. The latter spectral region also contains the He I line 4471 Å; however, this line was not very pronounced in any of the component spectra. The Hδ line did not have a sufficiently well-covered continuum behind its short-wavelength wing, so it had to be excluded from the analysis as well. The disentangled spectra revealed that the line profiles of the components were very similar, although the Mg ii line was stronger in the secondary component. The orbital parameters obtained from the disentangled spectra are displayed in Table 2. P_a is the anomalistic period, T_peri the time of periastron passage, |$\dot{\omega }$| the periastron advance, and q ≡ M₂/M₁ = K₁/K2 the mass ratio.

The solutions for individual lines were not mutually consistent and the errors of the mean values of the parameters were several times larger than the errors of individual sets of parameters. These discrepancies arose from the limited instrumental resolution and quality of the spectra, as confirmed by the outliers of the RVs (computed using Korel) revealed by the fit of individual exposures as a superposition of the disentangled line profiles. In a few instances, the outliers occurred in successive spectra obtained on the same night; it is these outliers that caused the large (especially in the Hγ region) r.m.s. errors ΔRV₁ and ΔRV₂ resulting from the differences (O − C)_{1, 2} between the disentangled RVs and the RV curve corresponding to the disentangled parameters (see Table 2 and Table A1 of Appendix  A).

The Korel RVs derived from lines in the Hβ region, together with the cross-correlation RVs derived from lines in the Mg ii 4481 Å region, are shown in the upper panel of Fig. 2. As shown in Table 2, the most precise values are obtained from the disentangling of the Hβ line. We provide the RVs found from this solution in Table A1 of Appendix  A. Notwithstanding the optimal precision of the fit derived from the Hβ lines, even in the Hβ region the precision of the semi-amplitudes of the RV-curves is not sufficient to establish which component is the more massive, whereas the fits obtained from the other two spectral regions (with lower precision) deliver contradictory results.

In order to fit the disentangled spectra of GV Nor, we selected models with log g  = 4. The best fit was obtained for T₁ = 12 000 K and T₂ = 12 250 K, rotational broadening of the first component equal to 180 km s⁻¹, systemic velocity equal to 45 km s⁻¹, and the mean damping of the line strengths equal to 0.47. Similarly, the best fit for the secondary component corresponded to a rotational broadening of 194 km s⁻¹, systemic velocity of 38 km s⁻¹, and mean damping of 0.63. The fits of these synthetic spectra to the observed Hβ line are shown in Fig. 3. However, these values should be treated with caution because of the relatively low resolution (40 km s⁻¹ per pixel; see Section 2.1) of the observed spectra. The systemic velocities derived from the spectra of both components also coincide within the precision of the measurements; consequently, there is no reliable evidence of any effects (for example, a stellar wind) that would cause these values to differ.

For V881 Sco, considering the relatively short time interval during which the available spectra were obtained (362 d), we fixed the orbital period and the rate of periastron advance to the values found from photometry, namely, the anomalistic period P_a = 2.491665304 d and |$\dot{\omega }=0.00559488^\circ$|/day. For final computations we also fixed the eccentricity e = 0.1068 and periastron longitude ω = 348.37468° (at the epoch HJD ≃ 2452126.5 – it should be noted that ω ≃ 21.5° at the epoch 2458506.3 in the middle of the spectroscopic observations), in accordance with the solution of the photometry. In the disentangling process, we converged only the periastron epoch T_peri, the semi-amplitude K₁ of the RV curve of the primary component, and the mass ratio q.

Notwithstanding the presence of many narrow metallic lines in the spectra of V881 Sco (see the lower panel of Fig. 1), which are usually the best tracers of RV, these lines did not contribute reliable results in our analysis. Taking into account the typical width of 3–4 Å of these lines and the width of 0.6 Å of a pixel (which corresponds to about 40 km s⁻¹), the narrow lines were covered by only 5–7 pixels of the detector. Consequently, the profiles of these lines could not provide precise radial velocities unless their depth significantly exceeded the noise. However, because the line Mg ii 4481 Å was shallower in the spectra of V881 Sco, for this system we also disentangled the spectral region containing the Mg I triplet 5167.3 Å, 5172.7 Å, and 5183.6 Å, in addition to the three spectral regions – centred on the Hβ, Hγ, and Mg ii 4481 Å lines – that were analysed for GV Nor. Each region was sampled in 2048 bins in RV; the step sizes are shown in Table 3, together with the wavelength ranges, the determined values of the parameters (and Bayesian estimates of their errors), the residual noise σ of the disentangled spectra, and the r.m.s. errors ΔRV_{1, 2} of the radial velocities of both components. The disentangled line-profiles of the Hβ line in V881 Sco are given in Fig. 4. The Korel RVs for the Hβ line region, together with the cross-correlation RVs from the metallic lines in the spectral region 4450–4600 Å, are shown in the lower panel of Fig. 2. As presented in Table 3, the most precise values were obtained from the disentangling of the Hβ line. Therefore, we provide the RVs found from this solution in Table A2 of Appendix  A. This solution was improved by also converging the eccentricity and periastron longitude (resulting in final values of e = 0.112 and ω = 345.8°); however, the improvement was insignificant, and we concluded that the values derived from the photometry were preferable.

The results derived from the regions containing the Hγ and Mg ii 4481 Å lines were associated with larger errors than the results derived from the Hβ line. However, the Mg ii 4481 Å line is a useful indicator of the temperature. In the disentangled spectra of both components, the lines of Ti II 4443 and 4468 Å appeared instead of He I 4471 Å, which would dominate at higher temperatures. The relative strengths of these lines are highly sensitive to the temperature, which can be determined by comparison with a set of synthetic spectra computed for different parameters of the stellar atmospheres. Such a comparison can either be conducted as a subsequent step in the analysis of the disentangled spectra of the component stars, as described in the next subsection, or performed directly using a template-constrained disentangling procedure. In both cases it is necessary to take into account that the disentangled spectra are normalized with respect to the overall continuum of the system; consequently, the depths of the lines are diminished compared to the model spectra of single stars, by an undetermined factor. In order to fit the observed spectra with the model, it is also necessary to establish the appropriate degree of broadening of the lines (rotational and/or instrumental) and the appropriate Doppler shift (the systemic velocity).

For the constrained disentangling procedure, we used a set of synthetic spectra computed using the code Spectrum v2.75 by Gray (1999) by referring to the solar-abundance LTE model atmospheres for a range of values of T_eff and log g published by Castelli & Kurucz (2003). The results are shown in Fig. 5. For each component star we disentangled the line broadening and Doppler shift by using variable line-strength factors for the synthetic templates, whereas the spectrum of the other component was fixed by the template derived from the free disentangling. It can be seen that the shallow features observed in the spectra are due to the blending of many weak lines. The minimal residual noise (with r.m.s. σ = 0.0088 of the continuum level) was found for the synthetic spectrum of the primary component that corresponded to T_eff = 8250 K and log g = 5.0. The estimated line broadening of the primary was 131 km s⁻¹ and the systemic velocity was 30 km s⁻¹. A quadratic fit to the residuals from neighbouring models produced the following estimated values of the parameters of the optimal model atmosphere: T_eff = 8330 K and log g = 5.07. Similarly, the best fit (with σ = 0.0083) of the secondary component was found for the synthetic spectrum with T_eff = 9250 K and log g = 4.5, rotationally broadened by 135 km s⁻¹ and Doppler-shifted by 29 km s⁻¹. The quadratic fit to the residuals from neighbouring models produced an estimated temperature for the secondary of 9225 K and log g = 4.42. In accordance with the low resolution of the spectrograph, it is clear that a non-negligible part of the broadening resulted from the instrumental profile.

2.3 Rotational velocities

We used the line profile analysis program Prof (Oláh et al. 1992; Budding & Zeilik 1994) to investigate rotational effects in these two eccentric systems. Prof convolves Gaussian and rotational broadening, and characterizes a line profile in terms of the following parameters: the local continuum intensity I_c, central line depth I_d at mean wavelength λ_m, rotational broadening parameter r, and Gaussian broadening parameter s for a given line, as well as a limb-darkening coefficient u, which we kept fixed at u = 0.5 in our analysis.

We fitted the Korel disentangled profiles of the Mg ii 4481 line for the components of GV Nor with Prof. The results are listed in Table 4 and are illustrated in Fig. 6. Using the values of the rotational broadening parameter r listed in Table 4, the derived mean observed projected rotation velocities of the primary and secondary components of GV Nor were identical within the errors: 133 ± 14 and 135 ± 13 km s⁻¹, respectively. The predicted rotation velocities obtained by using the orbital period and the derived radii of the component stars, assuming synchronization of both components’ rotation rates with the orbital period, were derived as 42 and 57 km s⁻¹, respectively. The discrepancy between these values and the rotation rates inferred from rotational broadening suggests that the components do not rotate synchronously with their orbital motion. The discrepancy may also arise from the wide instrumental line profile and should therefore be verified by spectroscopic observations obtained with a higher resolution. The values of the Gaussian broadening parameter s were derived as 38 ± 8 and 35 ± 6 km s⁻¹ for the primary and secondary components of GV Nor, respectively. This parameter is generally related to the thermal broadening and surface turbulence. Considering the effective temperatures of the components and the elements whose spectral lines are fitted here, thermal velocities were expected to be of the order of 10 km s⁻¹. Accordingly, turbulent effects on the surfaces of the component stars appear to be the dominant contributor to the Gaussian broadening.

We fitted three Korel disentangled lines (Ti ii 4443, Ti ii 4468, and Mg ii 4481) to investigate the rotational effects in V881 Sco. The profile fitting results for the disentangled lines of Mg ii 4481 are listed in Table 4 and are displayed in Fig. 7. Using the values of the rotational broadening parameter r listed in Table 4, the mean observed projected rotation velocities of the primary and secondary components of V881 Sco were derived as 125 ± 17 and 110 ± 16 km s⁻¹, respectively. Conversely, using the orbital period and the measured radii of the component stars, the synchronous rotation velocities of the primary and secondary components were derived as 51 and 47 km s⁻¹, respectively. Therefore, these values indicate (within the same observational uncertainty) that the components of V881 Sco also do not rotate synchronously with their orbital motion. The values of the Gaussian broadening parameter s in Table 4 were derived as 28 ± 8 and 22 ± 6 km s⁻¹ for the primary and secondary components of V881 Sco, respectively. These values suggest that the turbulence effect is dominant on the surfaces of the components of V881 Sco, in accordance with the case for GV Nor.

3 PHOTOMETRY

3.1 Transiting exoplanet survey satellite photometry

GV Nor (TIC 208301131) was observed by the TESS with 30-min cadence in sector 12 (2019 May 21st to June 18th) and with 10-min cadence in sector 39 (2021 May 27th to June 24th). V881 Sco was observed by TESS with 24-min cadence in sector 12 (2019 May 21st to June 18th). Photometry of these two eccentric systems is publicly accessible (Ricker et al. 2015). We downloaded and processed these data using the Lightkurve Python package (Lightkurve Collaboration 2018).

For the analysis of the TESS data, we used the numerical integration code (wd) of Wilson & Devinney (1971), combined with the Monte Carlo (MC) optimization procedure discussed by Zola et al. (2004). This method models the light curve of a given binary star, considering the ellipticity and proximity effects, and the equipotential surfaces of the components. The most important contribution of the MC search procedure in the WD+MC method of analysis is the use of a range of input values for the free parameters, rather than fixed input values. Thus, the WD + MC method searches for the model that provides the best fit of the observational light curve in the solution space, by conducting hundreds of thousands of iterations in the selected input ranges.

The WD + MC method was applied in three steps: (i) The TESS light curve was resolved by using the parameters and input ranges defined below; we used the mass ratio determined by Korel as an input value in this step. (ii) The RV curve was resolved by incorporating the tidal and rotational effects (in the form of the ellipticity and proximity effects) determined in the first step, generating a corrected value of the mass ratio, q_corr, in the process. Since the most precise values were obtained from the disentangling of the Hβ line, as shown in Tables 1, 2 and 3; the Korel RVs measured from Hβ lines of GV Nor and V881 Sco were solved during this step. (iii) In the last step, the procedure in step (i) was repeated, using q_corr as the new input value for the mass ratio.

The following procedure was followed to estimate the effective temperatures of the components: First, the spectral type of B5/8 for GV Nor allocated by Houk & Cowley (1975) was considered (see Section 1.2). Subsequently, we studied the observed spectra of GV Nor at phases close to the conjunction of the components. The spectral types of the components were estimated as B8V + B9V using the line matching method for the Hγ, Hβ, He I 4471, and Mg ii 4481 lines in the observed spectra of the components and the sample spectra of late B-type main sequence stars in Gray & Corbally (2008)’s book. The equivalent widths (EWs) of the He I 4471 line were measured with the program Prof. The means of the measured EWs for the primary and secondary components were derived as 0.15 ± 0.02 and 0.13 ± 0.02, respectively. The temperatures of the components were both estimated to be ∼11 000 K, by comparison with the calibrated data of Leone & Lanzafare (1998). Furthermore, the mass values of the components (M₁≅2.6 and M₂≅2.8 M_⊙) obtained from the solution of the RV curves in Table 1 imply spectral types close to B9V for both components, assuming that the components are main sequence stars. Based on these results, we initially assigned a temperature of 10700 ± 600 K, corresponding to a spectral type of B9V (as taken from the calibrated data of Pecaut & Mamajek (2013)) to both components. The preliminary solution of the TESS light curve indicated that the light contribution of the secondary component was dominant. Therefore, the temperature of the secondary component (T₂) was fixed at 10 700 K, and the temperature of the primary component (T₁) was adjusted in the range of 9000–12 000 K.

Three different methods were used to determine the effective temperatures of the components of V881 Sco. (i) The spectral types of the components were estimated as A1V + A2V using the line matching method for the Hγ, Hβ, Ca I 4226, Fe I (4271, 4383), Fe ii (4233, 4550, 4584, 4923), Ti ii (4468, 4550), Mg ii 4481, and Mg I (5167 - 5183) triplet lines (see Fig. 1 for V881 Sco) in the spectra observed at the conjunction phases of the components and sample spectra for early A-type main sequence stars in Gray & Corbally (2008)’s book. (ii) The solution of the RV curves obtained in this study implied that the masses of the components are approximately 2 M _⊙, which indicates that their spectral types are close to A0. (iii) As can be seen in Section 2.2, the code Spectrum produced T₁ = 8330 K and T₂ = 9225 K for the disentangled spectra of the components of V881 Sco covering the range from 4430 to 4490 Å. Based on these results, we adopted a temperature of 9100 ± 500 K, corresponding to a spectral type of A1/2V in the calibrated data of Pecaut & Mamajek (2013). This result fits the spectral type determined by Houk (1982) for V881 Sco, except for the luminosity class (see Section 1.3). Thus, the effective temperature of the primary component, which has a larger mass according to the RV curve solution, was taken as 9100 K and this value was fixed in the WD + MC iterations. We selected 8000–10 000 K as the input range for the effective temperature of the second component.

The input range of the orbital inclination for both eccentric systems was set to 75° < i < 90°, considering the light curve solution of Zasche (2012). For the changes in T₀ and P, the input range for the phase shift was set to −0.01 < Δϕ < 0.01. Based on the Korel solutions described in Section 2.2, the mass ratios (q) of GV Nor and V881 Sco were fixed at 1.03 and 0.98, respectively, in the first step. Accordingly, the range of the non-dimensional surface potential parameter of the components (Ω₁, Ω₂) was set to 5.0–8.0 for both systems.

Based on the RV solutions and O − C analyses of these two eccentric systems, the input ranges for the eccentricity (e) and longitude of the periastron (ω) were set to 0.09–0.15 and 140° − 180° for GV Nor, and to 0.09–0.13 and 20° − 40° for V881 Sco. Since both systems have eccentric orbits and exhibit apsidal motion, the synchronicity parameters² of the components were required for a conclusive analysis. Using the well-known formula for periastron-synchronization and assuming the eccentricity as 0.11 from the O − C analysis (Section 3.2), a synchronicity parameter F = 1.255 was adopted for both components of GV Nor, whereas the assumption of an eccentricity of 0.107 from the O − C analysis led us to adopt a synchronicity parameter F = 1.250 for both components of V881 Sco.

For both systems, a range from 0.3 to 0.8 was set for the fractional luminosity of the primary component. The input range for the third light was set to 0.01–0.20. This parameter was introduced to accommodate the large TESS pixels, which usually cause light contamination of the main target by nearby field stars appearing in the aperture, as well as potential third bodies in both the systems (see Section 3.2 for a discussion). A quadratic limb-darkening law was assumed; the limb-darkening coefficients (x_i and y_i) of the components were taken from Claret (2017). The limb-darkening law requires an effective wavelength, which we determined from the assigned surface temperatures T_eff and filter transmission data. The effective wavelength was close to 800 nm for the TESS photometry.

Our final model values for the TESS light curve and Korel RV curves of GV Nor and V881 Sco are listed in Table 5. The listed errors in the parameters correspond to a 90|${{\ \rm per\ cent}}$| confidence level. The first two rows of the table contain the light curve elements used for phasing both TESS LC data and RV data. These two parameters were taken from the O − C analysis in Section 3.2. The next five rows list the parameter values resulting from the second step of the WD+MC method, that is, RV curve solutions using ellipticity and proximity effects. In the second step of the WD + MC procedure for V881 Sco, the eccentricity converged to a value of 0.11; subsequently, the value was fixed at e = 0.105, the value generated in the first step. The remaining rows of the table list the photometric parameters resulting from the third/final step of the WD + MC method. The value of χ² was obtained from |$\sum _i({l}_{i,o}-{l}_{i,c})^{2}$|/|$\Delta {l}_{i}^{2}$| (Bevington 1969), where l_i, o and l_i, c are the observed and calculated light levels at a given phase, respectively, and |$\Delta {l}_{i}$| is an error estimate for the measured values of l_i, o. The reduced χ-squared value was obtained from |$\chi ^{2}_{red}$| = χ²/ν, where ν is the number of degrees of freedom of the data set.

The comparison of the Korel RV measurements and TESS photometric observations with the results generated by the WD + MC models are shown in Fig. 8 for GV Nor and in Fig. 9 for V881 Sco.

3.2 Apsidal motion

The respective apsidal motions of the two eccentric systems were investigated using the traditional O − C method, where O − C refers to the difference between observed and calculated times of eclipse / minimum.

In order to derive the apsidal motion of an eccentric orbit, numerous eclipse timings need to be accumulated over a substantial number of orbital cycles. We acquired these data from sources such as the ‘O − C gateway’³ (Paschke & Brat 2006) and ‘An Atlas of O − C Diagrams of Eclipsing Binary Stars’⁴ (Kreiner, Kim & Nha 2001).

Beyond these sources, the full data set we used in the present analysis of GV Nor included the very precise TESS light curves as well as the times of eclipses as derived from the ASAS-SN survey (Shappee et al. 2014; Kochanek et al. 2017), INTEGRAL-OMC (Mas-Hesse et al. 2003), the KELT (Kilodegree Extremely Little Telescopes) survey (Oelkers et al. 2018), and the older Bamberg photographic plates available online⁵ (Groote et al. 2014). For all these data sources (except the TESS data), we used our own automatic fitting procedure (AFP) Zasche et al. (2014) method for deriving the times of eclipses. The method uses the phased light curve as a template and the entire data interval divided into shorter parts.

For V881 Sco, in addition to the published data covering the past two decades that had been used previously to study the period, we re-calculated the older ASAS-3 data (Pojmanski 2002), together with the KWS survey (in both V and I_c filters, see Maehara 2014), and the ASAS-SN survey (Shappee et al. 2014; Kochanek et al. 2017). Furthermore, we also obtained new ground-based data with the FRAM telescope located in Argentina (The Pierre Auger Collaboration 2021). We incorporated the TESS data (Ricker et al. 2015) with the aforementioned sources when deriving the times of eclipses. The non-TESS data were used only for the derivation of times of eclipses and not for the LC modelling, because of their much poorer quality (both in terms of phase coverage and scatter) compared to the TESS data. Notwithstanding this approach, we were able to access additional data from the archive of old photographic plates from the Harvard observatory, DASCH (Digital Access to a Sky Century @ Harvard).⁶ When incorporating these data in conjunction with the AFP method, we obtained several eclipse timings of vital importance because of their locations in the time domain across the entire 20th century (see Fig. 11).

The aforementioned data were analyzed using standard methods (see, e.g. Giménez & García-Pelayo 1983; Giménez & Bastero 1995), yielding a set of parameters of the apsidal motion as listed in Table 6. A comparison of the observed data and the theoretical calculations is presented in Fig. 10 for GV Nor and in Fig. 11 for V881 Sco.

For GV Nor, a relatively slow period of approximately 224 yr was found for the apsidal motion; however, the uncertainty is very large, at approximately 60 yrs (see Fig. 10). Further observations would allow more precise estimates of the long-term evolution of the orbit. The on-going accumulation of new data would be expected to reveal the curvature of the O − C diagram in about 10 yr from now. Another question requiring additional (and precise) eclipse timing data is the validity of an apparent third-light component suggested by the data obtained since the year 2000.

For V881 Sco, the analysis implied an eccentricity e = 0.107 ± 0.017, in almost perfect agreement with the above-mentioned result obtained from the combined WD + MC fitting. The rate of change of ω was calculated as |$\dot{\omega }=0.014 \pm 0.003$| deg cycle⁻¹, producing a period of U = 175 ± 38 years for the apsidal motion. As illustrated in the bottom panel of Fig. 11, an additional eclipse time variation (ETV) is visible in the data after subtraction of the apsidal motion, with a suggested period longer than 20 yr.

The source of this variation (for example, an additional hidden component in the system or a simple artifact of the data) is an open question. Therefore, we also tested the significance of such an additional ETV variation in the O − C diagram using a Bayesian Information Criterion (Schwarz 1978) for both systems. However, for both stars this hypothesis failed to be proved and should therefore be rejected. The small ETV variation might be the result of scatter in the data, or some other undetermined effect. The nature of the third light derived from the TESS photometry (see the previous subsection) does not point unequivocally to a third body, primarily because of the large TESS pixel size. We leave this as an open question, subject to obtaining further data in the future to reveal the actual source.

4 RESULTS

4.1 Absolute parameters

Various physical parameters of eclipsing binary systems are obtained by combining photometric and spectroscopic observations and applying appropriate physical principles; the masses of the component stars (M₁, M₂) can then be derived by applying Kepler’s third law to the results obtained from the observational data. The fractional radius is given as r = R/a, where a is the semimajor axis of the binary system’s relative orbit and is determined from the solution of the RV curves (see Table 5). Thus, the fractional radii of the components (r₁, r₂), obtained from the photometric data, were used to derive the absolute radii (R₁, R₂). Subsequently, surface gravities (g₁, g₂) were derived directly. Using the effective temperatures listed in Table 7, the bolometric magnitudes (M_bol) and luminosities (L) of the component stars were computed. We used the nominal solar values adopted by IAU 2015 Resolutions B2 and B3 in our calculations. The absolute visual magnitudes M_V were derived from the bolometric correction formula, BC  = M_bol - M_V. Bolometric corrections for the components were taken from the online version of the color tabulation of Pecaut & Mamajek (2013), according to their effective temperatures.

For each system, its distance was calculated from the well-known adaptation of Pogson’s law: M_V = m_V  + 5 - 5log(d) - A_V. The interstellar absorption and intrinsic color index were computed using the following method: First, the total absorption towards the system in the galactic disk in the V band, A_∞(V), was taken from Schlafly & Finkbeiner (2011), using the NASA Extragalactic Database.⁷ Second, the interstellar absorption corresponding to the distance to the system, A_d(V), was derived from the formula given by Bahcall & Soneira (1980) (in their equation 8), using the system’s Gaia-DR3 parallax (Gaia Collaboration 2022). Finally, the colour excess for the system at the distance d was estimated as E_d(B−V)  = A_d(V)/3.1.

The distances to the systems – including the correction for interstellar absorption – were computed as 1100 ± 120 pc and 520 ± 80 pc for GV Nor and V881 Sco, respectively, using the distance modulus. When the interstellar absorption was ignored, the distances were derived as 1738 and 685 pc for GV Nor and V881 Sco. Considering the distances provided by Gaia-DR3 (Gaia Collaboration 2022), 1500 and 405 pc for GV Nor and V881 Sco, the interstellar absorption appears to be overestimated. The photometric parallax formulation of Popper (1998) produces distances of 1309 ± 80 pc for GV Nor (which is consistent with the distance of 1353 pc given by Gaia-DR2 [Gaia Collaboration 2018)] and 408 ± 40 pc for V881 Sco [which matches the distance of Gaia-DR3 (Gaia Collaboration 2022)].

The absolute parameters of GV Nor and V881 Sco as determined by our analysis, together with their errors, are listed in Table 7.

4.2 Evolutionary status

Granada evolution models (e.g. Claret 2006) were used to estimate the evolutionary status of the components of the two eccentric binaries. Specifically, the evolutionary tracks for the components were calculated by interpolating between models taken from Claret (2006) while applying the measured parameters listed in Table 7. Our analysis proceeded as follows:

• The metallicity Z of each system was estimated using the log (T_eff) – log (g) diagram. A range of metallicities with 0.007 ≤ Z ≤ 0.020 for main-sequence stars was considered. The optimal value of Z = 0.010 ± 0.001 determined for GV Nor corresponds to the tracks displayed in the upper panel of Fig. 12, whereas the optimal value of Z = 0.007 ± 0.001 determined for V881 Sco corresponds to the tracks displayed in the upper panel of Fig. 13. The solid black line in each figure represents the evolutionary track associated with the measured mass of the primary component and the dashed black line represents the track associated with the measured mass of the secondary component. In both diagrams, the primary components are located on their respective evolutionary tracks, whereas the locations of the secondary components relative to their tracks demonstrate that the observationally-determined values of the surface gravity are smaller than the theoretical values.

• The radial evolution associated with the measured masses of the components of GV Nor and V881 Sco for Z = 0.010 and Z = 0.007, respectively, were plotted in log (age) – radius diagrams, as shown in the centre panels of Figs 12 and 13. In Fig. 12, the measured radius of the primary component is attained at log (age) = 8.51, whereas that of the secondary component is attained at log (age)  = 8.56. The age of GV Nor was derived from the average of these values: 340 ± 35 Myr. As shown in Fig. 13, the measured masses and radii of both components of V881 Sco correspond to an age of 760 ± 15 Myr for the system.

• Finally, positions of GV Nor and V881 Sco in the H-R diagram are depicted in the bottom panels of Figs 12 and 13, for Z = 0.010 and Z = 0.007, respectively. Although the location of the primary component of GV Nor is consistent with its mass, temperature, and luminosity, the secondary component appears significantly more luminous than predicted for its measured mass. The primary component appears close to the middle of main sequence, whereas the secondary component is approaching the terminal-age main sequence (TAMS). Notwithstanding this discrepancy, the isochrone for an age of 340 Myr (Bressan et al. 2012), which was estimated from the log (age) – radius diagram, provides a good fit to the locations of the components of GV Nor in the H-R diagram. In the H-R diagram for V881 Sco, both components appear more luminous than predicted for their measured masses, within the error limits. Notwithstanding this discrepancy, the locations of the components of V881 Sco in the H-R diagram are in good agreement with the isochrone for an age of 760 Myr (Bressan et al. 2012), which was derived from the log (age) – radius diagram. Both components appear to be approaching the TAMS.

4.3 Internal structure constant

The second-order internal structure constant, k₂, is a parameter representing the internal density concentration of stars and is calculated by numerically integrating the Radau differential equation in theoretical stellar structure and evolution models (e.g. Kopal 1978). On the contrary, as can be seen from equation (2) below, the parameter k₂ is also dependent on the parameters |$c_{i}^{rot}$| and |$c_{i}^{tid}$|⁠, which play a role in the classical contribution (that is, the Newtonian contribution) to the apsidal motion of the common orbit in binary stars. Therefore, the analysis of observations of apsidal motion in eccentric eclipsing binary stars together with their spectral and photometric observations is the only way to directly obtain this internal structure constant. Thus, comparisons of k_{2, obs} derived from observational data with k_{2, theo} calculated from theoretical models provide tests of the validity of stellar models (e.g. Baroch et al. 2021; Claret et al. 2021).

The first step in calculating the mean observational internal structure constant |$\bar{k}_{2,obs}$| is to obtain the observed rate of apsidal motion |$\dot{\omega }_{obs}$| from the O − C analysis of the binary system. The relativistic term in the rate of apsidal motion |$\dot{\omega }_{rel}$| is calculated using equation (5). The classical term |$\dot{\omega }_{cl}$| is then derived from equation (1). The parameters c₁ and c₂ are calculated from equations (3) and (4), and |$\bar{k}_{2,obs}$| is obtained using equation (8). The mean theoretical internal structure constant, |$\bar{k}_{2,theo}$|⁠, is derived from equation (7). Here, the internal structure constants k₂₁ and k₂₂ of the component stars are taken from the theoretical evolution models. The observationally-determined and theoretically-derived values of the internal structure constants are compared as a test of the validity of stellar evolution models.

We predicted the theoretical internal structure constants |$\bar{k}_{2,theo}$| from the Granada evolution models (Claret & Giménez 1992; Claret 2006). The Granada models produce theoretical values of log k₂ for a grid of masses, surface gravities, luminosities, effective temperatures, ages, and metallicities. For a given binary system, values of k_2i for each component star were calculated by interpolating between the model values and the observed physical properties listed in Table 7.

Finally, parameters related to apsidal motion and observational and theoretical internal structure constants calculated for GV Nor and V881 Sco, following the path outlined in the three paragraphs above, are listed in Table 8. In this table, the metallicity and evolutionary ages determined for these two eccentric binary stars in Section 4.2 are also shown.

Table 8 indicates that the contribution of the relativistic term to the observed rate of apsidal motion for both systems was quite small, 5–6|${{\ \rm per\ cent}}$|⁠. On the conttrary, the difference between the observational and theoretical internal structure constants (Δlog k₂ = log k_{2, obs} − log k_{2, theo}) was approximately twice as large as the error of the observational internal structure constant for GV Nor (i.e. the difference was equal to approximately 2σ). For V881 Sco, this difference (Δlog k₂ = −0.128) was within the error of the observational internal structure constant (σ = 0.123). These results indicated that the theoretical internal structure constants predicted from the Granada models for these two eccentric systems were slightly larger than the observational values, almost within the margin of error. These small discrepancies imply that the predictions of the Granada models are compatible with the observed properties of GV Nor and V881 Sco.

5 REMARKS AND DISCUSSION

In previous studies of the apsidal motion of the system V881 Sco, the periastron longitude was quoted as 161° (Zasche 2012) and 186° (Kim et al. 2018), respectively, at the epoch HJD ≃ 2452128.4, whereas our Korel results implied a periastron longitude of 348° at the epoch HJD ≃ 2452126.5. When the epoch T₀ and orbital period P of V881 Sco were corrected, the light curve minima for V881 Sco were interchanged, so that the times of Min II recorded in previous studies became the times of Min I in this study and the times of Min I recorded in previous studies became the times of Min II in this study (in essence, the secondary component identified in previous studies is the primary component in this study).

High-precision TESS data combined with Korel RV data enabled the calculation of the masses and radii of the components with a precision of |$3{{\ \rm per\ cent}}$|⁠. The relative error (= error / value) was considered as a measure of precision. The highly precise absolute parameters calculated for GV Nor and V881 Sco showed good agreement with the predictions of Granada and Padova evolution models (see Figs 12 and 13). In particular, the fact that the radius growth curves of the components of GV Nor and V881 Sco with respect to their masses produced identical ages (within the errors) in the log (age) - radius diagram indicates that the Granada models represent the evolutionary stages of the components successfully. The Padova isochrones estimated from the locations of the components in the H-R plane matched the ages derived from the Granada models’ growth in component radii, within the errors.

Consequently, the successful prediction of the evolutionary status of the two eccentric systems by evolutionary models raises the question why the Granada models estimate a larger theoretical internal structure constant than the observational value (although within the margin of error). The most likely answer to this question is based on O–C analysis, from which apsidal motion parameters are derived. Considering Fig. 10 for GV Nor and Fig. 11 for V881 Sco, the sparseness and scatter of observational data in the first half of the cycle produced a large margin of error, especially in the calculation of the period of apsidal motion U and the rate of apsidal motion |$\dot{\omega }_{obs}$| (see Table 6). The high relative error of |$20{{\ \rm per\ cent}}$| in these parameters is also reflected in the calculation of the observational internal structure constant. Therefore, in order to be able to compare and interpret the parameters k_{2, obs} and k_{2, theo} more sensitively, it will be necessary to observe the photometric times of minimum of these two eccentric binary stars continuously for at least a few decades.

Another problem is that the evolutionary ages of these two binary stars, which have remarkably eccentric orbits (e ∼ 0.11), are quite large with respect to circularization dynamics. This raises the question whether the orbit circularization mechanism is indeed operating in these two systems. To investigate this problem, the circularization and synchronization timescales were calculated using Claret & Cunha (1997)’s formulation of the Zahn (1977) theory. The tidal torque constant E₂ and fractional gyration radius R_gyr in the formulas were determined by interpolating, according to the masses and surface gravities of the component stars, from the Granada evolution models (Claret 2006). By using log (E₂₁) = −7.05, R_{gyr, 1} = 0.20, log (E₂₂) = −8.14, and R_{gyr, 2} = 0.18 for the primary and secondary components of GV Nor and log (E₂₁) = −7.86, R_{gyr, 1} = 0.19, log (E₂₂) = −6.96, and R_{gyr, 2} = 0.18 for the components of V881 Sco, and the absolute parameters listed in Table 7, the equations of Claret & Cunha (1997) produced the values listed in Table 9 for the synchronization and circularization timescales of the component stars.

The following further two problems arise when considering the evolutionary ages and the critical ages of circularization and synchronization of GV Nor and V881 Sco: (i) Why do the components appear to rotate almost three times more rapidly than their synchronous rotation rates (see Section 2.3) even though their synchronization ages (calculated to be of the order of 10⁷ and 10⁶ yr, respectively) have been exceeded? (ii) The circularization timescales of the order of 10⁹ years for the primary components of both systems appear to be consistent with the evolutionary ages, which are of the order of 10⁸ years. However, the findings that the circularization timescales calculated with the contribution of primary and secondary components are of the order of 10⁸ years in both systems, shorter than their evolutionary ages, and that the orbits of both systems are still quite eccentric, are not compatible with the predictions of Zahn’s theory of circularization.

A possible explanation for this problem is that the dynamic effect of a possible third component in both systems maintains the eccentricities of the orbits of GV Nor and V881 Sco (e.g. Mazeh 1990). Unfortunately, the spectrometric and photometric analyses conducted in this study did not indicate the presence of a third body in these two eccentric systems. In the literature, there is also no astrometric study of the presence of a third body in these systems. Therefore, it will be necessary to carefully monitor the O − C changes of both eccentric systems and to investigate the light-time effect caused by the third body by using increasingly precise eclipse time data, in order to address these questions.

ACKNOWLEDGEMENTS

This research was supported by TÜBİTAK (Scientific and Technological Research Council of Türkiye) under Grant No. 121F203. PH acknowledges support by project RVO 67985815. CAE thanks the South African National Research Foundation and the University of Johannesburg for financial support. All the authors thank the SAAO for observing time. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. The authors thank the anonymous referee for guidance, which led to an improved paper.

AP is also acknowledged for sending us his photometric observations. The authors would also like to thank the Pierre Auger Collaboration for the use of its facilities. The operation of the robotic telescope FRAM is supported by the grant of the Ministry of Education of the Czech Republic LM2018102. The data calibration and analysis related to the FRAM telescope is supported by the Ministry of Education of the Czech Republic MSMT-CR LTT18004, MSMT/EU funds CZ.02.1.01/0.0/0.0/16|$\_$|013/0001402, CZ.02.1.01/0.0/0.0/18|$\_$|046/0016010, and CZ.02.1.01/0.0-/0.0/18|$\_$|046/0016007.

It is a pleasure to express the appreciation of the high quality and ready availability, via the Mikulski Archive for Space Telescopes (MAST), of data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. This research partly made use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund and NSF AST-1412587. This research was also partly based on data from the OMC Archive at CAB (INTA-CSIC), pre-processed by ISDC. The authors thank the ASAS, ASAS-SN, APASS, Bamberg, DASCH, FRAM, KELT, KWS, OMC, and TESS teams for making all of the observations readily accessible on public platforms. This work made partial use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research partly made use of the SIMBAD database, operated at CDS, Strasbourg, France, and of NASA’s Astrophysics Data System Bibliographic Services.

DATA AVAILABILITY

The majority of data included in this article are available as listed in the paper or from the online supplementary material it refers to. All the times of minima used for the apsidal motion analysis in this work will be made available at CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via https://cdsarc.unistra.fr/viz-bin/cat/J/MNRAS.

Footnotes

REFERENCES

APPENDIX A: TABLE OF RADIAL VELOCITY MEASUREMENTS