A STRICT TEST OF STELLAR EVOLUTION MODELS: THE ABSOLUTE DIMENSIONS OF THE MASSIVE BENCHMARK ECLIPSING BINARY V578 MON

The available time series photometry of V578 Mon covers nearly 40 yr and more than one full apsidal motion period. A summary of the various light curve epochs, including filters and observing facilities used, is presented in Table 2. Photometry from Heiser (2010) includes multi-band light curves spanning 1967–2006 from the 16 inch telescope at Kitt Peak National Observatory (KPNO) and from the Tennessee State University (TSU)–Vanderbilt 16 inch Automatic Photoelectric Telescope (APT) at Fairborn Observatory. The KPNO Johnson UBV light curves comprise 725 data points spanning 1967–1984 with average uncertainties per data point of 0.004 mag computed by Heiser (2010). The APT Johnson BV light curves span 1994–2006 and consist of 1783 data points with average uncertainties per data point of 0.001 mag for B and 0.002 mag for V (Heiser 2010). Light curves from Hensberge et al. (2000) span 1991–1994 from the 0.5-m Strömgren Automatic Telescope (SAT) at La Silla, with 248 data points in each of the uvby filters and average uncertainty per data point of 0.003 mag (Hensberge et al. 2000). We begin our light curve analysis with the observational errors originally estimated by Heiser (2010) and Hensberge et al. (2000). Table 2 lists these average uncertainties, σ₀, as reported by the original authors. However, from our light curve fits (see below) we found that these uncertainties were in most cases underestimated. Thus, we also report as σ in Table 2 the uncertainties that we ultimately adopted for each light curve.

A new series of high-resolution echelle spectra were secured in 2011 December (36 exposures) and 2012 February (8 exposures) with hermes, the fiber-fed high-resolution spectrograph on the Mercator telescope located at the Observatorio del Roque de los Muchachos, La Palma, Canary Islands. hermes samples the entire optical wavelength range (3800–9000 Å) with a resolution of R = 85,000 (Raskin et al. 2011). The observations listed in Table 3 cover the orbital cycle uniformly. Groups of two concatenated exposures allow us to obtain a robust estimation of random noise as a function of wavelength and a check on cosmic-ray events surviving the detection algorithm in the data reduction. In total, 44 exposures were obtained at 19 epochs, 16 of which are out of eclipse. One series of six exposures starts near the primary mid-eclipse. One series of two concatenated exposures taken around secondary mid-eclipse has a significantly lower exposure level, but another one consisting of four concatenated exposures starting around secondary mid-eclipse is available.

Exposure times close to 2100 s were used for most spectra, but in case of one out-of-eclipse epoch, the exposure time was significantly shorter, 1200 s. The signal-to-noise ratio of the spectra is 50–100 at 4000 Å, rapidly increasing to 120 to 200 at 5000 Å, and remaining close to this level at longer wavelengths. The numbers apply to the sum of two concatenated exposures. The reduction of the spectra has been performed using the heres pipeline software package. The spectra resampled directly in constant-size velocity bins (ln λ), which are very nearly the size of the detector pixels, were used. Normalization to the continuum is done separately.

The hermes spectra outnumber the caspec spectra used by Hensberge et al. (2000), but fall short with regard to signal-to-noise ratio. However, they cover a much larger wavelength region, include epochs in both eclipses, and cover the orbit more homogeneously. In the wavelength region covered by both sets, the reconstruction has better signal-to-noise ratio in the caspec set, but the risk of bias due to phase gaps might be higher with the caspec data. Both data sets were obtained in different parts of the apsidal motion cycle.

In the V578 Mon binary system, the eclipses are partial, which causes degeneracy in the light curve solution for the radii of the components. It was checked whether a spectroscopic light ratio has sufficient precision to reduce the degeneracy. This light ratio might be constrained either by the changing line dilution during the eclipse, and/or by constrained fitting of the reconstructed component spectra by theoretical spectra, simultaneously deriving the light ratio as well as the photospheric parameters (Tamajo et al. 2011). In the latter implementation, the light ratio is assumed identical in all observed spectra, hence eclipse spectra are not used.

