Modeling Apsidal Motion in Eclipsing Binaries Using ELC

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.

Zoom In Zoom Out Reset image size

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,

$$\begin{eqnarray}&&\dot{\omega }={\dot{\omega }}_{\mathrm{GR}}+{\dot{\omega }}_{N}+{\dot{\omega }}_{3}.\end{eqnarray} \tag{ 1 }$$

Analytic expressions exist with which the various $\dot{\omega }$ 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, ${\dot{\omega }}_{\mathrm{GR}}$ may dominate, in which case the binary could be a useful test of GR. In other situations, ${\dot{\omega }}_{\mathrm{GR}}\ll {\dot{\omega }}_{N}$, 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 $\dot{\omega }$ (e.g., Gimenez & Garcia-Pelayo 1983). For binaries where ${\dot{\omega }}_{N}$ 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 $\dot{\omega }$ 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 the ELC model by comparing model-computed values of $\dot{\omega }$ to the standard formula. We conclude with a short summary.

2.1. General Relativity Contribution

The rate of apsidal advance in an orbit in GR is given by Equation (10) in Barker & O'Connell (1975a),

$$\begin{eqnarray}&&\dot{\omega }=\displaystyle \frac{3{G}^{2/3}{n}^{5/3}{\left({M}_{1}+{M}_{2}\right)}^{2/3}}{{c}^{2}(1-{e}^{2})},\end{eqnarray} \tag{ 2 }$$

where n = 2π/P is the mean daily motion. Substituting, we obtain

$$\begin{eqnarray}&&\dot{\omega }=\left[\displaystyle \frac{3{G}^{2/3}{\left({M}_{1}+{M}_{2}\right)}^{2/3}}{{c}^{2}(1-{e}^{2})}\right]{\left(\displaystyle \frac{2\pi }{P}\right)}^{5/3},\end{eqnarray} \tag{ 3 }$$

where the units of $\dot{\omega }$ are radians per second. We multiply by the orbital period P to obtain radians per cycle, then by 180/π to obtain degrees per cycle,

$$\begin{eqnarray}&&\dot{\omega }=\left[\displaystyle \frac{540{\left(2\pi \right)}^{5/3}{G}^{2/3}}{\pi {c}^{2}}\right]\left(\displaystyle \frac{1}{1-{e}^{2}}\right){\left(\displaystyle \frac{{M}_{1}+{M}_{2}}{P}\right)}^{2/3}.\end{eqnarray} \tag{ 4 }$$

In the above equation, the units of the masses are kilograms and the units of the period are seconds. To convert into solar masses and days, respectively, we let M_⊙ be the solar mass in kg and s = 86,400 be the number of seconds in a day. We then have

$$\begin{eqnarray}&&\dot{\omega }=\left[\displaystyle \frac{540{\left(2\pi \right)}^{5/3}}{\pi {c}^{2}}{\left(\displaystyle \frac{{{GM}}_{\odot }}{s}\right)}^{2/3}\right]\left(\displaystyle \frac{1}{1-{e}^{2}}\right){\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{{P}_{d}}\right)}^{2/3}\end{eqnarray} \tag{ 5 }$$

in units of degrees per cycle, where m₁ and m₂ are the masses of the binary components in solar masses, and P_d is the period in days. The combination GM_⊙ is known more accurately than either G or M_⊙ are alone,

$$\begin{eqnarray}&&{{GM}}_{\odot }=\displaystyle \frac{{k}^{2}{A}^{3}}{{s}^{2}},\end{eqnarray} \tag{ 6 }$$

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 well-known formula for $\dot{{\omega }_{\mathrm{GR}}}$ (e.g., Gimenez 1985), where the units are in degrees per cycle,

$$\begin{eqnarray}\begin{array}{rcl}\dot{\omega } & = & \left[\displaystyle \frac{540{\left(2\pi \right)}^{5/3}{k}^{2/3}{A}^{2}}{\pi {c}^{2}{s}^{2}}\right]\left(\displaystyle \frac{1}{1-{e}^{2}}\right){\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{{P}_{d}}\right)}^{2/3}\\ & = & \,\left(5.447127276\times {10}^{-4}\right)\left(\displaystyle \frac{1}{1-{e}^{2}}\right){\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{{P}_{d}}\right)}^{2/3}.\end{array}\end{eqnarray} \tag{ 7 }$$

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.

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

