Modeling Apsidal Motion in Eclipsing Binaries using ELC

Apsidal motion is the precession of the line of apsides in the orbit of a binary star due to perturbations from General Relativity (GR), tides, or third-body interactions. The rate of precession due to tidal effects depends on the interior structures of the stars, and as a result, binaries in which this precession occurs are of great interest. Apsidal motion is observed through the analysis of eclipse times, which reveal small changes in the average interval between successive primary and secondary eclipses, taking all available observed times of eclipse and yielding an estimate of the apsidal rate. Given that this is a single observed quantity, various degeneracies are unavoidably present. Ideally, one would have a model that predicts eclipse times given the orbital and stellar parameters. These parameters for a given binary could then be computed using least squares, provided a suitably large number of eclipse times. Here we use the eclipsing light curve (ELC) program as such a model. The Newtonian equations of motion with additional force terms accounting for GR contributions and tidal distortions are integrated, yielding precise sky positions as a function of time. Times of mid-eclipse and instantaneous orbital elements are computed as a function of time. In this paper, we outline the method and compare numerically computed apsidal rates with standard formulae using a set of 15 binaries based on real systems. For our simulated systems, the derived apsidal rates agree with the standard formula.


Introduction
Apsidal motion is the precession of the line of apsides in an orbit, generally in the same direction as the motion of the secondary body. For a binary with an eccentric orbit, the phase difference between the primary and secondary eclipse depends on the argument of periastron ω. If ω changes over time, then the primary and secondary eclipse times will not follow simple linear ephemerides. When the primary and secondary eclipse times are fit to a common ephemeris (the free parameters to this fit are the common period, the common reference time of primary eclipse, and the average time difference between the primary and secondary eclipses), the residuals in the common period observed minus computed (CPOC, or O − C) diagram are roughly sinusoidal with opposite phases (see Figure 1 for a schematic example for a binary resembling Y Cyg). When eclipse times are available over a long enough time interval, then the apsidal period (e.g., the time it takes for ω to change by 360°) can be found from the CPOC diagram.
Mathematically, the apsidal motion can be computed as a result of perturbations of a Keplerian orbit that include a General Relativity (GR) contribution, tidal and rotational contributions, and possibly a third-body contribution. The combination of these results in the expression of the rate of change of the argument of periastron, Analytic expressions exist with which the various w  components can be computed. These expressions may depend on the orbital parameters (e.g., the mean period, the mean eccentricity, and the component masses) and on the stellar parameters (e.g., the stellar radii and the internal density profiles). Depending on the situation, GR w  may dominate, in which case the binary could be a useful test of GR. In other situations, N GR w w    , in which case the rate of apsidal motion can be useful as a probe of the internal structure of the stars.
Generally speaking, in previous studies of apsidal motion, all of the available eclipse times for a given binary have been condensed into a single measurement of w  (e.g., Gimenez & Garcia-Pelayo 1983). For binaries where N w  dominates, this leads to unavoidable degeneracies because each star contributes to the apsidal motion. Nevertheless, useful information has been found from binaries with roughly equal components, where the hope is that each star has similar interior structures (see the review paper by Torres et al. 2010).
In the ideal case, a model would predict the times of primary and secondary eclipses given the orbital and stellar parameters, which would include a parameterization of the internal density profiles. Then, for a given binary, the observed times of the primary and secondary eclipses (with their uncertainties) can be modeled, and the best-fitting orbital and stellar parameters can be found using standard least-squares techniques. In this paper, we discuss how we use the eclipsing light curve (ELC) program (Orosz & Hauschildt 2000) as such a model. The Newtonian equations of motion with additional force terms to account for GR contributions and tidal distortions are integrated, yielding precise sky positions as a function of time. Times of mid-eclipse and instantaneous orbital elements can then be computed as a function of time.
This paper is organized in the following manner: In Section 2 we discuss the analytic formula used to find w  for the GR, tidal, and third-body cases. In Section 3 we discuss how the ELC code is used to find eclipse times given the orbital and stellar parameters. This model can also be used to find the apsidal period for a given binary, and in Section 5 we validate ⎠ in units of degrees per cycle, where m 1 and m 2 are the masses of the binary components in solar masses, and P d is the period in days. The combination GM e is known more accurately than either G or M e are alone, where A = 149597870700.0 is the number of meters in an Astronomical Unit (an exact number by definition), and k = 0.01720209895 radians per day is the Gaussian gravitational constant (Clemence 1965). We finally arrive at the wellknown formula for GR w  (e.g., Gimenez 1985), The coefficient in this approximation is an exact expression from fundamental constants. From the formula, we see that systems with higher stellar masses and shorter periods will have faster rates of apsidal motion due to GR effects.