With partial eclipses of roughly 0.1 mag depth and less for the secondary eclipse at the epoch of the spectroscopy, the line depth in the composite spectrum is affected at the level of 0.5% of the continuum only when the two components have in their intrinsic spectra a line differing by 7% of the continuum depth. The similarity of the components and the rotational broadening in the spectra imply that no metal line approaches this level. Hence, using the changing line dilution to measure the light ratio precisely is challenging. Exceedingly large signal-to-noise ratios would be required to be able to use single or few lines. Including many lines, i.e., large stretches of spectrum offers the opportunity to reduce the requirements on the signal-to-noise ratio. However, bias in tracing the continuum is expected to put an upper limit on the precision with which the light ratio can be measured in a system with components with similar spectra and substantial rotation.

Therefore, we explored the alternative option of constrained fitting, although it is model-sensitive. Spectral disentangling (Hadrava 1995), further referred to as spd is performed in a spectral range of about 100–150 Å (of the order of 4000 bins) in the wavelength range 3900–5000 Å, centered on prominent lines of He i, He ii, and stronger metal lines. The apsidal motion study (Garcia et al. 2011) permitted us to fix the eccentricity e, the longitude of the periastron, ω, for the epoch of the spectra, and the time of periastron passage. The spd code used is FDBinary⁹ (Ilijic et al. 2004).

spd was applied to selected spectral regions of the hermes spectra, well distributed over the full range of Doppler shifts in the orbit (see orbital phases in Table 3), leads to radial velocity amplitudes K1 and K2 compatible with Hensberge et al. (2000) within better than 1 km s⁻¹. Thus, the spectra are reconstructed using the mean orbital elements (Table 4), now also including regions around H_γ and H_δ (H_β has a broad interstellar band centered on its red wing). For the constrained fitting, optimization was done for hydrogen and helium lines only, and for combinations of them. The reconstructed spectra for both out-of-eclipse and in-eclipse phases are shown in Figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The component spectra for different dilution factors can be obtained from a single disentangling computation, followed by an adequate re-normalization. As a starting point for the photospheric parameters, T_{eff, 1} = 30, 000 K and T_{eff, 2} = 26, 400 K are used based on the extensive study of Hensberge et al. (2000) and Pavlovski & Hensberge (2005). The surface gravities of the components are fixed to log g₁ = 4.133 ± 0.018 and log g₂ = 4.185 ± 0.021 as derived in this paper. This suppresses the degeneracy of line profiles of hot stars in the (temperature, gravity) plane. Calculations for a small grid in log g have shown that the effect of fixing log g might produce deviations of about a few tenths of the percentage in determining the light dilution factors.

Optimization of relative light factors includes a search through a grid of theoretical spectra, using a genetic algorithm. A grid of synthetic spectra was calculated assuming non-LTE line formation. The calculations are based on the so-called hybrid approach of Nieva & Przybilla (2007) in which model atmospheres are calculated in LTE approximation and non-LTE spectral synthesis with detailed statistical balance. Model atmospheres are constructed with atlas9 for solar metallicity, [M/H] = 0, and helium abundance by number density, N_He/N_tot = 0.089 (Castelli et al. 1997). Non-LTE level populations and model spectra were computed with recent versions of detail and surface (Butler & Giddings 1985). Further details on the method, grid, and calculations can be found in Tamajo et al. (2011) and Pavlovski et al. (2009).

Depending on the line(s) included, the primary is found to contribute 68%–72% of the total light, with hydrogen lines supporting the larger fractions. Hydrogen suggests a few percent lower temperature for the primary, compared to the starting values. This is compatible with the tendency seen in Figure 7 of Hensberge et al. (2000), that H and He lines for the primary only marginally agree in effective temperature (taking the minimum χ² at the relevant gravity, a 1000 K difference in temperature estimation occurs).

The inconsistency between different indicators underlines the importance of developing a more consistent atmosphere model for these stars. One way, following Nieva & Przybilla (2012), is to include more ionization equilibria by analyzing the full wavelength range covered by the new spectra. This work-intensive analysis is out of the scope of the present paper, but probably indispensable to better constrain the degeneracy in the determination of the radii. Its success might be limited by the rotational broadening in the spectra. Another point of attention is the need to take into account temperature and gravity variations over the surface, due to the slightly non-spherical shape of the stars. Our work shows that the purely photometrically estimated light factors (Table 5) lie within the broader range of light factors (primary to total light) derived from the hermes spectra, 0.68–0.72. However, one should be mindful that further improvement is needed—the spectroscopic estimates may be biased as different indicators are not yet fully compatible.

We use EB modeling software phoebe (Prša & Zwitter 2005) based on the Wilson–Devinney code (Wilson & Devinney 1971; Wilson 1979) for our light curve analysis. We fit light curves spanning 40 yr, covering one full apsidal motion cycle, in Johnson UBV and Strömgren uvby photometry.

