A Photometric, Spectroscopic, and Apsidal Motion Analysis of the F-type Eclipsing Binary BW Aquarii from K2 Campaign 3

BW Aqr (EPIC 205982900) was observed by Kepler from 2014 November 15–2015 January 23 during K2 Campaign 3 in both short cadence (1 minute) and long cadence (29.4 minute) exposures (Howell et al. 2014). We downloaded the extracted, long cadence light curve and the short cadence target pixel file from MAST to use in our analysis. The long cadence light curve was produced by the standard Kepler pipeline, so we used the PyKE code (Still & Barclay 2012) to normalize the Simple Aperture Photometry flux and remove exposures taken during thruster firing as noted in the quality flag. For the short cadence data, we used PyKE to extract the light curve from the target pixel file, subtract the background flux, normalize the light curve, and identify and remove exposures taken during thruster firing. Finally, we converted both the long cadence and short cadence light curves from normalized flux to Kepler magnitudes using K_p = 10.233 (Huber et al. 2016) for the out-of-eclipse magnitude.

We obtained 14 nights of data for BW Aqr using the CHIRON echelle spectrograph (Tokovinin et al. 2013) on the CTIO 1.5 m telescope during 2015 May 19–June 19. Three 800 s exposures were taken in fiber mode and averaged for each night, and thorium-argon lamp spectra were taken for wavelength calibration. The CTIO pipeline from Yale University (Tokovinin et al. 2013) was used for all data reduction. We then transformed each spectrum onto a heliocentric, logarithmic wavelength grid and used the reduced flat field spectra for continuum normalization. The CHIRON spectra cover 4500–8900 Å over 61 orders at an average resolving power of R ∼ 27000. The average signal-to-noise ratio of our spectra is S/N = 75 near the blaze peak of the Hα echelle order.

We used the short cadence K2 C3 light curve for eclipse timing because of its precise, one minute time step. To measure the times of mid-eclipse, we first used a spline interpolation to determine the times of ingress and egress at twenty evenly-spaced depths within each eclipse. We then calculated the time of mid-eclipse from the average times of ingress and egress at all depths. The errors in eclipse time correspond to the standard deviation in the center times from each depth and are on the order of 6–10 s, which is not unprecedented. This level of precision has been reached for many other short and long cadence Kepler targets with deep eclipses and high signal-to-noise photometry (e.g., Gies et al. 2012; Conroy et al. 2014; Orosz 2015). The times of mid-eclipse for the K2 C3 data are listed in Table 1, along with the errors and the eclipse type. The deeper, primary eclipses are type I and the shallower, secondary eclipses are type II.

We then combined the K2 eclipse times with the other eclipse times from the literature (Clausen 1991; Bulut 2009; Volkov & Chochol 2014). Weighted, least-squares fits to the primary and secondary eclipse times result in the following linear ephemerides,

$$\begin{eqnarray*}&&\mathrm{Min}\ {\rm{I}}=\mathrm{BJD}\ 2456990.46690(2)+6.7196988(2)E\end{eqnarray*}$$

$$\begin{eqnarray*}&&\mathrm{Min}\ \mathrm{II}=\mathrm{BJD}\ 2456993.65447(2)+6.7196831(2)E\end{eqnarray*}$$

where E is the integer epoch number. The difference in period for the primary and secondary eclipses, albeit small, is indicative of apsidal motion. The true orbital period, known as the sidereal period (P_s), corresponds to the average of the periods from the primary and secondary eclipses. For BW Aqr, we found P_s = 6.7196909 ± 0.0000002 days.

We used the Eclipsing Light Curve (elc) code by Orosz & Hauschildt (2000) to model the long cadence light curve of BW Aqr. We ran elc's genetic optimizer to fit for several parameters: eccentricity (e), longitude of periastron (ω) of the primary star, epoch of periastron (T), orbital inclination (i), relative radius (R/a) of each component, and the effective temperature ratio (${T}_{\mathrm{eff}2}/{T}_{\mathrm{eff}1}$). We held the orbital period fixed to the sidereal period, set the orbital semi-amplitudes from spectroscopy (Section 4.2), and took all limb darkening coefficients from van Hamme (1993). The model light curve also accounts for the time averaging over the 29.4 minute cadence of the K2 measurements.