Tidal (Newtonian) Contribution
In Newtonian gravity, the orbit of two bound point masses is given by the well-known Kepler equations. The orbit is closed, and the orientation of the semimajor axis (the line of apsides) remains fixed. Real stars are not point masses, and departures from spherical symmetry due to rotation and/or tides give rise to small nonradial forces that cause the orbital elements to change with time (see the Lagrange planetary equations). In most cases, the rate of change of the line of apsides (characterized by the so-called argument of periastron ω) is highest and therefore most readily observable.
where n 2 a 3 = G(m 1 + m 2 ) and χ = − nT, with n = 2π/P being the mean daily motion. The eccentricity, nodal angle, and inclination are defined by e, Ω N , and i, respectively. These expressions and the functional form of S therein determine which of the orbital elements will change over time, and how quickly they change.
Expressions for N w due to tidal distortions were first derived by Cowling (1938) and Sterne (1939). Following the revised formulation by Kopal (1978), the expression for the apsidal period U from tides is given by where the k 2i factors are the second-order (quadrupolar) internal structure constants. These internal structure constants are related to the density distribution within the star, and they can be computed from stellar evolution models.
The weighting coefficients c i are functions of the mass ratio, eccentricity, and relative radii with respect to the orbital separation, and they are given by where r i is the fractional radius scaled to the orbital separation (R i /a), the parameters γ i are the ratio of the angular velocity of the axial rotation to that of orbital motion, and g(e) and f (e) are functions of the eccentricity, originally estimated as a power series in e by Sterne (1939), and provided by Bulut et al. (2017) When we adopt pseudo-synchronous rotation in eccentric orbits, as done in Hut (1981), the maximum angular velocity at periastron can be well approximated as e e 1 1 .
A similar recent approach by Bulut et al. (2017) determines the c i coefficients through a combination of the contributions of rotational distortion and tidal effects, In the expression for c i (Equation (19)), r i is the fractional radius, m i is the mass, Ω r is the angular velocity of the axial rotation for each component i, Ω K is the Keplerian angular velocity, and e is the orbital eccentricity. Thus, when the Keplerian parameters are known, then the stellar parameters including the rotation rates and the apsidal period can be calculated.
The mean value of the internal structure constants can be derived from the observed value of w  using the expression where the c 2i coefficients are the same functions of the eccentricity, mass, radius, and separation as in Equation (19).
It is known that the observed mean value of the internal structure constants contains both contributions from the Newtonian and the relativistic effects of apsidal motion. When the constants are combined through the equation the weighted average coefficient k 2,theō can be determined. This weighted average is directly comparable with the observed value. However, this historical method fails to constrain the individual constants k 21 and k 22 as only the average value is returned, and this method does not work well with binaries with mass ratios q 1 » .

Stellar Spin Axes
We can break down the Newtonian component of apsidal motion into the contributions from tides and rotation: If the spin axes are misaligned, there is an additional contribution to the tidal term in the expression for apsidal motion. The pointing direction of the angular momentum vector of the stars will affect the apsidal motion of the system. This is parameterized in two values, where ∠α i is the deflection in the plane of the orbit, and ∠β i is the deflection in the plane of the sky. The combination of the two of these can result in any pointing direction for either the primary or secondary star.
A summary of the theory behind misaligned spin axes and the subsequent affect on the secular motion of the apse is given in Shakura (1985). In particular, his Equation (4) where the (dω/dt) E is the GR contribution. If the spin axes of the stars are aligned, then α 1 = α 2 = 0 and β 1 = β 2 = i, where i is the inclination of the orbital plane of the binary.