Figures 4–7 are the residuals (data model) for our global best-fit model to the light curves for every light curve epoch and filter in Table 2. Overall, the residuals are small—typically ≈0.005 mag. The residuals are significantly larger for light curve epochs 1970–1984 since error bars on the photometry data points measured using photometric plates are larger. We explore ranges for our light curve parameters as listed in Table 6. Our global best-fit matches observations well—the final light curve parameters Ω₁, Ω₂, i, and (T₂/T₁) are listed in Table 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7. Light Curve Analysis Results and Comparison

Light Curve Parameters	This Work & Garcia	H2000
	et al. (2011)
Primary surface potential, Ω1	4.88 ± 0.03	5.02 ± 0.05
Secondary surface potential, Ω2	4.89 ± 0.04	4.87 ± 0.06
Temperature ratio, (T2/T1)	0.858 ± 0.002	0.88 ± 0.02
Inclination, i (deg)	72.09 ± 0.06	72.58 ± 0.3
Eccentricity, e	$0.07755^{+0.00018}_{-0.00026}$	0.0867 ± 0.0006
Angle of periastron, w (deg)	159.8 ± 0.33	153.3 ± 0.6
Ephemeris, HJD0 (days)	2449360.6250	2449360.6250
Total apsidal motion,	$0.07089^{+0.00021}_{-0.00013}$
$\dot{\omega }_{\rm tot}$ (deg cycle−1)
Light curve filters	Strömgren uvby,	Strömgren uvby
Light curve filters	Johnson UBV
Total light curve points	3489	992

Notes. The uncertainties on light curve parameters Ω₁, Ω₂, i, and T₂/T₁ are determined from confidence intervals in Figure 8. Light curve parameters e, w, and $\dot{\omega }_{\rm tot}$ are taken from Garcia et al. (2011). This work utilizes photometry that span one full apsidal motion period (U = $33.48^{+0.10}_{-0.06}$ yr). In contrast to the Hensberge et al. (2000) analysis, this work incorporates apsidal motion in the light curve model. Finally, the temperature ratio from Hensberge et al. (2000) is measured from spectral disentangling.

Download table as:  ASCIITypeset image

4.2.1. Setup

4.2.2. Fitting Method

Our fitting method is adapted from Y. Gómez Maqueo Chew et al. (2014, in preparation). We determine our best-fit global light curve solution by finding a unique set of light curve parameters Ω₁, Ω₂, T₂/T₁, and i, that correspond to the minimum chi square, $\chi ^{2}_{\rm min}$, in a well mapped grid of parameter space. The chi square is a function of the light curve parameters, χ² = χ²(Ω₁, Ω₂, (T₂/T₁), i). We map parameter space by computing χ² for a grid of >10⁵ unique sets of these light curve parameters. We use our map of parameter space to compute the uncertainties on our light curve parameters using confidence intervals. Plots of Δχ² versus stellar radii R₁, R₂, temperature ratio T₂/T₁, and inclination i with confidence intervals are shown in Figure 8.

Zoom In Zoom Out Reset image size

The step-by-step procedure is as follows.

1. We sample a coarse grid of 10⁴ points defined by a range of potential Ω₁, potential Ω₂, inclination i, and temperature ratio T₂/T₁. The parameter ranges and spacings are given in Table 6.

For each grid point, we fit only for the "light levels" in phoebe which is equivalent to the total light contribution from each star in the photometric bandpass. We avoid using the WD2003 differential corrections (DC) fitting algorithm within phoebe to fit our light curve parameters. The DC algorithm can fall into local minima when fitting for many parameters. We compute the total chi square $\chi ^{2}_{k}$ for each light curve fit as the sum of the chi square $\chi ^{2}_{p}$ at each passband and epoch:

$$\begin{equation} \chi ^{2}_{k}\left(\Omega _{1k},\Omega _{2k},\frac{T_{2k}}{T_1},i_{k}\right)=\sum ^{15}_{p}\frac{\chi ^{2}_{p}}{\sigma _p^{2}}, \end{equation} \tag{ 1 }$$

where index k corresponds to a unique point in parameter space (Ω_1k, Ω_2k, T_2k/T₁, i_k). $\chi ^{2}_{k}$ is the total chi square over all light curves at a unique point k. Index p corresponds to a unique light curve passband epoch as specified in Table 2. The chi square at the specific passband $\chi ^{2}_{p}$ is computed as