Because BW Aqr does not show any variations outside of eclipse, such as ellipsoidal variations or reflection effects, the out-of-eclipse points do not hold any information about the orbital parameters of the system. We therefore used only the observed data points during primary or secondary eclipses, then kept only every third data point in order to reduce computation time. This resulted in 174 points used for fitting and 167 degrees of freedom (ν), so we rescaled the optimizer's output χ² values such that ${\chi }_{\min }^{2}=\nu$. The 1σ errors in each parameter were calculated by fitting a parabola to the projected χ² values and taking ${\chi }^{2}\leqslant {\chi }_{\min }^{2}+1$.

The best-fit orbital parameters are listed in the beginning of Table 2, along with the spectroscopic orbital elements found in Section 4.2. The best fit relative radii are R₁/a = 0.0840 ± 0.0010 and R₂/a = 0.0966 ± 0.0009, and the best-fit temperature ratio is ${T}_{\mathrm{eff}2}/{T}_{\mathrm{eff}1}=0.9930\pm 0.0016$. The surface gravity of each star was estimated in elc to be $\mathrm{log}{g}_{1}=4.069$ and $\mathrm{log}{g}_{2}=3.986$, which are consistent with the findings of Clausen (1991). Figure 1 shows the unbinned, folded K2 C3 long cadence light curve and the best-fit elc model.

Zoom In Zoom Out Reset image size

Table 2.  Orbital and Apsidal Motion Parameters

Parameter	Clausen (1991)	This Work
Ps (days)	6.719695 ± 0.000003	6.7196909 ± 0.0000002
T (BJD-2400000)	2444545.5215 ± 0.0006	57165.3471 ± 0.0012
e	0.17 ± 0.01	0.1758 ± 0.0012
ω (deg)	101.3 ± 0.7	103.04 ± 0.07
i (deg)	88.4 ± 0.1	88.37 ± 0.04
K1 (km s−1)	84.2 ± 0.5a	84.56 ± 0.27
K2 (km s−1)	78.4 ± 0.5a	77.83 ± 0.27
γ (km s−1)	9.9 ± 0.3a	9.03 ± 0.16
${M}_{2}/{M}_{1}$	1.07 ± 0.009	1.087 ± 0.004
$a\sin i$ (R⊙)	21.28 ± 0.14	21.23 ± 0.05
Pa (days)	6.719712 ± 0.000003	6.719714 ± 0.000003
$\dot{\omega }$ (deg cycle−1)	0.00090 ± 0.00010	0.00114 ± 0.00017
U (years)	7400 ± 900	5710 ± 830
$\mathrm{log}{\bar{k}}_{2\mathrm{obs}}$	−2.2 ± 0.1	−2.04 ± 0.08

Note.

^aSpectroscopic elements from Imbert (1987).

Download table as:  ASCIITypeset image

There seems to be an artifact in the residuals during the secondary eclipse that we could not model, even after testing several limb darkening laws and coefficients from van Hamme (1993). In the end, the logarithmic law produced the lowest residuals. We also did not account for possible third light from a tertiary companion, which would dilute the eclipses and cause the inclination to be underestimated. However, the inclination is very close to 90°, so this effect would be small.

The anomalistic period (P_a) of a binary system includes the apsidal motion ($\dot{\omega }$) in the observed motion of the binary and is related to the sidereal period by ${P}_{s}={P}_{a}(1-\dot{\omega }/360)$, where $\dot{\omega }$ is in units of deg cycle⁻¹ (Hilditch 2001). P_a and $\dot{\omega }$ can be determined using eclipse timing, because apsidal motion causes the observed eclipse times to deviate from a linear ephemeris over time (often written as observed minus calculated time, O − C). One would calculate the predicted O − C values for a given pair of P_a and $\dot{\omega }$, then compare these to the observed O − C values in order to test different apsidal motion parameters. For example, Khaliullin & Kozureva (1986) first estimated $\dot{\omega }=0.0013\pm 0.0001$ deg cycle⁻¹ for BW Aqr, while Clausen (1991) found $\dot{\omega }$ = 0.0009 ± 0.0001 deg cycle⁻¹ using the method of Giménez & Garcia-Pelayo (1983). Lacy (1992) published a new method of determining the predicted O − C values by using an iterative least-squares minimization to solve the apsidal motion equations numerically, without relying on the simplifications needed in previous studies. Using this method, Bulut (2009) found $\dot{\omega }$ = 0.00092 ± 0.00015 deg cycle⁻¹.