From Roy (1978), a form of the Lagrange planetary equations is

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{2}{{na}}\displaystyle \frac{\partial S}{\partial \chi }\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=\displaystyle \frac{1}{{{na}}^{2}e}\left[{\left(1-e\right)}^{2})\displaystyle \frac{\partial S}{\partial \chi }-{\left(1-{e}^{2}\right)}^{1/2}\displaystyle \frac{\partial S}{\partial \omega }\right]\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\chi }{{dt}}=-\displaystyle \frac{(1-{e}^{2})}{{{na}}^{2}e}\displaystyle \frac{\partial S}{\partial e}-\displaystyle \frac{2}{{na}}\displaystyle \frac{\partial S}{\partial a}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Omega }}}_{N}}{{dt}}=\displaystyle \frac{1}{{{na}}^{2}{\left(1-{e}^{2}\right)}^{1/2}\sin i}\displaystyle \frac{\partial S}{\partial i}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\omega }{{dt}}=\displaystyle \frac{{\left(1-{e}^{2}\right)}^{1/2}}{{{na}}^{2}e}\displaystyle \frac{\partial S}{\partial e}-\displaystyle \frac{\cot i}{{{na}}^{2}{\left(1-{e}^{2}\right)}^{1/2}}\displaystyle \frac{\partial S}{\partial i}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{di}}{{dt}}=\displaystyle \frac{1}{{{na}}^{2}{\left(1-{e}^{2}\right)}^{1/2}}\left[\cot i\displaystyle \frac{\partial S}{\partial \omega }-\csc i\displaystyle \frac{\partial S}{\partial {{\rm{\Omega }}}_{N}}\right],\end{eqnarray} \tag{ 13 }$$

where n² a³ = G(m₁ + m₂) 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 $\dot{{\omega }_{N}}$ 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

$$\begin{eqnarray}&&\displaystyle \frac{P}{U}=\displaystyle \frac{{\dot{\omega }}_{N}}{2\pi }={c}_{1}{k}_{21}+{c}_{2}{k}_{22},\end{eqnarray} \tag{ 14 }$$

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

$$\begin{eqnarray}&&{c}_{i}={r}_{i}^{5}\left(\displaystyle \frac{{m}_{3-i}}{{m}_{i}}[15{\boldsymbol{g}}({\boldsymbol{e}})+{\gamma }_{i}^{2}f(e)]+{\gamma }_{i}^{2}{\boldsymbol{f}}({\boldsymbol{e}})\right),\end{eqnarray} \tag{ 15 }$$

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) as

$$\begin{eqnarray}&&f(e)=\displaystyle \frac{1}{{\left(1-{e}^{2}\right)}^{2}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&g(e)=\displaystyle \frac{(8+12{e}^{2}+{e}^{4})f{\left(e\right)}^{2.5}}{8}.\end{eqnarray} \tag{ 17 }$$

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

$$\begin{eqnarray}&&{\gamma }_{i}^{2}=\displaystyle \frac{(1+e)}{{\left(1-e\right)}^{3}}=\displaystyle \frac{{\omega }_{P}^{2}}{{{\rm{\Omega }}}_{{\rm{K}}}^{2}}.\end{eqnarray} \tag{ 18 }$$

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,

$$\begin{eqnarray}&&{c}_{i}={r}_{i}^{5}\left[{\left(\displaystyle \frac{{{\rm{\Omega }}}_{r,i}}{{{\rm{\Omega }}}_{{\rm{K}}}}\right)}^{2}\left(1+\displaystyle \frac{{m}_{3-i}}{{m}_{i}}\right)f(e)+\displaystyle \frac{15{m}_{3-i}}{{m}_{i}}g(e)\right].\end{eqnarray} \tag{ 19 }$$

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 $\dot{\omega }$ using the expression

$$\begin{eqnarray}&&{\bar{k}}_{2,\mathrm{obs}}=\displaystyle \frac{1}{{c}_{21}+{c}_{22}}\displaystyle \frac{P}{U}=\displaystyle \frac{1}{{c}_{21}+{c}_{22}}\displaystyle \frac{\dot{\omega }}{2\pi },\end{eqnarray} \tag{ 20 }$$

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

$$\begin{eqnarray}&&{\bar{k}}_{2,\mathrm{theo}}=\displaystyle \frac{{c}_{21}{k}_{21,\mathrm{theo}}+{c}_{22}{k}_{22,\mathrm{theo}}}{{c}_{21}+{c}_{22}},\end{eqnarray} \tag{ 21 }$$

the weighted average coefficient ${\bar{k}}_{2,\mathrm{theo}}$ can be determined. This weighted average is directly comparable with the observed value. However, this historical method fails to constrain the individual constants k₂₁ and k₂₂ as only the average value is returned, and this method does not work well with binaries with mass ratios $q\rlap{/}{\approx }1$.

2.2.1. Stellar Spin Axes

We can break down the Newtonian component of apsidal motion into the contributions from tides and rotation:

$$\begin{eqnarray*}\begin{array}{rcl}{\dot{\omega }}_{\mathrm{total}} & = & {\dot{\omega }}_{\mathrm{GR}}+{\dot{\omega }}_{\mathrm{tidal},1}+{\dot{\omega }}_{\mathrm{tidal},2}\\ & & +\,{\dot{\omega }}_{\mathrm{rot},1}{\phi }_{1}+{\dot{\omega }}_{\mathrm{rot},2}{\phi }_{2}+{\dot{\omega }}_{\mathrm{LTT}}.\end{array}\end{eqnarray*}$$

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) gives the rate of change of ω,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d\omega }{{dt}}={\left(\displaystyle \frac{d\omega }{{dt}}\right)}_{E}+\bar{\omega }15g(e)\left[{k}_{1}{r}_{1}^{5}\displaystyle \frac{{m}_{2}}{{m}_{1}}+{k}_{2}{r}_{2}^{5}\displaystyle \frac{{m}_{1}}{{m}_{2}}\right]\\ -\,\displaystyle \frac{\bar{\omega }}{{\sin }^{2}i}f(e)\left\{{k}_{1}{r}_{1}^{5}{\left(\displaystyle \frac{{\omega }_{1}}{\bar{\omega }}\right)}^{2}\left(1+\displaystyle \frac{{m}_{2}}{{m}_{1}}\right)\right.\\ \times \,\left[\cos {\alpha }_{1}(\cos {\alpha }_{1}-\cos {\beta }_{1}\cos i)+\displaystyle \frac{1}{2}{\sin }^{2}i(1-5{\cos }^{2}{\alpha }_{1})\right]\\ +\,{k}_{2}{r}_{2}^{5}{\left(\displaystyle \frac{{\omega }_{2}}{\bar{\omega }}\right)}^{2}\left(1+\displaystyle \frac{{m}_{1}}{{m}_{2}}\right)\\ \times \,\left.\left[\cos {\alpha }_{2}(\cos {\alpha }_{2}-\cos {\beta }_{2}\cos i)+\displaystyle \frac{1}{2}{\sin }^{2}i(1-5{\cos }^{2}{\alpha }_{2})\right]\right\},\end{array}\end{eqnarray} \tag{ 22 }$$