$$\begin{equation} \chi ^{2}_{p} = \sum ^{N}_{i}\frac{(f -f_{m})^{2}}{\sigma _i^{2}}, \end{equation} \tag{ 2 }$$

where N = N_d − N_p = 3485 is the number of photometry data points N_d minus the number of parameters of interest N_p over all light curve epochs. Each data point has an error bar σ_i. Each light curve at a specific epoch and filter has a multiplicative factor σ_p which takes into account the systematic error. Multiplicative factor σ_p is used to normalize the χ² such that $\chi ^{2}_{\rm min}=N$ or reduced $\chi ^{2}_{\rm min,red}=1.0$. f is the total flux of the binary at an HJD, and flux f_m is the corresponding model. From our coarse grid, we find the minimum total chi square $\chi ^{2} = \chi ^{2}_{\rm min}$ in parameter space.

2. We adjust the error bars of the individual photometry data points for all light curves to take into account any systematic error. For the minimum $\chi ^{2}_{\rm min}$ solution, the passband σ_p is computed for each separate light curve epoch and filter using the following equation:

$$\begin{equation} \sigma _p = \sqrt{\frac{N}{\chi ^2_{\rm min}}}, \end{equation} \tag{ 3 }$$

where N = 3486 as in step 1, and $\chi ^{2}_{\rm min}$ is the minimum total χ² of the coarse grid. We choose to compute the multiplicative factor σ_p to weight each light curve such that the minimum reduced chi squared $\chi ^{2}_{\rm min,red} = 1.0$ for our global best-fit solution. We then rescale the χ² of all other light curve fits using the passband σ_p:

$$\begin{equation} \chi ^{2}_{k} = \sum ^{15}_{p}\frac{\chi ^{2}}{\sigma _p^{2}}, \end{equation} \tag{ 4 }$$

where χ² is un-scaled and $\chi ^{2}_{k}$ is the scaled chi square at a unique point in parameter space k.

3. We perform steps 1 and 2 for a fine grid of >10⁵ points in parameter space around the location of the minimum $\chi ^2_{\rm min}$. In this way, we carefully map out parameter space at the location of the $\chi ^{2}_{\rm min}$. We use multiple fine grids to precisely find our global best-fit minimum. The average grid spacings are 0.005, 0.005, 0.03, and 0.0005, respectively, for Ω₁, Ω₂, i, and T₂/T₁.

We find that the location of the minimum χ² moves slightly, and we recompute the multiplicative factor σ_p for each light curve to account for this, again making $\chi ^2_{\rm min,red} = 1.0$. Finally, we have a global best-fit solution within a finely sampled parameter space. Our global best-fit solution listed in Table 7 corresponds to the point in parameter space where chi square is scaled by σ_p such that $\chi ^{2}_{\rm min,red} = 1.0$.

4.2.3. A Comparison of Light Curve Models

In order to ensure that our light curve solution is robust and thus our light curve parameters are accurate, we compare our best-fit light curve model described above with several other models. As shown in Table 9, we find little effect on our best-fit light curve parameters from using different light curve models. All other models are not as favorable due to larger χ² or temperatures that do not agree with the analysis of the component spectra of V578 Mon from spectral disentangling of Hensberge et al. (2000).

Table 9. A Comparison of Light Curve Models

Model	Ω1	Ω2	i	(T2/T1)	χ2
			(deg)
Best-fit	4.88 ± 0.03	4.89 ± 0.04	72.09 ± 0.06	0.858 ± 0.002	3489.00
Fitting for LD coefficients	4.92	4.89	72.18	0.835	3299.11
Linear law	4.92	4.89	72.15	0.860	3480.01
Logarithmic law	4.91	4.88	72.14	0.858	3503.96
Fix T1 = 28,500	4.94	4.87	72.17	0.857	3460.16
Fix T1 = 31,500	4.92	4.89	72.15	0.867	3488.01
Light reflection	4.90	4.92	72.20	0.856	3522.57
Third light	4.94	4.87	72.24	0.855	3414.93

Notes. The best-fit model uses the square root limb darkening law, a fixed T₁ = 30,000 K, no light reflection, and no third light.

Download table as:  ASCIITypeset image

Table 10. Third Light

Observatory	Year	Filter	L3/Ltot
APT	2005–2006	Johnson B	0.0441
		Johnson V	0.0218