With all of the literature eclipse times and the new K2 eclipse times in hand, we used the method of Lacy (1992) to calculate the predicted O − C for a grid of P_a and $\dot{\omega }$ values. While the apsidal motion parameters do depend on the eccentricity and are usually solved together, we held e and ω fixed at the epoch of the K2 data because they are better constrained by the light curve through the eclipse durations and separations (Matson et al. 2016). We also fixed the inclination and the time of primary eclipse to the values from the light curve solution. We calculated the χ² for each pair of apsidal motion parameters to create the χ² contour map shown in Figure 2. The best-fit apsidal motion parameters are P_a = 6.719714 ± 0.000003 days and $\dot{\omega }$ = 0.00114 ± 0.00017 deg cycle⁻¹, as listed in Table 2. The corresponding O − C diagram is shown in Figure 3. Our $\dot{\omega }$ is slightly higher than that of Clausen (1991), whose errors are likely underestimated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Despite the eclipses measurements for BW Aqr spanning over 100 years, the apsidal motion is so slow that the observations cover only a few percent of the total apsidal motion cycle. This slow precession and the correlation between P_a and $\dot{\omega }$ result in a family of solutions and large uncertainties in these parameters, as is evident in Figure 2. Nonetheless, P_a and $\dot{\omega }$ can be used with the orbital parameters of an eccentric binary system to determine the system's internal structure constant as part of the apsidal motion test, which we discuss later in Section 5.3.

We used a two-dimensional cross-correlation algorithm (Matson et al. 2016) to measure the radial velocities (V_r) of BW Aqr. The template spectra for the primary and secondary components were created from Bluered model spectra of Bertone et al. (2008) based on the effective temperatures, surface gravities, and rotational velocities from Clausen (1991). Clausen et al. (2010) completed an abundance analysis for BW Aqr and found [Fe/H] = −0.07 ± 0.11, so we interpolated the Bluered models to $\mathrm{log}Z/{Z}_{\odot }=-0.07$ for our analysis. Note that Clausen et al. (2010) used the solar abundances from Grevesse et al. (2007) with Z_⊙ = 0.012, while the Bluered models use the solar abundances from Anders & Grevesse (1989) with Z_⊙ = 0.0189.

The cross-correlation algorithm first determines the best-fit velocity separation of the components; the model for the primary star is combined with models of the secondary shifted by different relative velocities, and the maximum correlation is recorded for each trial velocity. A parabolic fit to the resulting maximum correlations determines the relative velocity of the secondary star with respect to the primary. The template for this separation is then cross-correlated with the observed spectrum to determine the absolute velocities of each component.

This process was repeated for all echelle orders covering 4500–7000 Å, discarding any extreme outliers that were more than 3σ away from the mean or where the primary and secondary components' identities were interchanged. The radial velocities for each night were calculated using a weighted average of the velocities from the remaining echelle orders, and the errors correspond to the standard deviations of these values. The final radial velocities, re-derived with templates formed from the updated atmospheric parameters from Section 4.3, are listed in Table 3.

We used the orbit fitting code rvfit¹ by Iglesias-Marzoa et al. (2015) to determine the orbital parameters of the system. This code uses adaptive simulated annealing to fit for any combination of the orbital period (P), epoch of periastron (T), eccentricity (e), longitude of periastron (ω) of the primary star, systemic velocity (γ), and velocity semi-amplitudes (K₁, K₂). Two of our spectra were taken during an eclipse (at phases ϕ = 0.422, 0.958), so the radial velocities measured at these times are not accurate. Unfortunately, these points correspond to phases which influence the shape of the radial velocity curve and resulting fits for e and ω. These parameters, as well as the orbital period and epoch of periastron, are much better constrained by the light curve model, so we held P, T, e and ω fixed and fit only for γ, K₁, and K₂. The best-fit values are listed in Table 2, along with other derived quantities such as the mass ratio ($q={M}_{2}/{M}_{1}$) and the projected semimajor axis ($a\sin i$). We calculated the errors on the three fitted parameters using the Monte Carlo Markov Chain (MCMC) feature of rvfit, where the error in each parameter corresponds to the standard deviation of the Gaussian fit to each MCMC result. These errors are roughly $\sqrt{N}$ better than the individual errors in radial velocity, as expected. The best-fit model radial velocity curve is shown in Figure 4 and the residuals are listed in Table 3. The systemic velocity is very similar to that of Imbert (1987) and we found no periodicity in the radial velocity residuals, which likely excludes the presence of a nearby tertiary companion.

Zoom In Zoom Out Reset image size