Third-body Contribution
Although we focus our attention in this paper on binaries with apsidal motion due to tidal or GR effects, for completeness, we mention the possible contribution from a third body, which can lead to observable changes in the CPOC diagram. The addition caused by the third body has to be accounted for before the apsidal motion can be analyzed.
Because the binary orbits the barycenter of the triple system, the observed times of the primary and secondary eclipses can either be early or late because the distance between the observer and the binary changes periodically. The signal in the O − C diagram then superficially resembles a radial velocity curve. Irwin (1952) gives a formula to model this light travel time effect (LTTE) signal. In the massive binary DR Vul, the CPOC signals have been modeled by a combination of an LTTE orbit with an orbital period of ≈63 yr and apsidal motion of the inner binary with an apsidal period ≈36 yr (Wolf et al. 2019;Dimoff 2021).
If the third body is sufficiently close to the binary, the gravitational perturbations can lead to changes in the orbital period as well as in the precession of the orbit, and these changes can even dwarf the LTTE. Several binaries with large eclipse-timing variations (ETVs) due to a third body have been discovered using data from the Kepler and Transiting Exoplanet Survey Satellite (TESS) missions. Generally speaking, the ETV signal seen in the O − C diagram from dynamical interactions is complex. Borkovits et al. (2011) and Baycroft et al. (2023;their Appendix A) give expressions for w  including effects from a third body. Although these dynamical interactions can produce measurable effects in the O − C on both short and long timescales, the short-period low-amplitude variations are less important for apsidal motion studies. The long-period perturbations in the apsidal motion, however, can substantially alter the tidal and relativistic effects (see Naoz 2016, or Borkovits et al. 2020 for a recent example).

Modeling Eclipse Times Using ELC
We now discuss our forward model, which can produce the times of primary and secondary eclipses given the stellar parameters and initial orbital parameters. It is based on the ELC code (Orosz & Hauschildt 2000;Orosz et al. 2019). The code is general, and the light and velocity curves of a variety of binary and three-body systems can be modeled directly. In the mode that we describe here, the observed times of eclipse can be fitted, which is useful in situations without access to the light curves.
The basic outline of a photodynamical code like ELC is relatively straightforward. Given the masses of two (or more) bodies and initial conditions (positions and velocities) for these bodies, the equations of motion are integrated, yielding the sky positions as a function of time. From the time series of the positions, it is easy to find the times when any two bodies are at conjunction. It can easily be checked if an eclipse occurs at or near conjunction when given the radius of each body.
For convenience, ELC has a Keplerian-to-Cartesian converter (based on the algorithms given in Murray & Dermott 1999), where six orbital elements (the period P, the time of periastron passage T, the eccentricity e, the argument of periastron ω, the inclination i, and the nodal angle Ω) uniquely determine six phase-space coordinates for each star (e.g., x, y, z, v x , v y , and y z ). The x, y plane is the sky plane, and the z-axis points toward the observer. The inverse transformation (the Cartesian-to-Keplerian converter) is available, where six phasespace coordinates give a unique set of orbital elements.
The equations of motion are the usual force equations for the two-body problems, plus additional force terms that arise from tidal distortions and force terms from the full GR treatment (Mardling & Lin 2002), The acceleration due to the quadrupole moment of body 1 is a combination of its spin distortion and tidal distortion produced by the presence of the companion, where k 21 is the apsidal motion constant of body 1, Ω 1 is the ratio of the rate of stellar rotation and the orbital rotation, and r is a unit vector in the direction of r. A similar expression exists for the tidal distortion of body 2 due to body 1. The orbital acceleration due to the relativistic potential of the binary (Kidder 1995) is given by To find the positions of all the bodies at any given time, the equations of motion are integrated using a symplectic twelfthorder Gaussian Runge-Kutta (GRK) integration routine based on the methods from Hairer et al. (2006). In our implementation, one can solve only the purely Newtonian equations (the forces are only from point masses), the Newtonian equations plus the tidal forces, the Newtonian equations plus the GR forces, or the equations with all three contributions. As noted previously, if the additional force terms due to tides or due to GR effects are included, the orbit will not be closed, and the orbital elements will change with time. The Cartesian-to-Keplerian converter can be used to compute time series of the orbital elements. Of interest here is the time series for the argument of periastron ω(t).