APT	1999–2000	Johnson B	0.0158
		Johnson V	0.0080
APT	1995–1996	Johnson B	−0.0037
		Johnson V	0.0104
APT	1994–1995	Johnson B	0.0059
		Johnson V	0.0046
SAT	1991–1994	Strömgren u	−0.0116
		Strömgren v	−0.0004
		Strömgren b	0.0013
		Strömgren y	−0.0045
KPNO	1967–1984	Johnson U	0.0163
		Johnson B	0.0467
		Johnson V	−0.0100

Notes. Our best-fit light curve model includes no third light. The small amount of third light varies as a function of epoch.

Download table as:  ASCIITypeset image

Table 11. V578 Mon Absolute Dimensions

Parameter	Primary	Secondary
Orbital period, P (days)	2.4084822
Mass, M (M☉)	14.54 ± 0.08	10.29 ± 0.06
Radius, R (R☉)	5.41 ± 0.04	4.29 ± 0.05
Effective temperature, Teff (K)	30{,}000 ± 500	25{,}750 ± 435
Surface gravity, log g (cm s−2)	4.133 ± 0.018	4.185 ± 0.021
Surface velocity, vrot (km s−1)	123 ± 5	99 ± 3
Luminosity, $\log {\frac{L}{L_{\odot }}}$	4.33 ± 0.03	3.86 ± 0.03
Synchronicity parameter, $F=\frac{w}{w_{{\rm orb}}}$	1.08 ± 0.04	1.10 ± 0.03
Apsidal period, U (yr)	$33.48^{+0.10}_{-0.06}$
Observed Newtonian internal	−1.975 ± 0.017
structure constant, log k2, newt

Download table as:  ASCIITypeset image

For all the tests described below, we start at our best-fit solution, then fit all light curves in phoebe for primary potential Ω₁, secondary potential Ω₂, temperature ratio T₂/T₁, and inclination i. Our global best-fit uses a fixed primary temperature T₁ = 30,000 K, no light reflection, and no third light. Furthermore, our global best-fit uses fixed square root law limb darkening coefficients, which are found to work best for hot (T_eff > 9000 K) stars (Diaz-Cordoves & Gimenez 1992; van Hamme 1993). We discuss the different light curve models in the order in which they appear in our summary in Table 9.

• 1.  
  Fitting for Limb Darkening Coefficients. We test the effect of fitting for square root law limb darkening coefficients, finding a lower chi square due to a larger number of free parameters. We find little effect on Ω₁, Ω₂, or i. However, we do find a much lower T₂ = 25, 049. We reject this light curve model since T₂ = 25, 049 is significantly outside of the acceptable range for T₂ = 26, 400 ± 400 from the spectral disentangling of Hensberge et al. (2000). We therefore perform another test: we keep T₂/T₁ fixed to our best-fit value, and fit for the limb darkening parameters, Ω₁, Ω₂, and i. We again find little effect on Ω₁, Ω₂, or i.
• 2.  
  Using a Different Limb Darkening Law. We test the linear cosine and logarithmic limb darkening laws, finding little effect on our light curve parameters. The linear cosine law has a lower χ² = 3480.01 than our best-fit model χ² = 3489.00. The light curve model with logarithmic limb darkening has a larger χ² = 3503.96-, and we therefore reject this model. See Table 8 for a list of the theoretical limb darkening coefficients for each light curve model that we test.
• 3.  
  Changing the Assumed Primary Star Temperature. We test the effect of changing our adopted primary star effective temperature T₁. Our adopted primary temperature for our best-fit solution is T₁ = 30, 000 ± 500 K. Once again, we find little effect on Ω₁, Ω₂, i, or T₂/T₁.We start with our best-fit global solution, but set T₁ = 31,500 K and T₁ = 28,500 K, 3σ above and below our adopted primary star effective temperature. Fits with lower primary temperature T₁ result in a better χ², however, T₁ < 29,000 K does not agree with the spectral disentangling analysis from Hensberge et al. (2000). This may be due to the fact that the phoebe light curve analysis constrains the temperature ratio and not the individual temperatures themselves. Further light curve tests at lower preferred temperatures T₁ and T₂ confirmed that changing effective temperatures has little effect on the geometric parameters, Ω₁, Ω₂, and i.
• 4.  
  Light Reflection. We fit our light curve model with one light reflection. We find an inclination i larger by 2σ. However, the χ² = 3522.57 is higher than our best-fit χ² = 3489.00. We reject this model on this basis.