In order to determine the atmospheric parameters for each component of BW Aqr, we used the Doppler tomography algorithm of Bagnuolo et al. (1992) to reconstruct the individual component spectra. We adopted Bluered template spectra of each component with atmospheric parameters from Clausen (1991) as starting estimates in the algorithm. Doppler tomography also requires an input for the flux ratio (${F}_{2}/{F}_{1}$), which determines the strengths of the absorption lines in the reconstructed spectra and affects any resulting fits for the effective temperature (${T}_{\mathrm{eff}}$).

In order to determine the flux ratio of BW Aqr, we created a χ² contour across a grid of flux contributions from each component (F₁, F₂) as follows. The Balmer lines are the absorption lines most sensitive to temperature for F-type stars, so we ran Doppler tomography on the Hα echelle order to create reconstructed spectra of each component for our grid of F₁ and F₂ values. At each grid point, we fit for ${T}_{\mathrm{eff}}$ using mpfit (Markwardt 2009) to compare Bluered models of various ${T}_{\mathrm{eff}}$ to the reconstructed spectra of the Hα order. We then used the best-fit model spectrum to calculate χ² across the order, creating a χ² contour for our grid of flux values. We fit a parabola to the projected χ² curve in each dimension to determine the best flux contributions to be F₁ = 0.36 ± 0.03 and F₂ = 0.54 ± 0.04, corresponding to a flux ratio of F₂/F₁ = 1.45 ± 0.08 at 6563 Å.

We created new reconstructed spectra with this flux ratio to use in determining the final effective temperatures of BW Aqr. Because F₁ + F₂ ≤ 1, we think that there is some extra background flux that was not removed during reduction process. To avoid the issue of unreliable absolute line depths in the resulting temperatures fits, we instead used line ratios. We chose eight pairs of absorption lines with varying dependencies on temperature and well defined continuum levels, which compared an Fe i line to either an Fe ii line or the core of Hα. For each pair of absorption lines, we measured the line depth ratios in the reconstructed spectra and model spectra of various ${T}_{\mathrm{eff}}$, then interpolated between the model line ratios to find the effective temperature corresponding to the observed ratio. We took the mean and standard deviation of the results from all line pairs to calculate the final effective temperatures for BW Aqr, ${T}_{\mathrm{eff}1}=6370\pm 270$ K and ${T}_{\mathrm{eff}2}=6320\pm 220$ K. The temperature ratio is then ${T}_{\mathrm{eff}2}/{T}_{\mathrm{eff}1}=0.992\pm 0.090$, which is very similar to the ratio from the light curve analysis. These effective temperatures correspond to spectral types of about F6 and F7 using the temperatures from Gray (2008). Our results are also consistent with the results from Clausen et al. (2010), though they achieved smaller errors from their abundance analysis. Example reconstructed spectra of BW Aqr for the Hα and Hβ echelle orders are shown in Figure 5, along with the final, best-fit Bluered model templates.

Zoom In Zoom Out Reset image size

Finally, we determined the projected rotational velocity ($V\sin i$) of each component from the individual, reconstructed spectra. We identified 12 strong, well separated metal absorption lines to use in comparing the reconstructed spectra to Bluered model spectra of varying $V\sin i$. For each absorption line, we calculated χ² of the models and fit a parabola to the curve to find the best-fit $V\sin i$. We calculated the 1σ errors from the velocities corresponding to ${\chi }^{2}\leqslant {\chi }_{\min }^{2}+1$. The final $V\sin i$ for each component were then calculated from the weighted averages of the results from all twelve absorption lines, which we found to be ${V}_{1}\sin i=12.6\pm 1.9$ km s⁻¹ and ${V}_{2}\sin i\,=14.6\pm 2.0$ km s⁻¹. Both components are rotating near the projected synchronous velocities of 13.4 and 15.4 km s⁻¹, which is rather common, as Lurie et al. (2017) found that 72% of Kepler eclipsing binaries with orbital periods between 2 and 10 days are rotating at the synchronous velocities.

Combining the results from the spectroscopic and light curve analyses, the absolute parameters for BW Aqr are M₁ = 1.365 ± 0.008 M_☉, M₂ = 1.483 ± 0.009 M_☉, R₁ = 1.782 ± 0.021 R_☉, and R₂ = 2.053 ± 0.020 R_☉. The surface gravities of each component are then $\mathrm{log}{g}_{1}=4.071\pm 0.010$ and $\mathrm{log}{g}_{2}=3.985\pm 0.009$. A summary of our results are listed in Table 4 and have errors in mass and radius of about 0.6% and 1%, respectively. They are also consistent with the results of Clausen (1991).