Comparing Derived Apsidal Rates with Analytic Formula
The analytic formula discussed in Section 2 gives the apsidal rate, namely w  . For a model that can produce a time series for ω(t), the rate of change of ω (e.g., w  ) should be easy to compute and can be compared to the analytic formula. However, in this case, there are some subtle issues. Figure 2 displays the argument of periastron (ω) over time for a system resembling Y Cyg considering tides and GR independently. For this calculation, ω 0 = 1°. The system is precessing, and over the course of 1000 days, ω changed by about 20°due to tidal effects. However, when the graph is examined more closely, small small-scale oscillations become visible with a period equal to the binary orbital period. The amplitude of these oscillations depends on the apsidal period, and in the limit where there is no apsidal motion (e.g., when 0 w =  ), the oscillations must vanish. At each time step, six phase-space coordinates give a unique set of orbital parameters (including ω).
However, because the underlying motion is not Keplerian, higher-order oscillations are noticeable. The equations of motion with the corrections provide an orbit-averaged force. The effect of that force depends on the instantaneous separation of the bodies, so the actual speed of the star changes in a slightly different manner relative to the average Keplerian. Given these oscillations, the value of w  from a given model is computed by using a linear fit to the ω(t) curve. More reliable values of omega dot are obtained from averaging over a longer time-span that contains many orbital periods.
From a collection of eclipsing binaries exhibiting rapid apsidal motion (Claret & Willems 2002), we select 15 prototypical binary systems to use as a test set for our model. We simulate the orbital motion of these systems each for a set number of days relative to an apsidal cycle, taking into account their physical parameters including the mass, radius, apsidal constants, and orbital spin-axis orientations of both components. A summary of the input parameters for each system is presented in Table 1 (for assumed rotationally aligned systems). We do not claim that these are the actual physical parameters of the selected systems, but they are similar enough to their physical values to represent realistic binaries.
For each of the 15 binaries, we ran 200,000 forward models for three situations: (i) GR perturbations only, (ii) tidal perturbations with aligned spin axes, and (iii) tidal perturbations with misaligned spin axes. The initial value of ω was always 1°, and the values of the other orbital elements and the stellar masses and radii were drawn from normal distributions that represent typical measurement uncertainties on these quantities. The model values of w  for each binary system are determined from a linear regression of ω(t), and are compiled in Table 2. We compare these values to results from the respective analytic formula as a fractional percent difference, and discuss each of these three cases in turn.

Accuracy of the General Relativity Contribution
A histogram of the error distributions for the GR case for each of our selected 15 systems is shown in Figure 3, limited to the range −0.002% to 0.002%. While the combined distributions are not centered exactly on zero, each individual distribution is approximately symmetric. Systems with lower eccentricity (e.g., CW Cep, EM Car, V478 Cyg, or U Oph) exhibit wider histograms, as well as peaks or horns close to the ends of the distributions. This is the result of the nearambiguity in the value of the argument of periastron ω in a nearly circular orbit.
When only GR is taken into account, the model apsidal rates agree quite well with the formula given by Equation (7). The percent errors are ∼5 × 10 −4 on average. We note at this level, there are not enough digits in the coefficients given in Equation (7) from Gimenez (1985) for a meaningful comparison with our models.  (2002)

Accuracy of the Tidal Contribution
Histograms of the error distributions for each of our selected 15 systems considering the tidal contribution are shown in Figure 4 for the aligned case and in Figure 5 for the misaligned case. The analytic approach to modeling tides taken by Sterne (1939;e.g., Equation (14)) works well. In this case, we find that the relative errors are acceptably small, within 0.4%. For the misaligned case (Equation (22)), the percent errors are larger, ∼2%.