• 5.  
  Third Light. We test the possibility of third light and its effect on our best-fit parameters. We fit for a third light parameter starting from our best-fit light curve solution. The third light model has a lower χ² due to a larger number of free parameters. We find Ω₁ and i to be larger by 2σ and 2.5σ, respectively, from our best-fit model.However, the third light parameter L₃ varies on the order of an apsidal period of the system. As shown in Table 10, we find at max a small contribution of third light (L₃/L_tot) ≈ 0.045 for the Johnson B filter of light curve epochs 1967–1984 and 2005–2006. This is likely due to phoebe using the L₃ parameter to minimize the small systematic error of 0.005 mag in the residuals of the 1967–1984 and 2005–2006 light curve epochs. Furthermore, the systemic velocity measured with the hermes spectra and the caspec spectra in Hensberge et al. (2000) does not give any evidence for a large third body in the system that would contribute significantly to the light. This is consistent with the third light tests performed here.

4.2.4. Uncertainties on Light Curve Parameters

4.2.5. Consistency of Light Fractions

As mentioned by Torres et al. (2010), an important consistency check of our light curve solution is that the light fractions l_{f, 1} = (l₁/l₁ + l₂) determined from spectroscopy and photometry agree. Given the small degeneracy between R₁ and R₂ as seen in Figure 8, we compare our photometrically determined light fraction with the light fraction from the HERMES spectral disentangling and a previous combined light curve and spectral disentangling analysis from Hensberge et al. (2000). We find that all three light fractions agree with each other to within 1.2σ. A comparison of light fractions is shown in Table 5.

For each of the ≈10⁵ light curve fits to our 40 yr of photometry data, we compute the light fraction at each of the passbands, Johnson UBV, and Strömgren uvby photometry, l_{f, 1}(λ) = (l₁(λ)/l₁(λ) + l₂(λ)), where l₁(λ) and l₂(λ) are the contribution of the primary and secondary star to the total light at a specific passband out of eclipse. The distribution of light fractions l_{f, 1} for light curve models with confidence intervals of 1σ and 2σ are shown in Figures 9 and 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Hensberge et al. (2000) uses an iterative, combined light curve and spectral disentangling analysis using the Wilson–Devinney light curve modeling program to compute their light curve parameters. We find that R₁ = 5.23 ± 0.06 R_☉ from Hensberge et al. (2000) is 2.5σ discrepant from our best-fit R₁ = 5.41 ± 0.04 R_☉. We find that our inclination i = 72.09 ± 0.06 deg is 1.6σ discrepant from i = 72.58 ± 0.30 deg from Hensberge et al. (2000). These discrepancies are likely due to the addition of apsidal motion and an updated eccentricity determined in Garcia et al. (2011). Apsidal motion and eccentricity can affect the potentials Ω₁ and Ω₂, and hence the determination of the radii at a low level. The potential Ω for a non-circular orbit is a function of eccentricity (see Wilson 1979). The addition of more light curve epochs may also play a role. Hensberge et al. (2000) only use the 1991–1994 light curve epoch with Strömgren uvby photometry. As a check, we also recover the Hensberge et al. (2000) light curve solution when we fit only the 1991–1994 light curve epoch. Finally, simply the addition of more photometry data points may play a role. We use 3489 photometry data points in our light curve solution, whereas Hensberge et al. (2000) use 992. Our best-fit secondary radius R₂ = 4.29 ± 0.05 R_☉ is in agreement with 4.32 ± 0.07 R_☉ from Hensberge et al. (2000). Our best-fit temperature ratio (T₂/T₁) = 0.858 ± 0.002 is in agreement with the temperature ratio of 0.88 ± 0.020 from an analysis of the disentangled component spectra (Hensberge et al. 2000).

In Figure 11, we place the primary and secondary star on mass–radius and log g − log T_eff isochrones for each set of models. For the stellar evolution models to pass the isochrone test, the models should predict a common age for both components of V578 Mon within the uncertainty. Given how different the masses of the primary and secondary stars of V578 Mon are, the isochrone test provides a stringent test of stellar evolution models. We also match all evolution models to the rotational velocities of the primary and secondary star.

Zoom In Zoom Out Reset image size

We find several Geneva13, Utrecht11, and Granada04 models predict masses, radii, and temperatures for the components of V578 Mon that fall within the 1σ uncertainty of the measured absolute dimensions. Therefore, we estimate an age range for each star as shown in Table 12. The age difference for the Geneva13, Utrecht11, and Granada04 models is given as the smallest possible difference between the ages of the two stars given the age range of each star.