Table 4.  Astrophysical Parameters

Parameter	Primary	Secondary
Mass (M☉)	1.365 ± 0.008	1.483 ± 0.009
Radius (${R}_{\odot }$)	1.782 ± 0.021	2.053 ± 0.020
${T}_{\mathrm{eff}}$ (K)	6370 ± 270	6320 ± 220
$\mathrm{log}g$ (cgs)	4.071 ± 0.010	3.985 ± 0.009
$V\sin i$ (km s−1)	12.63 ± 1.86	14.56 ± 2.01
$[\mathrm{Fe}/{\rm{H}}]$	−0.07 ± 0.11a
${F}_{2}/{F}_{1}$ (at 6563 Å)	1.45 ± 0.08

Note.

^aFixed from Clausen et al. (2010).

Download table as:  ASCIITypeset image

We compared the observed mass and radius to the predictions of several stellar evolution models: the Yonsei–Yale Y² models of Demarque et al. (2004), the Geneva models of Mowlavi et al. (2012), the Granada models of Claret (2004, 2006), the MESA code of Paxton et al. (2011, 2018), and the Victoria-Regina models of VandenBerg et al. (2006). Nonrotating models with scaled solar abundance were used throughout. Note that each model uses a different solar metallicity prescription, which causes a slight scatter in the zero age main sequence positions of each mass track. For each evolutionary model, we estimated the ages of each component star by interpolating the age of the evolutionary tracks at the observed radii.

For the Yonsei–Yale Y² models,² we used the evolutionary track interpolator provided to create tracks at the observed masses and metallicity, shown in Figure 6. These models use the step-function method to characterize convective core overshooting based on an overshoot parameter Λ_ov (sometimes written as α_ov). The amount of overshooting therefore corresponds to $d={{\rm{\Lambda }}}_{\mathrm{ov}}{H}_{p}$, where H_p is the pressure scale height at the convective boundary (Demarque et al. 2004). In the Y² models, Λ_ov is a function of mass and metallicity and is about 0.13 and 0.20 for the components of BW Aqr. These models predict ages of 2.72 and 2.19 Gyr for the primary and secondary components and a mean age of 2.45 Gyr. The difference in age between each component is about 21% of the mean age and is more easily seen in the Y² isochrones shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the Geneva models, we used the online model interpolator³ to produce the evolutionary tracks shown in Figure 6. The Geneva code also uses the step-function method to increase α_ov with mass, but with a slower increase than the Y² models. For both components of BW Aqr, α_ov = 0.05. These models predict ages of 2.80 and 2.20 Gyr and a mean age of 2.50 Gyr (24% difference in age).

The Granada models of Claret (2004, 2006)⁴ cover a grid of mass and metallicity values with a fixed overshoot parameter, α_ov = 0.20. We interpolated between models to the observed mass and metallicity values, but the red hook is very difficult to interpolate across so only the main sequence parts of the evolutionary tracks are shown in Figure 6. These models predict ages of 2.87 and 2.22 Gyr with a mean age of 2.55 Gyr (26% age difference).

The MESA code⁵ computes models for any given mass and metallicity, which the user can specify in the input files. To characterize convective core overshooting, MESA uses the diffusion method based on the overshoot parameter, f_ov, which the user can also specify. Claret & Torres (2017) created a semi-empirical calibration of f_ov based on well-studied eclipsing binaries, so we estimated f_ov = 0.002 and 0.008 for BW Aqr from their calibration. Additionally, MESA can employ the solar abundance prescriptions from several different sources, so we chose to use the solar abundances from Anders & Grevesse (1989) to match the Bluered models. The MESA models predict ages of 2.31 and 1.89 Gyr and a mean age of 2.10 Gyr (20% age difference).

The Victoria-Regina models⁶ also cover a grid of mass and metallicity values, so we interpolated between models to the observed values of mass and metallicity, using only the main sequence portions of the evolutionary tracks. The Victoria-Regina models use a different implementation of convective overshooting, based on the Roxburgh criterion (Roxburgh 1978, 1989). The free parameter F_ov is a function of mass and metallicity, which corresponds to F_ov ∼ 0.2 and 0.4 for BW Aqr. These models predict ages of 2.50 Gyr and 2.29 Gyr with a mean age of 2.30 Gyr (18% age difference).