Discussion
Our implementation of GR is very accurate when compared to the analytic formula (Equation (7)). The differences are small but systematically positive, where the median percent difference is ≈0.00067. For comparison, Weinberg (1972) gives the rate of precession as MG L 6 radians revolution , 27 where L = (1 − e 2 )a is a measure of the eccentricity of the orbit, the so-called semilatus rectum. This formula is exact in    the Schwarzschild metric, and agrees with Barker & O'Connell (1975b) when appropriate substitutions are made. The agreement between Equation (5) and our model indicates that the force equations account quite well for the GR perturbations.
While not as good as GR, our modeled apsidal rates for the tidal contribution with an aligned stellar spin axes are accurate compared to the analytic formula (Equation (14)). Again, the differences are systematically positive, with a median percent difference value of ≈0.013. Unlike the GR case, there is an approximation in the formula for w  . Sterne (1939) uses the n = 2 (second-order) term in the tidal potential, which is proportional to (r 2 /a 3 ). The additional force terms in our model are presumably derived from a treatment that uses the n = 2 term in the disturbing potential. Apparently, in this case, the solutions to the differential equations do not reproduce the analytic results as faithfully as they do in the GR case. According to Sterne (1939), when the third-order term is included in the tidal potential, w  scales with (r/a) 7 . For the binaries considered here, (r/a) 7 is only a few percent of (r/a) 5 in Equation (19).
When the axial misalignment of the stars is taken into account, the accuracy is not as good, although the histograms of the percent differences are closer to symmetric about zero. For each star, we establish an axis deflection α i as the spinorbit inclination, and axis β i as the inclination with respect to the plane of the sky (Equation (22)). These axial inclination parameters for each star can define any spin-axis orientation for a binary. For a system where the spin axes of the stars are aligned with the angular momentum axis of the orbit to within about 15°, Equation (22) is accurate to within a few percent.
As seen in Figure 6, the largest relative errors occur at mutual inclinations of ∼45°, representing the maximum allowed asymmetry as the tidal bulge of the companion changes hemispheres. Furthermore, the relative errors are minimized at mutual inclinations of 0°and 90°, where in the former case, the spin axis is perpendicular to the plane of the orbit, and in the latter case, the spin axis is within the plane of the orbit. One simulated system each of these regimes is plotted in Figure 7.
We fit a linear regression to the line of ω versus t and compute the residuals. In systems close to the aligned case (mutual inclination ≈0), a linear trend is recovered, visible in the top left panel of Figure 7. For extreme cases, the ω(t) is highly nonlinear. In these extreme cases, where the mutual inclination 0 » (where the computed relative error compared to the analytic formula is large), we find a second-order trend with the ω(t) curve. This can be seen in the top right panel of Figure 7. Mardling & Lin (2002) point out that if the spin axis of one of the stars is not aligned, there will be a torque on that star. If there is a torque on that star, then the direction of the spin axis will change ( r f r QD W µ ). The expression for the force given in Equation (24) contains Ω r , so that in principle, if that vector changes with time, it could add complications to the system. However, our numerical experiments show that r W  is very small in most practical cases. Furthermore, Shakura (1985) shows that the inclination i and nodal angle Ω N also change Figure 7. Residuals from a linear fit to the advance of the argument of periastron in two cases of axial misalignment for a Y Cyg-like system. The left panels represents a system close to alignment; the axial deflections result in a slightly different pattern of the oscillations, but maintain a linear trend. The right panels represent a system far from axial alignment; the oscillations have more structure, and the overall trend is nonlinear.
with time, and gives corresponding formulae for di/dt and dΩ N /dt. Given the changes in the other orbital elements, it is not surprising that the ω versus t curve becomes nonlinear on a shorter timescale. Indeed, Equation (22) for dω/dt includes a term for the inclination. If di/dt is nonzero, there are additional contributions to dω/dt.

Conclusion
The precession of the line of apsides in a binary orbit is driven by perturbations from GR, tidal, and third-body interactions. These perturbations depend on the physical and orbital properties of the systems, including the orbital period, eccentricity, stellar mass, internal structure constants, and relative rotation speeds of each star, and it is affected by the angular deflection of the stellar rotation axes. We investigate the effectiveness of the ELC software in modeling the apsidal motion of 15 realistic binary systems. Modifying Newton's equations of orbital motion to include dynamical perturbations from tides, GR, axial misalignment, and possible third-body effects, we run a set of forward models and compute the rates of apsidal motion, and compare them to the corresponding analytical formulae. Our results indicate an extremely good agreement between for the GR contribution, and a good agreement between our modeling and theoretical predictions for the tidal contribution to the apsidal motion. Despite this agreement, inconsistencies remain in the case considering tides with misaligned rotation axes. These functional degeneracies between derived orbital and physical parameters will be addressed in a follow-up study.
This approach to numerically modeling the tidal and GR forces is a fast and precise way to model apsidal motion in binary (or higher-order multiple) stellar systems. We plan to apply this technique to real eclipse-timing data and compare our results in a future work, and we aim to compute refined orbital parameters, internal structure constants, and stellar rotation axis alignments. Note. For all systems except for DI Her, the rate of advance due to GR is much slower than the tidal rate.