For the Geneva models, we use isochrones with initial rotational velocities of (v_i/v_crit) = 0.30 and (v_i/v_crit) = 0.35, which allow us to match the observed rotational velocities for each star. We interpolate the model evolution tracks for the primary and secondary star using the online interactive tool provided by the Geneva group.¹¹ Attempts to match the observed rotational velocities of V578 Mon with lower ((v_i/v_crit) < 0.30) or higher ((v_i/v_crit) > 0.40) initial velocities for either star were unsuccessful. Attempts to find a single initial rotational velocity to reproduce the current observed rotational velocities for both stars with reasonable predicted radii and masses were also unsuccessful. However, given that V578 Mon is very near synchronization with the orbital period (F₁ = 1.08 ± 0.04, F₂ = 1.10 ± 0.03), the rotational history of V578 Mon could be different from the best matched v_i/v_crit found here. If the initial velocities of the components of V578 Mon were larger at the ZAMS than the orbital velocity, the stars could spin down to synchronize with the orbital velocity. Conversely, if v_i/v_crit was smaller than the orbital velocities, then the components of V578 Mon could spin up (Song et al. 2013). From Table 13, we find an age difference of 1.6 Myr for mass–radius isochrones and an age difference of only 0.1 Myr for log g − log T_eff isochrones. It is easier to find consistency for the latter isochrones given our uncertainty in the effective temperatures of the two stars. We find that a primary radius of R₁ = 5.50 R_☉ and a secondary star radius of R₂ = 5.20 R_☉ yields common ages for the Geneva13 models. However, these radii are 3σ larger and 3σ smaller than our best-fit model, respectively.

For the Utrecht11 models, we use isochrones that match the observed surface velocities of the components of V578 Mon, v_{1, rot} = 123 ± 5 km s⁻¹ and v_{2, rot} = 99 ± 3 km s⁻¹. The Utrecht11 models are computed at very small steps of mass and initial rotational velocity, such that interpolating between model tracks is unnecessary. From Table 13, we see a marginally common age (age difference 0.4 Myr) for mass–radius isochrones, and a common age of 3.5 ± 1.5 Myr for log g − log T_eff isochrones. The models were computed at solar metallicity by Dr. I. Brott (2014, private communication).

We compute the Granada04 models at the masses of the primary and secondary stars and chose rotational velocities to match the observed rotational velocities of V578 Mon. We attempt to match the absolute dimensions of V578 Mon to log g − log T_eff or, alternatively, mass–radius isochrones for V578 Mon. We find an age gap of 1.5 Myr for mass–radius isochrones and a marginally common age for log g − log T_eff isochrones when both stars have an overshoot of α_ov = 0.2. Again, finding a match on the log g − log T_eff isochrones is easier given the greater uncertainty in the effective temperatures.

In an attempt to match the ages of the two stars on a mass–radius isochrone, we also compute Granada04 models for α_ov = 0.4 and α_ov = 0.6. Figure 12 demonstrates the time evolution of the radii for V578 Mon for these different models. We find a near match on a single mass–radius isochrone with an age difference of only 0.2 Myr if we assume that the primary star has a convective overshoot of α_ov = 0.6 and the secondary star has a convective overshoot of α_ov = 0.2. We also find a common age of 5.5 ± 1.0 Myr for the log g − log T_eff isochrone. This does not mean that an α_ov = 0.6 for the primary star is correct for V578 Mon, merely that a higher convective overshoot allows for compatible ages between the two stars. High convective overshoot has been found to work in matching other EBs on a single isochrone (Claret 2007).

Zoom In Zoom Out Reset image size

In general, we find younger ages by ≈1 Myr for the Utrecht11 models of V578 Mon and similar ages for the Geneva13 and Granada04 models. This can be attributed to the larger convective overshoot of α_ov = 0.355 included in Utrecht11 models compared with Geneva13 models (α_ov = 0.2). While the primary star for the Granada04 models does have an even higher convective overshoot of α_ov = 0.6, the models do not include rotational mixing, which also extends the MS lifetime of the stars.

Measurement of apsidal motion in eccentric binary systems allows for a stringent test of the internal structure constant k_{2, theo} predicted from stellar evolution models (e.g., Claret & Giménez 2010). It is not possible to separate out each individual star's contribution to the apsidal period U from Newtonian apsidal motion.