As seen in Figure 6, all evolutionary models are able to fit the observed temperatures and radii individually, but none of the models are able to reproduce the observed properties with a single isochrone. There seems to be a difference of about 20% between the ages of the primary and secondary stars, with the secondary star predicted to be younger in all cases. This is most easily seen for the Y² isochrones shown in Figure 7, where the observed slope between components is shallower than the slope of the isochrones.

Absolute stellar parameters can be combined with the apsidal motion parameters to probe the internal structure of a binary system because the rate of apsidal motion depends on the internal mass distributions of the component stars. However, one cannot calculate the contributions of each star's internal structure to the observed apsidal motion individually; one can calculate only a mean observed internal structure constant ($\mathrm{log}{\bar{k}}_{2\mathrm{obs}}$) to compare to stellar evolution theory as part of the apsidal motion test (Claret & Giménez 2010).

The observed apsidal motion rate has a classical contribution (${\dot{\omega }}_{\mathrm{clas}}$) and a relativistic contribution (${\dot{\omega }}_{\mathrm{rel}}$). The relativistic component can be calculated from the orbital elements and component masses as described in Appendix A. We found $\dot{\omega }$_rel = 0.00032(5) for BW Aqr, which constitutes about 27% of the observed apsidal motion rate. Once the total observed and relativistic apsidal motion rates are known, one can calculate the classical apsidal motion rate and the corresponding $\mathrm{log}{\bar{k}}_{2\mathrm{obs}}$ value using the equations in Appendix A. Using this method, we found $\mathrm{log}{\bar{k}}_{2\mathrm{obs}}=-2.02\pm 0.08$ for BW Aqr.

From certain evolutionary models, one can predict log k₂ for each star in the binary system ($\mathrm{log}{k}_{21\mathrm{theo}}$, $\mathrm{log}{k}_{22\mathrm{theo}}$). Because rapidly rotating stars are more centrally condensed and have lower log k₂ than slowly rotating stars, these theoretical log k₂ values must be corrected for rotation (Claret 1999) using the equations in Appendix B. For BW Aqr, the corrections due to rotation are small, ${\rm{\Delta }}\mathrm{log}{k}_{21}=-0.00073\pm 0.00011$ and ${\rm{\Delta }}\mathrm{log}{k}_{22}=-0.00097\pm 0.00015$. Additionally, time-dependent tidal distortions due to nonsynchronous rotation affect the predicted apsidal motion rate. Claret & Willems (2002) calculated the necessary corrections for several binary systems and found ${{\rm{\Delta }}}_{\mathrm{dyn}}=0.00054$ for BW Aqr. One can then calculate an average theoretical internal structure constant ($\mathrm{log}{\bar{k}}_{2\mathrm{theo}}$) using the equations described in Appendix B to compare with the value from observations.

We calculated $\mathrm{log}{\bar{k}}_{2\mathrm{theo}}$ from both the Granada and MESA evolutionary models. The Granada models provide theoretical log k₂ values for a grid of masses, $\mathrm{log}g$, and [Fe/H], so we interpolated between these values to determine $\mathrm{log}{k}_{21\mathrm{theo}}=-2.30$ and $\mathrm{log}{k}_{22\mathrm{theo}}=-2.42$ for each component of BW Aqr. After correcting for rotation and dynamic tides, taking the weighted average and then the logarithm, we found $\mathrm{log}{\bar{k}}_{2\mathrm{theo}}=-2.37$ for the Granada models. For the MESA models, log k₂ is not a direct output, but can be calculated from the density and interior mass profiles. Using the procedure in Section 1 of Cisneros-Parra (1970) and described in Appendix C, we integrated the Radau equation to calculate $\mathrm{log}{k}_{21\mathrm{theo}}\,=-2.27$ and $\mathrm{log}{k}_{22\mathrm{theo}}=-2.36$. We corrected for rotation and dynamic tides, took the weighted average and logarithm, and found $\mathrm{log}{\bar{k}}_{2\mathrm{theo}}=-2.33$ for the MESA models.

Both sets of models predict lower log k₂ than the observed value. However, the first condition in the apsidal motion test is that the models must be able to match the observed surface properties of the stars at the same age, because the internal structure constant is highly dependent on the radius of the stars. Neither the Granada nor MESA models were able to fit the temperatures and radii of both components with a single isochrone, so it is not surprising that the theoretical internal structure constants do not agree with the observed value.