Issue 
A&A
Volume 649, May 2021



Article Number  A64  
Number of page(s)  14  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202040004  
Published online  12 May 2021 
Analysis of apsidal motion in eclipsing binaries using TESS data
I. A test of gravitational theories^{⋆}
^{1}
Institut de Ciències de l’Espai (ICE, CSIC), Campus UAB, c/ Can Magrans s/n, 08193 Bellaterra, Barcelona, Spain
email: baroch@ice.cat
^{2}
Institut d’Estudis Espacials de Catalunya (IEEC), c/ Gran Capità 24, 08034 Barcelona, Spain
^{3}
Centro de Astrobiología (CSICINTA), 28692 Villanueva de la Cañada, Madrid, Spain
^{4}
International Space Science Institute (ISSI), Hallerstrasse 6, 3012 Bern, Switzerland
^{5}
Instituto de Astrofísica de Andalucía (IAACSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain
Received:
27
November
2020
Accepted:
2
March
2021
Context. The change in the argument of periastron of eclipsing binaries, that is, the apsidal motion caused by classical and relativistic effects, can be measured from variations in the difference between the time of minimum light of the primary and secondary eclipses. Poor apsidal motion rate determinations and large uncertainties in the classical term have hampered previous attempts to determine the general relativistic term with sufficient precision to test general relativity predictions.
Aims. As a product of the TESS mission, thousands of highprecision light curves from eclipsing binaries are now available. Using a selection of suitable wellstudied eccentric eclipsing binary systems, we aim to determine their apsidal motion rates and place constraints on key gravitational parameters.
Methods. We compute the time of minimum light from the TESS light curves of 15 eclipsing binaries with precise absolute parameters and with an expected general relativistic contribution to the total apsidal motion rate of greater than 60%. We use the changing primary and secondary eclipse timing differences over time to compute the apsidal motion rate, when possible, or the difference between the linear periods as computed from primary and secondary eclipses. For a greater time baseline we carefully combine the highprecision TESS timings with archival reliable timings.
Results. We determine the apsidal motion rate of 9 eclipsing binaries, 5 of which are reported for the first time. From these, we are able to measure the general relativistic apsidal motion rate of 6 systems with sufficient precision to test general relativity for the first time using this method. This test explores a regime of gravitational forces and potentials that had not been probed before. We find perfect agreement with theoretical predictions, and we are able to set stringent constraints on two parameters of the parametrised postNewtonian formalism.
Key words: binaries: eclipsing / gravitation / relativistic processes / techniques: photometric
Full Table 3 is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/649/A64
© ESO 2021
1. Introduction
Eclipsing binaries have demonstrated to be a basic source of fundamental information about stellar properties such as masses and radii (Andersen 1991; Torres et al. 2010). The comparison of the observed values with theoretical predictions has been used profusely to perform critical tests of stellar structure and evolution models (Pols et al. 1997; Ribas et al. 2000; Torres & Ribas 2002; Lastennet & VallsGabaud 2002; Feiden & Chaboyer 2012; Higl & Weiss 2017; Tkachenko et al. 2020). But eccentric eclipsing binaries offer further opportunities to gain indirect insight into the internal structure of stars through the measurement of the precession rate of the line of apsides of the orbit, that is, the apsidal motion rate (). Such precession motion is found to arise from two different contributions: a general relativistic () term arising from general relativity (GR), and a classical or Newtonian () term. These two contributions are additive, so that (Shakura 1985). The general relativistic term of the apsidal motion rate, when only considering quadratic corrections, can be calculated with the equation given by LeviCivita (1937), in the form presented by Giménez (1985) in degrees per orbital cycle as:
where e is the orbital eccentricity, M_{1} and M_{2} are the component masses in solar units, and P_{a} is the anomalistic period in days, which measures the time between two consecutive periastron passages, and is related to the sidereal period, P_{s}, through:
where is in degrees per orbital cycle. The classical contribution to the apsidal motion is produced by perturbations in the gravitational potential due to the lack of spherical symmetry in the shape of the components, which are distorted due to rotational flattening and tidal oblateness, the socalled quadrupole effect, and for the most part depends on the degree of mass concentration towards the centre (Shakura 1985). The term , when only considering the contributions arising from the secondorder harmonic distortions of the potential, can be described by the expressions given by Sterne (1939) and Kopal (1959):
where the index i refers to the component and k_{2, i} are the secondorder internal structure constants characterising the internal mass distribution of the stars, which can be derived from theoretical models by numerically integrating the Radau differential equation (Kopal 1978; Hejlesen 1987; Schmitt et al. 2016). The parameters and are the rotational and tidal contributions to the apsidal motion, respectively, which, assuming that the stellar and orbital rotation axes are aligned, can be expressed (Kopal 1978; Shakura 1985) as:
where Ω_{m} = 2π/P_{s} is the mean angular velocity of the orbital motion, Ω_{i} = v_{rot, i}/R_{i} is the angular velocity of the rotation of each component, and r_{i} = R_{i}/a are the relative component radii. In Eq. (3) we do not use terms with k_{3} and k_{4} because they usually produce negligible contributions given the uncertainties of the lower order terms (Claret & Giménez 1993; Rosu et al. 2020). As all the terms appearing in Eq. (1) can be obtained from observations, the GR apsidal motion rate is most often calculated analytically and subsequently subtracted from the measured rate to compare with stellar model predictions and provide constraints on interior structure (e.g., Claret & Giménez 1993; Khaliullin & Khaliullina 2007; Rauw et al. 2016; Zasche & Wolf 2019).
Direct GR tests using measured apsidal motion rates have been hampered by the typical large relative uncertainty due to the error propagation from the dominant classical term. Eclipsing binaries suitable to test GR are those with relatively long orbital periods, where the general relativistic term makes up most of the apsidal motion contribution, that are still sufficiently short to produce an apsidal motion rate larger than the expected detection threshold (see Eqs. (16) and (17) in Giménez 1985, and their Table 1 to see the dependence of the period limits with the eccentricity). The most famous case is that of DI Her, whose detailed initial analyses showed disagreement with GR predictions (e.g., Guinan & Maloney 1985). It was later shown (Albrecht et al. 2009) that this mismatch was actually caused by a strong spin axis misalignment of the component stars, which leads to a different classical term, with a negative contribution of the rotational term to the total apsidal motion rate. There have been some attempts to further test GR using larger samples of eclipsing binaries but these failed to reach conclusive results (e.g., De Laurentis et al. 2012). This was caused by the large relative uncertainties resulting from the subtraction of the dominant classical contributions to the total apsidal motion term. The systems selected had a small fractional contribution of the GR term. With a carefully selected sample and the use of longer time baselines and more precise eclipse timings it should be possible to perform more accurate analyses and reach conclusive tests of GR. Precise determinations and even simply the detection of apsidal motion based on eclipse timings requires in most cases a dedicated longterm monitoring, generally spanning several decades. This is particularly difficult for systems whose orbital periods are long or have very slow apsidal motion rates, which happen to be the systems with the largest relativistic contributions due to the decrease of the classical terms with the orbital period as shown by Giménez (1985).
Data from the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2015) aimed at detecting exoplanets through photometric transits provide the opportunity to obtain denselycovered light curves of eclipsing binary systems at very high precision. Stellar eclipse events are generally much deeper than exoplanet transits and therefore TESS data are of excellent quality for such studies. Largely uninterrupted monitoring of the light curves is possible from space without the disturbing effect of the daynight cycle, therefore minimising the existence of phase gaps. In addition, most eclipsing binary systems were observed by TESS at the high cadence rate of 2 minutes. Thanks to all these circumstances, accurate eclipse timings can be obtained for numerous events. Furthermore, already during its initial 2 years of mission, TESS has covered a large fraction of the sky for a duration of at least 27 days, and therefore the chances of having observed eclipse events of most wellstudied eclipsing binary systems are very high. Of course, the time baseline of the typical TESS observations is rather modest, but the highprecision measurements can be combined with past timings to enlarge the covered timespan significantly.
In this paper we perform apsidal motion determinations in eccentric eclipsing binaries with relatively long orbital periods and for which general properties such as mass, radius, and orbital period are accurately determined. We restrict our analysis to those systems with a literaturepredicted relativistic contribution to the total apsidal motion rate of at least 60%. In some cases, we are able to detect apsidal motion for the first time. Particular care is put in the estimation of the classical terms to determine the GR term of the apsidal motion rate to the best possible accuracy and with realistic uncertainties, meaning that a comparison with theory can be performed. We analyse the individual times of eclipse for 15 eclipsing binaries, from which we aim to determine a precise value of the GR apsidal motion term for each system. We then compare these terms with predictions from different gravitation theories. In Sect. 2 we present the analysed sample and its general properties. In Sect. 3 we explain the methodology used to determine the times of eclipse, while in Sect. 4 we describe the methodology followed to determine the apsidal motion for each of the systems using TESS and archival timings. In Sect. 5 we present the resulting apsidal motion rates, the theoretical prediction of their classical contributions, and compare the differences with GR predictions. Finally, we use our measurements to test different gravitational theories in Sect. 6, and provide our concluding remarks in Sect. 7.
2. The eclipsing binary sample
For any useful interpretation of the observed apsidal motion in eccentric eclipsing binaries, it is essential to have good knowledge of the general properties of the component stars, mainly their masses and radii. Relevant equations such as Eqs. (4) and (5) have terms with high powers (up to the fifth) of the relative radii for example, thus enhancing their potential uncertainties. For this reason, we have limited our dynamical study to cases with wellstudied and precise fundamental properties. The basic source is the compilation of wellstudied detached eclipsing binaries by Torres et al. (2010). This has been complemented with systems from the DEBCAT catalogue (Southworth 2015), which is permanently updated, and includes eclipsing binaries with mass and radius measurements with 2% accuracy or better, following the criterion in Andersen (1991).
From the compilations above and imposing the restriction of having TESS measurements for both primary and secondary eclipses, we selected the sample of eccentric eclipsing binaries in Table 1. Given our objective of testing the theoretically predicted GR contribution, as given by Eq. (1), with observed apsidal motion rates, we limited the sample to only systems with a literaturepredicted fractional contribution greater than 60%. We used the general properties, as referenced in Table 1, and an estimated classical contribution given by Eq. (3), with the assumption of rotational synchronisation at periastron and coaligned rotational axes. Table 1 lists the adopted general properties of the systems such as the sidereal period, eccentricity, component masses, and radii of the 15 selected eclipsing binaries, together with the predicted GR contribution to the apsidal motion (%). It should be noted that in all tables presented throughout this paper the values in parentheses indicate the uncertainties affecting the last digits.
Properties of 15 eclipsing binary systems with eccentric orbits, accurate absolute dimensions, and literaturepredicted relativistic apsidal motion relative contribution greater than 60%.
Three additional eclipsing binary systems not included in Table 1 deserve specific attention. DI Her is known to have misaligned rotational axes (Albrecht et al. 2009). The computation of the classical term of the apsidal motion therefore requires the use of the general form of Eq. (3) given by Shakura (1985) and an analysis of those angles constrained but not directly observed with the RossiterMcLaughlin effect (Claret et al. 2010), which complicates the interpretation. EP Cru has component stars with rotational velocities 5.8 times higher than synchronisation at periastron (Albrecht et al. 2013). When such rotational velocities are considered, the classical term becomes dominant and the system does not meet our limiting criterion on the fractional GR contribution to the total apsidal motion rate to be included in the present study. Finally, a detailed analysis of BF Dra revealed a trend in the primary and secondary eclipse residuals of the highprecision TESS timings. Such a trend may likely be produced by the presence of a third body orbiting the system. The unconstrained nature of the third companion makes the determined apsidal motion rate illsuited for comparison with theoretical predictions. Eclipsing binaries with a relative contribution of the relativistic term of below 60%, together with DI Her, EP Cru, and BF Dra, will be analysed in a subsequent paper focusing on the classical terms (tidal and rotational).
3. Determination of times of minimum light
We used TESS data from sectors 1–29 to compute precise timings of the primary and secondary eclipses for all of our targets. We gathered the 2min shortcadence simple aperture photometry (SAP) produced by the Science Process Operation Centre (SPOC, Jenkins et al. 2016) available at the Mikulski Archive for Space Telescopes^{1}. In the cases of V530 Ori and V501 Mon, for which the 2min cadence photometry was not available, we extracted the 30min cadence simple aperture photometry from the TESS full frame images (FFIs) using the public TESS aperture photometry tool ELEANOR^{2} (Feinstein et al. 2019). Table 2 lists the TESS sectors and cadence at which each system has been observed, the timespan between the first and last eclipse considered for each target, and their Gaia G magnitudes (Gaia Collaboration 2016, 2018). Possible systematic deviations present in TESS data may include activityinduced modulations or other geometric effects that could bias the measurement of the time of minimum. To mitigate these effects, we employed the python package george (ForemanMackey 2015) to model the outofeclipse photometry using a Gaussian process correlatednoise model with a squaredexponential covariance function (see e.g., Gibson et al. 2011; Aigrain et al. 2016, for more details). The resulting model was used to normalise the entire light curve, including the eclipses. The lengthscale hyperparameter was constrained to values above twice the duration of the eclipse to avoid adding spurious highfrequency noise inside the eclipse region.
TESS sectors and cadence at which the sample of studied eclipsing binary systems have been observed, the timespan between the first and last eclipse included in this work, and their Gaia G magnitudes.
The Kwee & van Woerden (1956) method (hereafter, KvW) was adopted to determine all the times of minimum light. For consistency, we used the same orbital phase interval for all primary and secondary eclipses, and we only considered eclipses for which ingress and egress is well sampled. The median uncertainty in the obtained times of minimum are 2.5 and 12.9 s for the 2min and 30min cadence photometry, respectively. The precision of the KvW method for an equally sampled ideal eclipse is inversely proportional to the square root of the number of points used. Given that typically the number of TESS photometric points in the 2min cadence eclipses is 15 times larger than in the 30min cadence eclipses, we expect an error ratio between the two cadences of about , as approximately observed. In Fig. 1 we show two examples of eclipses observed with TESS using 2min and 30min cadence observations and we mark the time of minimum computed using the KvW method.
Fig. 1. Examples of TESS data for an eclipse of V459 Cas (a) with 2min cadence, and V501 Mon (c) with 30min cadence. In panels (b) and (d) the black line shows the bisector corresponding to each eclipse. The computed time of minimum for each system, which is used as reference time in the plots, is marked as a blackdashed line in panels a and c. Its uncertainty is illustrated in panels b and d as a read shaded area. The resulting timing uncertainty is 3.3 s for V459 Cas and 22.0 s for V501 Mon. 
We computed the bisector of each eclipse to diagnose local asymmetries that could be caused by stellar activity for example. Bisector points were determined as the time average between two symmetric points of the ingress and egress branches. Cubic spline interpolation was used to define points of equal flux. From the calculated values, we defined two indicators. The bis_{odd} parameter is defined as the difference of the average bisector between 50–100% and 0–50% of the eclipse depth (i.e. ⟨bis_{50−100}⟩−⟨bis_{0−50}⟩), while the bis_{even} parameter is computed as the difference between the average of the two intervals at the top 75–100% and the bottom 0–25%, and the average of the middle 25–75% of the eclipse depth (i.e. ⟨bis_{75−100}, bis_{0−25}⟩−⟨bis_{25−75}⟩). These bisector parameters represent first and secondorder measures of distortions in the eclipse shapes. We used them as signposts for the potential presence of biased timing determinations, for example those due to issues with the data or effects on the stellar surface such as spots. As a result, we did not consider for further analysis those times of minimum with one of these indices deviating by more than 3σ from their mean. Table 3 provides all the measured times of minimum light, the type of eclipse, and the bisector indicators bis_{odd} and bis_{even}. Also, the right panel of Fig. 1 shows two examples of the bisectors of two eclipses and the uncertainty in the determination of their time of minimum.
In the eclipses analysed in this work the final measurement errors are 1–11 times lower than the dispersion of the individual points used to measure the bisector. Given that the number of photometric points comprised within the eclipses ranges from 8 to 440, that the bisector points are computed from two photometric points, and that the timing error for an ideal eclipse scales with N^{−1/2}, we expect errors that are 2–s15 times lower than the bisector dispersion (i.e. and ). When applied to our data, this methodology yields uncertainties that are larger than those expected for an ideal eclipse, which may be taken as an indication that the errors are not underestimated. Potential systematic deviations from an ideal eclipse that are not accounted for by the KvW method caused several works to consider the error estimates coming from this method as exceedingly optimistic (Breinhorst et al. 1973; Mikulášek et al. 2014). However, with the novel use of the bisectors as indicators for such deviations, we were able to detect the affected eclipses and discard them for further use, therefore ensuring that our measurement errors are a good representation of the true statistical precision.
4. Apsidal motion determinations
Several methods can be employed to estimate apsidal motion rates in eclipsing binaries, such as the changing position of the eclipses with respect to each other over time (e.g., Giménez et al. 1987; Wolf et al. 2006, 2010; Kim et al. 2018), the change in the shape of the radial velocity curve over time (e.g., Ferrero et al. 2013; Schmitt et al. 2016; Rauw et al. 2016), or the change in the shape of the light curve over time (Bakış et al. 2008; Harmanec et al. 2014; Torres et al. 2017). Here we use the equations of Giménez & Bastero (1995) defining the times of eclipse:
where is expressed in deg cycle^{−1}, j is 1 or 2 for the primary or secondary eclipse, respectively, T_{0} is the time of eclipse at epoch N = 0, P_{a} is the anomalistic period, and ω is the argument of periastron expressed in radians at time T, which depends on and the argument of periastron at time T_{0}, ω_{0}, as . Finally, the coefficients A_{i} depend on the eccentricity e and the inclination of the orbit, and their full expressions are given in Eqs. (16)–(21) of Giménez & Bastero (1995).
Equation (6) is complete up to 𝒪(e^{5}) (Giménez & Bastero 1995), and is therefore only valid for orbital eccentricities that do not reach extremely high values, that is, below 0.7. Fitting this equation to the actual individual measurements allows the determination of the apsidal motion rate, , as well as the orbital eccentricity, e, the orbital period, the reference time of eclipse, T_{0}, and the argument of periastron at that time, ω_{0}. It should be noted that there are significant differences in the calculated when only expressions considering low orders in e are used. The severity of the deviations is a function of the eccentricity but we use terms up to e^{5} for all cases, implying relative errors well below 1%.
Strong degeneracies between e and ω_{0} can limit the use of Eq. (6) to derive all parameters simultaneously unless a significant fraction of the apsidal motion period is covered with observational data. Attempts at fiveparameter simultaneous fits using a poorly covered apsidal period may lead to biased solutions because of resulting correlations between e and ω. Potential issues are further exacerbated by the possible presence of lighttime effect variations caused by an orbiting third body. To obtain nondegenerate results when the timespan of the observations is limited, the value of the eccentricity, which is directly related to the amplitude of the variations, should at least be adopted independently, for example from modelling the light and radial velocity curves. When the eccentricity is fixed, the instantaneous shape of the apsidal motion curve can be described by the tangent at ω_{0}, which can also be expressed in terms of the difference between the periods corresponding to primary and secondary eclipses, as described by Giménez & Bastero (1995), considering that the value of ω at the time of the eclipses does not differ from ω_{0}. This is just the timederivative of the difference between the times of secondary (T_{2}) and primary (T_{1}) eclipses:
with expressed in deg cycle^{−1}.
In our study, we use the values of T_{2} − T_{1} from pairs of close primary and secondary minima of highprecision data. This is done to ensure maximum robustness of the apsidal motion determinations and to minimise the possible effects over the measured orbital period that might be produced by systematic errors arising from different methodologies or astrophysical variability. We favour this method over the use of the difference between periods resulting from linear fits to individual primary and secondary eclipses. This is equivalent to the right term of Eq. (7). Using T_{2} − T_{1} from close eclipses has several advantages, namely that these are more likely to be obtained by the same observer, that both timings correspond to the same epoch of stellar activity, and that it avoids effects from lighttravel time caused by a potential third body. However, it should be noted that the presence of a third body may induce other perturbations such as the eccentric LidovKozai effect (Lidov 1962; Kozai 1962; Naoz 2016) that could alter the observed apsidal motion rates. A full characterisation of the system is needed to take into account the effects of a potential third body. In the absence of such a detailed analysis, a comparison with the theoretically predicted apsidal motion rate could be affected.
The downside our adopted approach is that the number of T_{2} − T_{1} measurements (TESS plus archival) for some of the studied systems presented below is rather low. Depending on the time baseline and the impact of effects such as stellar activity, this may lead to some unaccountedfor source of additional error. However, our entire procedure has been devised to minimise such effects (strict selection of minimum timings, bisector correction criteria, etc.) and therefore we do not expect systematic deviations to significantly affect our measurements. Also, most of the eclipsing binaries that we analyse here do not show signs of any stellar activity, either because they lack a convective envelope or because of their low rotational velocities.
The eclipsing binaries in our sample have been typically monitored between TESS and archival data for less than 1% of their apsidal motion period, of the order of thousands of years. Therefore, the observed T_{2} − T_{1} values should accurately define the tangent to the curve described by Eq. (7), and relying on a previously determined orbital eccentricity, Equation (7) is fully applicable. The approximation of the derivative is of course sensitive to the argument of periastron and its precision decreases for values of ω_{0} near 0 or 180 degrees, when the variation of T_{2} − T_{1} over time and the corresponding period differences become close to zero. We assume the eccentricity derived from the best light and radial velocity curves, with the corresponding argument of periastron, in order to obtain the apsidal motion rate. We had to restrict the analysis to systems with highprecision data, as already discussed, and we focus our analysis on the systems listed in Table 1, which are discussed individually below.
All obtained secondary minus primary eclipse timings, T_{2} − T_{1}, for each of the available primary eclipses are given in Table 4 with the corresponding uncertainty computed from the propagation of errors of individual timings, given in parentheses. We also list the orbital cycle (N) computed with respect to the reference time T_{0} given in the first column. For those primary eclipses without any secondary counterpart on the same orbital cycle, we computed the secondary minus primary eclipse timings as T_{2} − T_{1} − dN × P_{s}, where P_{s} is the sidereal period listed in Table 1 and dN is the difference in orbital cycles between the two minima. We list dN in the last column in Table 4.
T_{2} − T_{1} computed from TESS light curves.
Below we discuss the studied systems in detail and describe the methodology applied to determine their apsidal motion rates in those nine cases where such a measurement was possible. For 6 out of the 15 systems in Table 1, namely, RW Lac, V530 Ori, HP Dra, TZ Men, LV Her, and RR Lyn, we could not find a measurable apsidal motion rate in spite of having highprecision TESS data available, but we have nevertheless included their time of minima and T_{2} − T_{1} values in Tables 3 and 4, which could be valuable for future determinations. For the systems in which we could determine an apsidal motion rate, additional values needed for their analysis, such as the relative radii or the projected rotational velocities, are listed in Table 5, together with the values of ω_{0} and the reference time T_{0}.
Other properties used to compute the apsidal motion of the systems analysed.
4.1. KX Cnc
The general properties of this highly eccentric system (e = 0.4666 ± 0.0003), such as mass, radius, eccentricity, and period, shown in Table 1, are the result of the combined analysis of the light and radial velocity curves carried out by Sowell et al. (2012), who reported no apsidal motion. The TESS data listed in Table 4 yield T_{2} − T_{1} = 20.077329 ± 0.000036 days. This measurement can be compared with T_{2} − T_{1} values from the literature to search for the presence of apsidal motion. Some published minima by Davies (2007) have insufficient precision and we therefore considered the light curve obtained by Sowell et al. (2012). The authors give a phase for the secondary eclipse of 0.64325, which is not measured directly but results from the joint analysis of the light curve and the radial velocity curve. We therefore decided to retrieve the original photometric data, b and yband in the Stromgren system, and calculated the position of the eclipses using the same method as for the TESS measurements, yielding a value of T_{2} − T_{1} = 20.08076 ± 0.00015 days, which is equivalent to a phase of the secondary of 0.64321. Comparing the value of T_{2} − T_{1} with the TESS results, we observe a decrease of −0.00343 ± 0.00015 days after 177 orbital cycles, indicating a measurable apsidal motion. With only two data points available, an analytical derivation of the statistical uncertainty in the slope could be imprecise and biased. Therefore, we performed 10^{6} Monte Carlo simulations of T_{2} − T_{1}, following a Gaussian distribution with the corresponding uncertainty, and computed the standard deviation of all the solutions. From the resulting slope, and combining it with the eccentricity given by Sowell et al. (2012) from their joint solution, an apsidal motion rate of deg cycle^{−1} is derived, including the uncertainty in the adopted orbital parameters.
4.2. AL Dor
The general properties given in Table 1 are derived from the radial velocity curve of Graczyk et al. (2019), combined with the astrometric solution by Gallenne et al. (2019). No indication of apsidal motion was reported. The adopted parameters, including the sidereal period, are not based on a lightcurve solution. Fortunately, our TESS eclipse timings, shown in Table 3, yield a good number of accurate minima values indicating a slow apsidal motion without any need for archival measurements, as illustrated in Fig. 2. We tried to enlarge the time baseline to improve the apsidal motion determination, but we could not find earlier data of sufficient precision. The linear weighted fit to the T_{2} − T_{1} values listed in Table 4 corresponds to a slope of −(5.0 ± 0.5) × 10^{−6} days cycle^{−1}. The uncertainty of the sidereal period is large compared to the errors of periods based on a lightcurve solution. For this reason we only used values of T_{2} − T_{1} within the same orbital cycle to eliminate the effect of the uncertainty in the sidereal period. Complementarily, we considered all the individual timings to fit the periods corresponding to primary and secondary eclipses, obtaining a difference between them of −(5.08 ± 0.19) × 10^{−6} days cycle^{−1}, in agreement with our value based on measurements within the same orbital cycle. The eccentricity given by Gallenne et al. (2019), e = 0.1952 ± 0.0002, produces a final value of deg cycle^{−1}.
Fig. 2. Ephemeris curve for AL Dor as a function of the orbital cycle. The dotted and dashed lines represent fits to the primary (blue circles) and secondary (red triangles) eclipses, respectively. For a better visualisation, an arbitrary shift to the primary and secondary minima has been applied. 
4.3. V541 Cyg
Apsidal motion in V541 Cyg was measured by Khaliullin (1985) and revised by Volkov & Khaliullin (1999). Wolf et al. (2010) used the individual times of minimum light to obtain = 0.00032 ± 0.00006 deg cycle^{−1}, leaving the eccentricity as a free parameter, while Kim et al. (2018) obtained deg cycle^{−1}, fixing the orbital eccentricity to the value of 0.479 given by Lacy (1998). Torres et al. (2017) obtained a new radial velocity curve and reanalysed the Vband light curve of Khaliullin (1985). For the apsidal motion determination, Torres et al. (2017) carried out a global analysis of all the data, photometric and spectroscopic, with variable argument of periastron. In this way, the authors determined the optimal value of the eccentricity to be e = 0.4684 ± 0.0014 and an apsidal motion rate of deg cycle^{−1}.
For our analysis, we included in Table 6 the T_{2} − T_{1} values derived from Table 4 of Torres et al. (2017) which we complement with the values computed from the TESS data in Table 4. This was done using precise primary and secondary eclipses and the closest possible counterpart when calculating the timing differences. As mentioned by Torres et al. (2017), some values had to be excluded because of their unusually large residuals. In Fig. 3, the variation of T_{2} − T_{1} over time is shown and clearly indicates a welldefined increase of (3.12 ± 0.03) × 10^{−5} days cycle^{−1}. For the value corresponding to orbital cycle 19, the authors did not report any error bar and we assumed the same uncertainty as the largest one in Table 6, namely 0.003 days. As a simple consistency test, we checked that the observed slope yields T_{2} − T_{1} = 6.9872 ± 0.0003 days at the time of the photographic light curve of Karpowicz (1961), while the actual measurement given by Torres et al. (2017) in their Table 4, after a careful analysis of the original data, is 6.988 ± 0.004 days. Both values are in very good agreement in spite of the extrapolation by about 550 orbital cycles. Furthermore, the predicted position of the secondary eclipse in phase at epoch 248 is 0.45798 ± 0.00002, in excellent agreement with the lightcurve solution by Torres et al. (2017) in their Table 5. Adopting e = 0.4684 ± 0.0014 (Torres et al. 2017), we obtain an apsidal motion rate of deg cycle^{−1}, more precise but in good agreement with the value of deg cycle^{−1} derived by Torres et al. (2017) using a completely different method, based on a simultaneous fit to photometric and spectroscopic data.
Fig. 3. T_{2} − T_{1} as a function of the orbital cycle for V541 Cyg. 
T_{2} − T_{1} values used to compute the apsidal motion rate of V541 Cyg.
4.4. V459 Cas
Lacy et al. (2004) determined the general properties of V459 Cas and estimated the apsidal motion rate using all available minima at that time, obtaining a rather uncertain value of deg cycle^{−1}. This could not be improved by Dariush et al. (2006) because of the narrow timespan of their reliable data. Torres et al. (2010) nevertheless quote an apsidal motion rate of 0.00057 ± 0.00006 deg cycle^{−1} as provided preliminarily to the authors by M. Wolf. This value was not confirmed later by Wolf et al. (2010), who reported 0.00071 ± 0.00008 deg cycle^{−1}, using again all available individual timings at that time.
The low eccentricity of 0.0244 ± 0.0004 (Lacy et al. 2004) of the system makes it difficult to determine the apsidal motion rate using the observed T_{2} − T_{1} values. We searched in the list provided by Wolf et al. (2010), in their Table A.1, and in Table 1 of Lacy et al. (2004) for the closest photoelectric pair of primary and secondary timings, which we list in Table 7. The best linear fit to these values, together with the TESS T_{2} − T_{1} measurements, is shown in Fig. 4, and yields a slope of (1.30 ± 0.19) × 10^{−6} deg cycle^{−1}, although with a large dispersion of the older timings due to the small variation in T_{2} − T_{1} in the considered timespan. Using the orbital eccentricity of Lacy et al. (2004), we obtain an apsidal motion rate of deg cycle^{−1}, which agrees with the value given by Wolf et al. (2010), although with a slightly larger uncertainty.
Fig. 4. T_{2} − T_{1} as a function of the orbital cycle for V459 Cas. 
T_{2} − T_{1} values used to compute the apsidal motion rate of V459 Cas.
4.5. V501 Her
The absolute parameters of V501 Her are given in Table 1 as derived by Sandberg Lacy & Fekel (2014) from their own lightcurve and radial velocity measurements, from which no apsidal motion could be reported. The large relative radii indicate the evolved nature of the components and the difficulty in obtaining accurate times of eclipse due to their long duration. The TESS measurements given in Table 4 present an internal dispersion that is larger than the estimated errors of the individual measurements. The weighted mean value is T_{2} − T_{1} = 4.1208 ± 0.0004 days, equivalent to a phase of the secondary eclipse of 0.47929 ± 0.00005 using the orbital period given by Sandberg Lacy & Fekel (2014). These authors give a secondary eclipse phase of 0.4791 ± 0.0002, implying that no significant variation is observed. From the most complete and precise light curve in Sandberg Lacy & Fekel (2014) and using the elements in their Table 4, we estimated the phase of the secondary eclipse at 0.47915 ± 0.00005, compatible with but more precise than the one given by their times of eclipse, and equivalent with a T_{2} − T_{1} of 4.1196 ± 0.0004 d. Comparing this value with that from the TESS observations, we measure a change in T_{2} − T_{1} of 0.0012 ± 0.0006 days over 354 orbital cycles, yielding a slope of (3.4 ± 1.7) × 10^{−6} deg cycle^{−1}. This shows the presence of apsidal motion in V501 Her but with a poorly determined rate. The corresponding apsidal motion rate is deg cycle^{−1}.
4.6. KW Hya
The only accurate masses and radii for this system date back to Andersen & Vaz (1984), resulting from the analysis of both light and radial velocity curves. Such results were confirmed more recently by Gallenne et al. (2019) with astrometric observations, and we adopted those values as listed in Table 1. Apsidal motion was not reported for KW Hya, despite the large timespan between the two orbital studies, because of the uncertainties involved in the determination of the argument of periastron.
The weighted mean of the TESS measurements in Table 4 provides an accurate value of the T_{2} − T_{1} difference of 3.54987 ± 0.00004 days. This can be compared with the corresponding timing difference measured by Andersen & Vaz (1984), 1740 orbital cycles before, of 3.5454 ± 0.0007 days. No other value of sufficient precision could be found in the literature. These two available measurements yield a slope of 2.6 ± 0.4 × 10^{−6} days cycle^{−1}. Adopting the eccentricity given by Gallenne et al. (2019), e = 0.094 ± 0.004, gives an apsidal motion rate of days cycle^{−1}.
4.7. V501 Mon
A detailed analysis of the light and radial velocity curves of V501 Mon was carried out by Torres et al. (2015) and the corresponding absolute parameters are given in Table 1. In their study, an analysis of the available times of eclipse at that time was carried out together with the radial velocities in order to determine the orbital elements and the possible variation in the argument of periastron simultaneously. The authors obtained an estimate of the apsidal motion rate of 0.00045 ± 0.00024 deg cycle^{−1}, which we now try to improve by combining the TESS measurements given in Table 4 with the data used by Torres et al. (2015). The weighted mean of the TESS T_{2} − T_{1} values gives a separation between primary and secondary eclipses of 3.14676 ± 0.00020, which does not show significant variation with respect to the solution by Torres et al. (2015).
We then fitted all the individual timings corrected to BJD in order to identify a possible difference between the linear periods of the primary and secondary eclipses. Using all the times of eclipse available in Table 1 of Torres et al. (2015; adopting also the same scaling for the photoelectric timings) and the TESS individual eclipses, we computed the periods resulting from primary and secondary eclipses, and obtained a difference of ΔP = 3.85 ± 0.89 × 10^{−6} days cycle^{−1}, which corresponds to the lefthand side term in Eq. (7) (Giménez & Bastero 1995). With the value of the eccentricity e = 0.1339 ± 0.0006, this ΔP value yields an apsidal motion rate of deg cycle^{−1}, which is in excellent agreement with the value published by Torres et al. (2015) but with a reduced uncertainty. Figure 5 displays the linear best fit to primary and secondary individual timings. The data points on the left, showing a large dispersion, correspond to the photographic measurements in Table 1 of Torres et al. (2015), for which we assumed uncertainties of 0.025 d.
4.8. GG Ori
The general properties of GG Ori were determined by Torres et al. (2000), as listed in Table 1, together with an apsidal motion rate of deg cycle^{−1}. Torres et al. (2000) followed a method based on the combination of eclipse timings with radial velocities, the same used in their analysis of V501 Mon mentioned above. Later, Wolf et al. (2010) revised the analysis of the individual eclipse timings available at the time and obtained a more precise deg cycle^{−1}, with an orbital eccentricity of e = 0.220 ± 0.001.
In order to improve the apsidal motion rate determination combining early timings with TESS, we compiled T_{2} − T_{1} measurements from the best minima given by Wolf et al. (2010) acquired within less than ten orbital cycles, but only three pairs of primary and secondary eclipses met this criterion. We derived deg cycle^{−1}. We then performed linear fits to primary and secondary individual timings, respectively, obtaining the linear period for primary and secondary eclipses and obtained a difference of ΔP = −8.33 ± 0.16 × 10^{−6} days cycle^{−1} using only those eclipses with weight ten in Wolf et al. (2010). Furthermore, using the eccentricity obtained by Torres et al. (2000), e = 0.2218 ± 0.0022, we obtain an apsidal motion rate of 0.00061 ± 0.00003 deg cycle^{−1}, in good agreement with the previous determinations albeit much more precise. The ephemeris curve of the individual primary and secondary eclipse timings is shown in Fig. 6.
4.9. EY Cep
This binary system was studied by Lacy et al. (2006), who obtained the general properties given in Table 1. Given the short timespan of their photometric and spectroscopic observations, the authors were not able to detect any indication of apsidal motion. The new TESS data listed in Table 4 include two sectors and the T_{2} − T_{1} variation already indicates the presence of apsidal motion, though the modest timespan does not allow a precise rate determination. We further used the eclipses given by Lacy et al. (2006) to obtain the time differences in Table 8.
T_{2} − T_{1} values used to compute the apsidal motion rate of EY Cep.
The graphical representation of the variations is shown in Fig. 7 and the linear leastsquares fit yields a welldefined slope of (−1.97 ± 0.06) × 10^{−5} days cycle^{−1}. Using the eccentricity given by Lacy et al. (2006), e = 0.4429 ± 0.0014, this yields an apsidal motion rate of deg cycle^{−1}, which is very precise thanks to the accurate TESS data, the high orbital eccentricity, and a time coverage of more than 800 orbital cycles.
Fig. 7. T_{2} − T_{1} as a function of the orbital cycle for EY Cep. 
5. Comparison between theory and observations
The values of obtained from the apsidal motion analysis using data from TESS and archival times of minima are listed in Table 9. Of the nine systems studied for which we were able to perform apsidal motion rate determinations, five are reported for the first time. The precision of the determinations varies widely but they are well above the 3σ threshold, except for the case of V501 Her, which has an apsidal motion rate determination with a significance of only 2σ. Table 9 provides the values of and for the systems with detected apsidal motion, computed using Eqs. (1) and (3). Calculating the classical term requires employing two additional parameters for each component, namely the internal structure constant, k_{2}, and the rotational velocity.
log k_{2} values and apsidal motion rates, both observed and theoretically predicted, for the eclipsing binaries with apsidal motion measurement.
The values of k_{2} given in Table 9 were specially calculated by A. Claret using theoretical models based on the Modules for Experiments in Stellar Astrophysics package (MESA; Paxton et al. 2011, 2013, 2015) and following the methodology described in the series of papers by Claret & Torres (2017, 2018), 2019). A coarsegrid search was performed over evolutionary tracks calculated for the measured masses of each component, allowing the convective core overshooting parameter (f_{ov}) and the mixing length parameter (α_{MLT}) to vary freely, with a variable metallicity (Z) common to both components. The log k_{2} values were determined from the best match between the grid of evolutionary tracks and the observed masses, radii, and effective temperatures. The theoretical internal structure constants, k_{2}, for the models fitting the observed parameters, were integrated using the differential equations of Radau as given by Eqs. (1)–(3) in Claret & Giménez (2010).
The other parameter used to estimate is the rotational velocity. The values listed in Table 5 were taken from the same spectroscopic analyses in the literature from which we adopted the general properties in Table 1. Nevertheless, in the case of AL Dor, no rotational velocity was reported for the component stars and we used the predicted values under the assumption of pseudosynchronisation, as described by Hut (1981).
All the parameters needed to apply Eqs. (1) and (3) are listed in Tables 1 and 5. The errors of the input parameters were propagated to obtain the uncertainty in and . The uncertainties of the stellar rotation and k_{2} dominate the error budget in the classical term, while the uncertainties in the component masses are the main source of error in the GR term. In Fig. 8 we compare the observed apsidal motion rates with the total calculated theoretical value , as given in Table 9. We excluded V501 Her from Fig. 8 given the poor significance of the apsidal motion determination (below 2σ), although the observed value agrees with theory within the uncertainties.
Fig. 8. Comparison between total computed and total observed apsidal motion rates. The black dashed line indicates the 1:1 relation. 
The comparison of the observed apsidal motion rates with those calculated using Eqs. (1) and (3), and the parameters in Tables 1 and 5, show good agreement and no systematic deviations within their uncertainties. An equivalent approach to test GR effects was employed before by Torres et al. (2010), who also considered systems with accurate masses and radii and limited the comparison between observed and theoretical apsidal motion rates to those having a relative contribution of the GR term of at least 40%. Figure 8 therefore complements Fig. 11 in Torres et al. (2010), including additional systems with a relative contribution of the GR term above 60%. The values used by Giménez (1985, 2007) for the system’s general properties were of insufficient precision to allow for a meaningful and unbiased test of the GR term. Claret & Giménez (1993) searched for systematic deviations in the comparison between observed and theoretical values of log k_{2} to falsify the predictions of the Moffat (1986) theory of gravitation. Claret (1997) and Wolf et al. (2010) studied the complementary problem, where the GR term was theoretically estimated and subtracted from the observed value. The resulting classical term was then compared with model predictions of the k_{2} parameter.
6. A test of gravitational theories
For a more indepth analysis, the two apsidal motion components should be analysed separately. As explained above, we carefully selected the sample to ensure a maximum contribution of the GR term, so that, after subtraction of the (modelcalculated) quadrupole Newtonian contribution, a comparison with gravitational theories could be carried out. De Laurentis et al. (2012) attempted a test of gravitational theories using a similar method but their eclipsing binary sample was not appropriate, with over half of their sample having a GR term contribution below 10%. Their results were therefore inconclusive.
Thanks to the improved precision in the apsidal motion determinations that we have attained, we can perform a comparison of the predicted theoretical relativistic apsidal motion rate with the difference between the observed values and the computed classical terms, as listed in Table 9. For this procedure to succeed, the uncertainty in the classical contribution needs to be well below this latter difference and therefore this limits the comparison to systems with a small nonrelativistic contribution to the observed apsidal motion. In Table 9 we also include the measured GR rate calculated as , where is the observed apsidal motion from precise minima timings derived in this work. Given the stringent requirements on the system characterisation (to ensure good estimates of the internal structure parameters and classical apsidal motion contribution), we excluded the cases of V459 Cas and V501 Mon in spite of their good agreement shown in Fig. 8.
Figure 9 shows the comparison between predicted , according to Eq. (1), and . The measured values for all systems are compatible with the GR predictions within their errors. A more general form of Eq. (1), which also allows to test other possible gravitation theories in addition to GR is given by Eq. (66) of Will (2014). This uses the parametrised postNewtonian (PPN) formalism, with several parameters whose values depend on the gravitation theory chosen:
Fig. 9. Measured GR apsidal motion as a function of the theoretically computed GR apsidal motion. The black dashed line and grey shadow area correspond to the best fit to Eq. (8), assuming α_{i} = ζ_{2} = 0, and its 1σ uncertainty. The colour code is chosen such that redder colours correspond to systems with the smaller relative error, and therefore are dominating the fit. 
Here, is the value of the apsidal motion predicted by GR, which is given by Eq. (1), η ≡ M_{1}M_{2}/(M_{1} + M_{2})^{2} is the dimensionless reduced mass, and γ, β, ζ_{2}, and α_{i} are the PPN parameters. In any fully conservative theory of gravity, α_{i} = ζ_{2} ≡ 0, while for GR, γ = β ≡ 1 and α_{i} = ζ_{2} ≡ 0, recovering the expression in Eq. (1). The PPN formalism is appropriate for weak gravitational fields and slow motions, such as the conditions in stellar eclipsing binary systems (see Will 2014, for more details).
Therefore, a test of GR can be performed by checking that the measured and theoretical GR apsidal motion rates are compatible with values of α and β predicted by the theory, assuming a fully conservative model (i.e. α_{i} = ζ_{2} ≡ 0). The black dashed line shown in Fig. 9 corresponds to the best fit to the measured and predicted values of the GR apsidal motion, assuming only a varying slope. This slope, as deduced from Eq. (8), corresponds to A ≡ (2 + 2γ − β)/3. We derive a value of A = 1.002 ± 0.012, which is fully compatible with the value predicted by GR of 1. It is clear from Fig. 9 that the stronger constraints are set by the systems with the smaller relative uncertainties, that is, AL Dor and V541 Cyg. However, it should be noted that the fit is still compatible with GR when removing these two targets, obtaining A = 1.015 ± 0.036.
From our determination of A, the validity of alternative gravitational theories in the PPN formalism can be assessed. For example, constraints can be put on the coupling constant ω_{BD} of the BransDicke theory (Estabrook 1969), which takes A = (4 + 3ω_{BD})/(6 + 3ω_{BD}). Furthermore, our measurements can also be used to place bounds on the standard individual PPN parameters. Adopting the limit of (γ − 1) from the Shapiro timedelay measurements using the Cassini spacecraft, of 2.3 × 10^{−5} (Bertotti et al. 2003), our measurement yields a constraint on β given by (β − 1) = − 0.005 ± 0.035. Equivalently, adopting the limit of (β − 1) from the perihelion shift of Mercury, of 8 × 10^{−5} (Verma et al. 2014), we find a constraint to γ of (γ − 1) = 0.002 ± 0.017.
A test of nonconservative models can be carried out by measuring the value of B ≡ (2α_{1} − α_{2} + α_{3} + 2ζ_{2})/6 in Eq. (8). Assuming now γ = β = 1, we can determine the value of B by fitting the relation . Such an assumption is justified by the measurements of the perihelion shift of Mercury, which correspond to a very small value of η, and determine A = 1 at high significance. In the last column of Table 9 we list the ratio between the measured and predicted GR apsidal motion rates, which we plot as a function of η in Fig. 10. The black dashed line and the grey shadow region represent the best fit of B and the 1σ uncertainty. We obtain B = 0.01 ± 0.05, which again is fully consistent with the predicted null value from GR. Only one system, KW Hya, has a value of η that is significantly different from the bulk of the sample at η = 0.250. This clustering of measurements at a small range in η limits the effectiveness of the sample at placing strong constraints on the slope B. For this reason, we cannot perform a fit varying simultaneously A and B from Eq. (8). This would be enabled by considering systems with very unequal component masses, which are difficult to discover and measure with high precision.
Fig. 10. Ratio of measured and predicted GR apsidal motion rates as a function of the dimensionless reduced mass η. The black dashed line and grey shadow area correspond to the best fit of Eq. (8), assuming γ = β = 1, and its 1σ uncertainty. We highlight the broken xaxis and the different scales. The colour code is chosen such that redder colours correspond to systems with the smaller relative error, and therefore are dominating the fit. 
Our analysis provides weaker constraints to PPN parameters compared to for example those from the perhielion shift of Mercury (Verma et al. 2014), the Shapiro timedelay measured by the Cassini mission (Bertotti et al. 2003), or the spin precession of millisecond pulsars (Shao et al. 2013). Nevertheless, the results obtained from the eclipsing binaries analysed in this work perform a test to gravity in a yet unexplored regime of the potentialcurvature diagram. Figure 11 reproduces this latter diagram from Baker et al. (2015) where we add the systems studied here. Stellar eclipsing binaries are one of the many probes of gravity in different environments, in this case covering the range between Solar System planets and binary pulsars.
Fig. 11. Adaptation of Fig. 1 of Baker et al. (2015) showing the parameter space for gravitational fields, including the regime probed by different astrophysical and cosmological systems. We add the regime tested by the eclipsing binaries in this work as a yellow square. See Baker et al. (2015) for more information about the labels. 
7. Conclusions
The goal of the present study is to employ precise and broad photometric coverage provided by the TESS mission to determine precise eclipse timings of a number of binaries and determine their apsidal motion rates. We carefully selected a sample of eccentric eclipsing binaries with accurately determined general properties such as masses and radii, and relatively long orbital periods to ensure a dominant contribution of the GR term over the classical term to the apsidal motion rate. The continuous monitoring and the excellent precision of the TESS new times of minima allowed us to derive determinations of T_{2} − T_{1} and therefore investigate possible variations. From our analysis we were able to detect apsidal motion in 9 out of the 15 systems included in our sample. In some cases, we improved previous determinations, but for 5 of the systems (KX Cnc, AL Dor, V501 Her, KW Hya, and EY Cep) we present the first measurement of apsidal motion. Furthermore, we were able to determine apsidal motion in one system, AL Dor, using data from TESS alone, which cover a very small fraction of the apsidal motion orbital cycle, thus highlighting the excellent performance of the mission for this kind of study.
A comparison between observed apsidal motion rates and theoretically predicted values yields excellent agreement. The good precision of the observed apsidal motion rates, the accurate stellar properties employed, and the relatively modest contribution from the classical term made it possible to calculate reliable estimates of the observed GR term. In some cases, when the classical contribution is especially insignificant, the GR apsidal motion could be measured with a relative precision approaching 1 per cent. This allows a test of GR for the first time using this method. Our results strongly favour the predictions of GR, with no deviations observed at the level of 10^{−2}. Furthermore, we set constraints on two of the PPN parameters. These latter are not quite as stringent as those resulting from other methods (which can be three orders of magnitude more constraining; Will 2014) yet they probe a regime of gravitational forces and potentials that has not been explored before (Baker et al. 2015). Given the excellent agreement with GR, the parameter space for alternative theories continues to narrow.
The results we present are mildly modeldependent as they may be affected by deviations of the k_{2} values. However, because of the small contribution of the classical term, variations of the k_{2} values within their quoted uncertainties induce relative variations in for the two targets dominating our fits, V541 Cyg and AL Dor, of lower than 1%.
With new TESS data still to come, observing new systems and increasing the data coverage of those already observed, promising systems still without a detection of apsidal motion could be added to further extend our study. In order to perform even more stringent tests of gravitation theories, an increase in the number of systems with a GR apsidal motion relative contribution larger than the threshold used in this work is needed. We therefore encourage spectroscopic monitoring and precise determination of absolute properties of longperiod, eccentric eclipsing binaries with TESS photometric data.
Acknowledgments
Raül Vera (EHU, Spain) is gratefully acknowledged for his guidance and discussions regarding the PPN formalism and observational parameter constraints. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018098153BC33, the support of the Generalitat de Catalunya/CERCA programme, and the Agència de Gestió d’Ajuts Universitaris i de Recerca of the Generalitat de Catalunya, with additional funding from the European FEDER/ERF funds, L’FSE inverteix en el teu futur. This work has been carried out within the framework of the PhD programme in Physics of the Universitat Autònoma de Barcelona.
References
 Agerer, F., Busch, H., Kleikamp, W., & Moschner, W. 1994, Berliner Arbeitsgemeinschaft fuer Veraenderliche Sterne  Mitteilungen, 69, 1 [Google Scholar]
 Aigrain, S., Parviainen, H., & Pope, B. J. S. 2016, MNRAS, 459, 2408 [NASA ADS] [Google Scholar]
 Albrecht, S., Reffert, S., Snellen, I. A. G., & Winn, J. N. 2009, Nature, 461, 373 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Albrecht, S., Setiawan, J., Torres, G., Fabrycky, D. C., & Winn, J. N. 2013, ApJ, 767, 32 [CrossRef] [Google Scholar]
 Andersen, J. 1991, A&ARv, 3, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Andersen, J., & Vaz, L. P. R. 1984, A&A, 130, 102 [NASA ADS] [Google Scholar]
 Andersen, J., Clausen, J. V., & Nordstrom, B. 1987, A&A, 175, 60 [Google Scholar]
 Baker, T., Psaltis, D., & Skordis, C. 2015, ApJ, 802, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Bakış, V., Bakış, H., Demircan, O., & Eker, Z. 2008, MNRAS, 384, 1657 [Google Scholar]
 Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Breinhorst, R. A., Pfleiderer, J., Reinhardt, M., & Karimie, M. T. 1973, A&A, 22, 239 [NASA ADS] [Google Scholar]
 Claret, A. 1997, A&A, 327, 11 [NASA ADS] [Google Scholar]
 Claret, A., & Giménez, A. 1993, A&A, 277, 487 [Google Scholar]
 Claret, A., & Giménez, A. 2010, A&A, 519, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A., & Torres, G. 2017, ApJ, 849, 18 [Google Scholar]
 Claret, A., & Torres, G. 2018, ApJ, 859, 100 [Google Scholar]
 Claret, A., & Torres, G. 2019, ApJ, 876, 134 [Google Scholar]
 Claret, A., Torres, G., & Wolf, M. 2010, A&A, 515, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dariush, A., Mosleh, M., & Dariush, D. 2006, Ap&SS, 305, 85 [Google Scholar]
 Davies, D. 2007, Peremennye Zvezdy Prilozhenie, 7, 16 [Google Scholar]
 De Laurentis, M., De Rosa, R., Garufi, F., & Milano, L. 2012, MNRAS, 424, 2371 [Google Scholar]
 Diethelm, R. 1992, Bulletin der BedeckungsveraenderlichenBeobachter der Schweizerischen. Astronomischen Gesellschaft, 99, 10 [Google Scholar]
 Estabrook, F. B. 1969, ApJ, 158, 81 [Google Scholar]
 Feiden, G. A., & Chaboyer, B. 2012, ApJ, 757, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Feinstein, A. D., Montet, B. T., ForemanMackey, D., et al. 2019, PASP, 131 [Google Scholar]
 Ferrero, G., Gamen, R., Benvenuto, O., & FernándezLajús, E. 2013, MNRAS, 433, 1300 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D. 2015, George: Gaussian Process regression [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallenne, A., Pietrzyński, G., Graczyk, D., et al. 2019, A&A, 632, A31 [EDP Sciences] [Google Scholar]
 Gibson, N. P., Pont, F., & Aigrain, S. 2011, MNRAS, 411, 2199 [NASA ADS] [CrossRef] [Google Scholar]
 Giménez, A. 1985, ApJ, 297, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Giménez, A. 2007, in Binary Stars as Critical Tools & Tests in Contemporary Astrophysics, eds. W. I. Hartkopf, P. Harmanec, & E. F. Guinan, IAU Symp., 240, 290 [Google Scholar]
 Giménez, A., & Bastero, M. 1995, Ap&SS, 226, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Giménez, A., Kim, C.H., & Nha, I.S. 1987, MNRAS, 224, 543 [Google Scholar]
 Graczyk, D., Pietrzyński, G., Gieren, W., et al. 2019, ApJ, 872, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Guinan, E. F., & Maloney, F. P. 1985, AJ, 90, 1519 [NASA ADS] [CrossRef] [Google Scholar]
 Guinan, E. F., Maley, J. A., & Marshall, J. J. 1996, IBVS, 4362, 1 [Google Scholar]
 Harmanec, P., Holmgren, D. E., Wolf, M., et al. 2014, A&A, 563, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hejlesen, P. M. 1987, A&AS, 69, 251 [NASA ADS] [Google Scholar]
 Higl, J., & Weiss, A. 2017, A&A, 608, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hubscher, J. 2015, IBVS, 6152, 1 [Google Scholar]
 Hubscher, J., & Lehmann, P. B. 2015, IBVS, 6149, 1 [Google Scholar]
 Hubscher, J., Paschke, A., & Walter, F. 2005, Inf. Bull. Var. Stars, 5657, 1 [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Software and Cyberinfrastructure for Astronomy IV, SPIE Conf. Ser., 9913, 99133E [Google Scholar]
 Karpowicz, M. 1961, Acta Astron., 11, 51 [NASA ADS] [Google Scholar]
 Khaliullin, K. F. 1985, ApJ, 299, 668 [Google Scholar]
 Khaliullin, K. F., & Khaliullina, A. I. 2007, MNRAS, 382, 356 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, C. H., Kreiner, J. M., Zakrzewski, B., et al. 2018, ApJS, 235, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Kopal, Z. 1959, Close binary systems (New York: John Wiley & Sons Inc.) [Google Scholar]
 Kopal, Z. 1978, Dynamics of Close Binary Systems (Netherlands: Springer) [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
 Kwee, K. K., & van Woerden, H. 1956, Bull. Astron. Inst. Netherlands, 12, 327 [Google Scholar]
 Lacy, C. H. S. 1998, AJ, 115, 801 [Google Scholar]
 Lacy, C. H. S. 2002, Inf. Bull. Var. Stars, 5357, 1 [Google Scholar]
 Lacy, C. H. S. 2004, Inf. Bull. Var. Stars, 5577, 1 [Google Scholar]
 Lacy, C. H. S., Hood, B., & Straughn, A. 2001, Inf. Bull. Var. Stars, 5067, 1 [Google Scholar]
 Lacy, C. H. S., Claret, A., & Sabby, J. A. 2004, AJ, 128, 1340 [NASA ADS] [CrossRef] [Google Scholar]
 Lacy, C. H. S., Torres, G., Claret, A., & Vaz, L. P. R. 2005, AJ, 130, 2838 [NASA ADS] [CrossRef] [Google Scholar]
 Lacy, C. H. S., Torres, G., Claret, A., & Menke, J. L. 2006, AJ, 131, 2664 [Google Scholar]
 Lastennet, E., & VallsGabaud, D. 2002, A&A, 396, 551 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LeviCivita, T. 1937, Am. J. Math., 59, 225 [Google Scholar]
 Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [Google Scholar]
 Mikulášek, Z., Chrastina, M., Liška, J., et al. 2014, Contrib. Astron. Obs. Skalnate Pleso, 43, 382 [Google Scholar]
 Milone, E. F., KurpińskaWiniarska, M., & Oblak, E. 2010, AJ, 140, 129 [Google Scholar]
 Moffat, J. W. 1986, Can. J. Phys., 64, 178 [Google Scholar]
 Naoz, S. 2016, ARA&A, 54, 441 [Google Scholar]
 Nelson, R. H. 2003, Inf. Bull. Var. Stars, 5371, 1 [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Pols, O. R., Tout, C. A., Schroder, K.P., Eggleton, P. P., & Manners, J. 1997, MNRAS, 289, 869 [Google Scholar]
 Rauw, G., Rosu, S., Noels, A., et al. 2016, A&A, 594, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ribas, I., Jordi, C., & Giménez, Á. 2000, MNRAS, 318, L55 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1 [Google Scholar]
 Rosu, S., Noels, A., Dupret, M. A., et al. 2020, A&A, 642, A221 [EDP Sciences] [Google Scholar]
 Sandberg Lacy, C. H., & Fekel, F. C. 2014, AJ, 148, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Sandberg Lacy, C. H., Zakirov, M., Arzumanyants, G., et al. 1995, IBVS, 4194, 1 [Google Scholar]
 Schmitt, J. H. M. M., Schröder, K. P., Rauw, G., et al. 2016, A&A, 586, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I. 1985, Sov. Astron. Lett., 11, 224 [NASA ADS] [Google Scholar]
 Shao, L., Caballero, R. N., Kramer, M., et al. 2013, CQG, 30 [Google Scholar]
 Smith, A. B., & Caton, D. B. 2007, IBVS, 5745, 1 [Google Scholar]
 Southworth, J. 2015, in Living Together: Planets, Host Stars and Binaries, eds. S. M. Rucinski, G. Torres, & M. Zejda, ASP Conf. Ser., 496, 164 [Google Scholar]
 Sowell, J. R., Henry, G. W., & Fekel, F. C. 2012, AJ, 143, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Sterne, T. E. 1939, MNRAS, 99, 662 [NASA ADS] [Google Scholar]
 Tkachenko, A., Pavlovski, K., Johnston, C., et al. 2020, A&A, 637, A60 [CrossRef] [EDP Sciences] [Google Scholar]
 Tomkin, J., & Fekel, F. C. 2006, AJ, 131, 2652 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., & Ribas, I. 2002, ApJ, 567, 1140 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., Lacy, C. H. S., Claret, A., & Sabby, J. A. 2000, AJ, 120, 3226 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., Sandberg Lacy, C. H., & Claret, A. 2009, AJ, 138, 1622 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., Andersen, J., & Giménez, A. 2010, A&ARv, 18, 67 [Google Scholar]
 Torres, G., Sandberg Lacy, C. H., Pavlovski, K., et al. 2014, ApJ, 797, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., Sandberg Lacy, C. H., Pavlovski, K., Fekel, F. C., & Muterspaugh, M. W. 2015, AJ, 150, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., McGruder, C. D., Siverd, R. J., et al. 2017, ApJ, 836, 177 [Google Scholar]
 Verma, A. K., Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2014, A&A, 561, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Volkov, I. M., & Khaliullin, K. F. 1999, IBVS, 4680, 1 [Google Scholar]
 Will, C. M. 2014, Liv. Rev. Relativ., 17, 4 [Google Scholar]
 Wolf, M., Kučáková, H., Kolasa, M., et al. 2006, A&A, 456, 1077 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolf, M., Claret, A., Kotková, L., et al. 2010, A&A, 509, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zasche, P., & Wolf, M. 2019, AJ, 157, 87 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Properties of 15 eclipsing binary systems with eccentric orbits, accurate absolute dimensions, and literaturepredicted relativistic apsidal motion relative contribution greater than 60%.
TESS sectors and cadence at which the sample of studied eclipsing binary systems have been observed, the timespan between the first and last eclipse included in this work, and their Gaia G magnitudes.
log k_{2} values and apsidal motion rates, both observed and theoretically predicted, for the eclipsing binaries with apsidal motion measurement.
All Figures
Fig. 1. Examples of TESS data for an eclipse of V459 Cas (a) with 2min cadence, and V501 Mon (c) with 30min cadence. In panels (b) and (d) the black line shows the bisector corresponding to each eclipse. The computed time of minimum for each system, which is used as reference time in the plots, is marked as a blackdashed line in panels a and c. Its uncertainty is illustrated in panels b and d as a read shaded area. The resulting timing uncertainty is 3.3 s for V459 Cas and 22.0 s for V501 Mon. 

In the text 
Fig. 2. Ephemeris curve for AL Dor as a function of the orbital cycle. The dotted and dashed lines represent fits to the primary (blue circles) and secondary (red triangles) eclipses, respectively. For a better visualisation, an arbitrary shift to the primary and secondary minima has been applied. 

In the text 
Fig. 3. T_{2} − T_{1} as a function of the orbital cycle for V541 Cyg. 

In the text 
Fig. 4. T_{2} − T_{1} as a function of the orbital cycle for V459 Cas. 

In the text 
Fig. 5. Same as Fig. 2 for V501 Mon. 

In the text 
Fig. 6. Same as Fig. 2 but for GG Ori. 

In the text 
Fig. 7. T_{2} − T_{1} as a function of the orbital cycle for EY Cep. 

In the text 
Fig. 8. Comparison between total computed and total observed apsidal motion rates. The black dashed line indicates the 1:1 relation. 

In the text 
Fig. 9. Measured GR apsidal motion as a function of the theoretically computed GR apsidal motion. The black dashed line and grey shadow area correspond to the best fit to Eq. (8), assuming α_{i} = ζ_{2} = 0, and its 1σ uncertainty. The colour code is chosen such that redder colours correspond to systems with the smaller relative error, and therefore are dominating the fit. 

In the text 
Fig. 10. Ratio of measured and predicted GR apsidal motion rates as a function of the dimensionless reduced mass η. The black dashed line and grey shadow area correspond to the best fit of Eq. (8), assuming γ = β = 1, and its 1σ uncertainty. We highlight the broken xaxis and the different scales. The colour code is chosen such that redder colours correspond to systems with the smaller relative error, and therefore are dominating the fit. 

In the text 
Fig. 11. Adaptation of Fig. 1 of Baker et al. (2015) showing the parameter space for gravitational fields, including the regime probed by different astrophysical and cosmological systems. We add the regime tested by the eclipsing binaries in this work as a yellow square. See Baker et al. (2015) for more information about the labels. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.