where the (d ω/dt)_E is the GR contribution. If the spin axes of the stars are aligned, then α₁ = α₂ = 0 and β₁ = β₂ = i, where i is the inclination of the orbital plane of the binary.

2.3. Third-body Contribution

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 phase-space 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),

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}=-\displaystyle \frac{G({m}_{1}+{m}_{2})}{{r}^{3}}{\boldsymbol{r}}+{{\boldsymbol{f}}}_{\mathrm{QD},1}+{{\boldsymbol{f}}}_{\mathrm{QD},2}+{{\boldsymbol{f}}}_{\mathrm{rel}}+{{\boldsymbol{f}}}_{3}.\end{eqnarray} \tag{ 23 }$$

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,

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{f}}}_{\mathrm{QD},1} & = & \displaystyle \frac{{R}_{1}^{5}(1+{m}_{2}/{m}_{1}){k}_{21}}{{r}^{4}}\\ & & \times \,\left(\left[5{\left({{\boldsymbol{\Omega }}}_{1}\cdot \hat{{\boldsymbol{r}}}\right)}^{2}-{{\rm{\Omega }}}_{1}^{2}-\displaystyle \frac{6{{Gm}}_{2}}{{r}^{3}}\right]\hat{{\boldsymbol{r}}}\right.\\ & & -\ \left.2({{\boldsymbol{\Omega }}}_{1}\cdot \hat{{\boldsymbol{r}}}){{\boldsymbol{\Omega }}}_{1}\Space{0ex}{3.5ex}{0ex}\right),\end{array}\end{eqnarray} \tag{ 24 }$$