The apsidal motion for V578 Mon was measured by Garcia et al. (2011). The observed apsidal motion of V578 Mon, $\dot{\omega }_{\rm tot} = 0.07089^{+0.00021}_{-0.00013}$ deg cycle⁻¹, has contributions from both Newtonian and general relativity components (Claret & Giménez 2010):

$$\begin{equation} \dot{\omega }_{\rm obs} = \dot{\omega }_{\rm newt} + \dot{\omega }_{\rm GR}, \end{equation} \tag{ 6 }$$

where $\dot{\omega }_{\rm GR}$ is given by

$$\begin{equation} \dot{\omega }_{\rm GR} = 0.002286\frac{M_1+M_2}{a(1-e^{2})}. \end{equation} \tag{ 7 }$$

We find that $\dot{\omega }_{\rm GR}=$ 0.002589 ± 0.000015 which is only 4% of the Newtonian apsidal motion $\dot{\omega }_{\rm newt}=$ 0.06830 ± 0.00017.

Both the Newtonian and general relativistic observed apsidal motions $\dot{\omega }_{\rm newt}$ and $\dot{\omega }_{\rm GR}$ have associated observed internal structure constants k_{2, obs}. The internal structure constant is twice the tidal Love number (Kramm et al. 2011), and is related to the density profiles, degree of sphericity, orbital parameters, masses, and rotation rates of both components of a binary star. Specifically, the internal structure constant is related to the solution of the Radau differential equation as in Equation (3) of Claret & Giménez (2010). Importantly, constant k_{2, obs} is one the few ways to directly constrain the internal structure of stars.

From the precise observed apsidal motion, we compute the observed internal structure constant, k_{2, obs} = k_{2, obs}(M₁, M₂, R₁, R₂, P, U, F₁, F₂, e), where U is the apsidal period, given by the equations (adopted from Claret & Giménez 2010)

$$\begin{equation} k_{\rm 2,obs} = \frac{1}{c_{21} + c_{22}}\frac{P}{U} \end{equation} \tag{ 8 }$$

$$\begin{equation} c_{2i} = \left[(F_i)^{2}\left(1+\frac{M_{3-i}}{M_i}\right)f(e)+15\frac{M_{3-i}}{M_i}g(e)\right]\left(\frac{R_i}{a}\right)^{5} \end{equation} \tag{ 9 }$$

$$\begin{equation} f(e)=(1-e^{2})^{-2} \end{equation} \tag{ 10 }$$

$$\begin{equation} g(e) = \frac{(8+12e^{2}+e^{4})f(e)^{2.5}}{8}. \end{equation} \tag{ 11 }$$

We compute the internal structure constant due to the Newtonian apsidal motion, log k_{2, newt} = −1.975 ± 0.017, and due to general relativity, log k_{2, GR} = −3.412 ± 0.018. The Newtonian apsidal motion is much larger than the general relativistic component, and therefore the internal structure constant is also much larger.

We compute the theoretical internal structure constant, k_{2, theo} using the methods of Claret & Giménez (2010). The theoretical k₂ constant was corrected for by rotation (Claret 1999) and dynamical tides (Willems & Claret 2002). The theoretical internal structure constant is a combination of the internal structure constants for both stars, such that

$$\begin{equation} k_{\rm 2,theo} = \frac{c_{21}k_{21}+c_{22}k_{22}}{c_{21}+c_{22}} \end{equation} \tag{ 12 }$$

which can then be compared to observations.

We find the predicted Newtonian apsidal motion to be $\dot{\omega }_{\rm theo}=$ 0.06883 ± 0.00017 and consequently the predicted Newtonian internal structure constant to be log k_{2, theo} = −2.005 ± 0.025. This is in very good agreement with the observed log k_{2, obs} = −1.975 ± 0.017. From Equation (9), the parameter c₁₂ is about 67% larger than c₂₂. Therefore, the weighted contribution of the primary dominates the theoretical apsidal motion. V578 Mon is a relatively young system; therefore, log k_{2, theo} is almost constant during the early phases of stellar evolution. The apsidal motion test is therefore complementary to the isochrone test. Claret & Giménez (2010) compile a list of eclipsing binaries with apsidal motion, demonstrating good agreement between observed and predicted apsidal motions. V578 Mon continues this trend of agreement between theoretical and observational internal structure constants. For this relatively young system, matching the radii, temperatures, and masses of isochrones is key given that we have so few young massive EBs with non-equal mass ratios.