where k₂₁ is the apsidal motion constant of body 1, Ω₁ is the ratio of the rate of stellar rotation and the orbital rotation, and $\hat{{\boldsymbol{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

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{f}}}_{\mathrm{rel}}=-\displaystyle \frac{{{Gm}}_{12}}{{r}^{2}{c}^{2}}\\ \quad \times \,\left(\left[(1+3\eta )\dot{{\boldsymbol{r}}}\cdot \dot{{\boldsymbol{r}}}-2(2+\eta )\displaystyle \frac{{{Gm}}_{12}}{r}-\displaystyle \frac{3}{2}\eta {\dot{r}}^{2}\right]\hat{{\boldsymbol{r}}}\right.\\ \quad -\ \left.2(2-\eta )\dot{r}\dot{{\boldsymbol{r}}}\Space{0ex}{3.5ex}{0ex}\right),\end{array}\end{eqnarray} \tag{ 25 }$$

where m₁₂ = m₁ + m₂, c is the speed of light, and $\eta ={m}_{1}{m}_{2}/{\left({m}_{12}\right)}^{2}$. The force due to a third body (where applicable) is given by

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{3}={{Gm}}_{3}({{\boldsymbol{\beta }}}_{23}-{{\boldsymbol{\beta }}}_{13}),\end{eqnarray} \tag{ 26 }$$

where ${{\boldsymbol{\beta }}}_{{ij}}={{\boldsymbol{r}}}_{{ij}}/{{\boldsymbol{r}}}_{{ij}}^{3}$ is a ratio of the relative positions of the bodies in the system.

To find the positions of all the bodies at any given time, the equations of motion are integrated using a symplectic twelfth-order 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).

The analytic formula discussed in Section 2 gives the apsidal rate, namely $\dot{\omega }$. For a model that can produce a time series for ω(t), the rate of change of ω (e.g., $\dot{\omega }$) 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, ω₀ = 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 $\dot{\omega }=0$), the oscillations must vanish. At each time step, six phase-space coordinates give a unique set of orbital parameters (including ω).

Zoom In Zoom Out Reset image size

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 $\dot{\omega }$ 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 $\dot{\omega }$ 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.

Table 2. Representative Apsidal Rates of Advance and Apsidal Periods as a Result of Their Physical and Orbital Parameters Used in Our Modeling

System	${\dot{\omega }}_{\mathrm{GR}}$ (deg cycle−1)	UGR (yr)	${\dot{\omega }}_{\mathrm{Tidal}}$ (deg cycle−1)	UTidal (yr)
AG Per	0.00161	1240.8	0.03714	53.8
DI Her	0.00068	15358.6	0.00069	14999.0
EM Car	0.00300	1120.9	0.17029	19.8
IQ Per	0.00113	1515.2	0.03220	53.4
QX Car	0.00147	2987.2	0.01346	327.8
V1647 Sgr	0.00076	4233.2	0.01168	276.9
V526 Sgr	0.00091	2073.1	0.01277	148.3
Y Cyg	0.00293	1008.1	0.05303	55.7
CW Cep	0.00238	1128.8	0.07803	34.5
DR Vul	0.00277	800.0	0.08160	27.2
GG Lup	0.00130	1397.5	0.02499	72.9
U Oph	0.00179	926.4	0.15210	10.9
V478 Cyg	0.00276	1027.6	0.20082	14.1
V760 Sco	0.00170	1000.6	0.07117	23.9
ζ Phe	0.00134	1223.9	0.06609	24.9

Note. For all systems except for DI Her, the rate of advance due to GR is much slower than the tidal rate.

Download table as:  ASCIITypeset image

4.1. 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 near-ambiguity in the value of the argument of periastron ω in a nearly circular orbit.

Zoom In Zoom Out Reset image size

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⁻⁴ 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.

4.2. 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%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\rm{\Delta }}\varphi =\left(\displaystyle \frac{6\pi {MG}}{L}\right)(\mathrm{radians}\,{\mathrm{revolution}}^{-1}),\end{eqnarray} \tag{ 27 }$$

where L = (1 − e²)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 $\dot{\omega }$. Sterne (1939) uses the n = 2 (second-order) term in the tidal potential, which is proportional to (r²/a³). 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, $\dot{\omega }$ scales with (r/a)⁷. For the binaries considered here, (r/a)⁷ is only a few percent of (r/a)⁵ 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 spin–orbit 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 $\rlap{/}{\approx }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 ($\dot{{{\boldsymbol{\Omega }}}_{r}}\propto {\boldsymbol{r}}\times {{\boldsymbol{f}}}_{\mathrm{QD}}$). 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 $\dot{{{\boldsymbol{\Omega }}}_{r}}$ is very small in most practical cases. Furthermore, Shakura (1985) shows that the inclination i and nodal angle Ω_N also change 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.

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.

This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 101008324 (ChETEC-INFRA). Support in funding comes from the State of Hessen within the Research Cluster ELEMENTS (Project ID 500/10.006). This work made use of the EMMY supercomputer, provided by the North German Supercomputing Alliance (Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen—HLRN), hosted at Georg-August-Universität Göttingen.

Software: ELC (Orosz & Hauschildt 2000), (Orosz et al. 2019), Numpy (Harris et al. 2020), Matplotlib (Hunter 2007).