A peanut-shaped W UMa contact binary is a very interesting system where two stars contact each other physically and a common gas envelope surrounds them (Eggleton 2012). Contact binaries consist of F~K late-type stars with orbital periods of 0.2~1 day and mass ratio (q) ranging from 0.06~1.0. They have light curves with the depths of the primary and secondary eclipses being almost equal or differing slightly, indicating nearly equal temperatures of both components. In some cases, the light level of phase 0.25 (Max I) differs from that of phase 0.75 (Max II), the so called O’Connell effect and the difference varies over time. These phenomena are usually explained by cool or hot spots on either or both components and the time variation of the spots. W UMa contact binaries are classified into A-subtype and W-subtype. When the more massive star is eclipsed by the less massive one during the primary eclipse, it is called A-subtype. On the other hand, when the more massive star obscures the less massive one at the primary eclipse, it is called W-subtype (Binnendijk 1970, 1977). At present, there are two conflicting arguments that W-subtypes evolve into A-subtypes (Flannery 1976; Lucy 1976; Wilson 1978; Mochnacki 1981; Hilditch 1989; Awadalla & Hanna 2005; Eker et al. 2006; Gazeas & Niarchos 2006) and that there is no evolutionary link between them (Maceroni & Van’t Veer F 1996; Yildiz & Doğan 2013).

Although there is no general agreement among scientists on the origin, structure, and evolution of W UMa contact binaries, an implicitly accepted evolution scenario is that late-type detached binaries has started evolving, via the evolutionary stage of semi-detached binaries like Algol, into current contact binaries and then their components will finally merge into a single star (Huang 1966; Van’t Veer 1979; Vilhu & Rahunen 1979; Mochnacki 1985; Webbink 1985; Guinan & Bradstreet 1988; Eggleton & Kiseleva-Eggleton 2001; Yakut & Eggleton 2005; Stepien 2007; Eggleton 2012). Further discussion on this topic will be addressed in the final section of this paper.

In accordance with the brief investigation described above, W UMa contact binaries are considered to be important celestial bodies for the understanding of formation and evolution of close binaries. Thus, it is necessary to make various observations on W UMa binaries to obtain detailed characteristics of the binaries. Although many new W UMa binaries have been discovered through all sky surveys with various purposes, the number of contact binaries whose photometric solutions are well known is less than 400 (Pribulla et al. 2003; Csizmadia et al. 2007) and the number of contact binaries whose absolute physical dimensions are accurately determined through both photometric and spectroscopic observations is only around 100 (Yildiz & Doğan 2013). Therefore, it is essential to generate accurate solutions of contact binaries that are as spectrophotometrically accurate as possible.

Recently, we have been making photometric observations of close binary stars whose photometric characteristics are not well known (Kim & Jeong 2012; Kim et al. 2014). As one task of the project, we have observed NSVS 1461538 (2MASS J23102148+5859170, SWASP J231021.48+585917.2) and its neighboring PV Cas at the same time. Hoffman et al. (2009) first found out that NSVS 1461538 is a variable. They classified NSVS 1461538 as an Algol/beta Lyr-type eclipsing variable that has an orbital period of 0.3912 day with a light variability of 0.636 mag. However, our multi-band light curves of NSVS 1461538 indicate that this star is a typical W UMa contact binary, unlike the previous classification. In this paper, the light curves of NSVS 1461538 and the orbital period are analyzed to determine photometric properties, absolute physical quantities, and the evolutionary status of the binary for the first time.

We have performed photometric observations of NSVS 1461538 for 45 days from 2009 to 2011 using a 61-cm reflecting telescope located at Mt. Sobaek and a CCD camera with a Johnson BVRI filter. A FLI 2K CCD for the observing seasons of 2009-2010 and a PIXIS 2K CCD for 2011 were used and the field of views (FOV) for the two CCD systems are 20'×20' and 17.6'×17.6', respectively. Both were cooled by liquid nitrogen. All the raw CCD images obtained are pre-processed by compensating for bias, dark, and flat using the IRAF/CCDPRO package and post-processed using IRAF/DAOPHOT. Further details of raw data processing are described in Kim et al. (2014).

A total of 9,781 (B:2,286, V:2,494, R:2,519, I:2,482) photometric measurements were obtained through the three year-observations. Two stars within the FOVs, which have stellar magnitudes and color indices similar to those of NSVS 1461538, were chosen as a comparison star (GSC 04010-01525) and a check star (TYC 4010-1432-1) in order to generate differential light curves of NSVS 1461538. The coordinates and the magnitudes of the variable star, the comparison star, and the check star are listed in Table 1. As shown in Fig. 1, NSVS 1461538 (V), the comparison (C), and the check (K) stars are indicated as ‘V’, ‘C’, and ‘K’, respectively, together with PV Cas, within the FOVs of our two CCD camera systems. All the photometric data of NSVS 1461538 relative to the comparison star are archived and available at our local website (http://binary.cbnu.ac.kr/index.php?mid=pot_data).

The differential light curves of the variable and the check star relative to the comparison star are plotted in the top and the bottom panels of Fig. 2, respectively. The standard deviations of the filtered light curves of the check star are 0^m.020, 0^m.017, 0^m.018, and 0^m.020 in B, V, R, and I, respectively. Thus, the mean standard deviation of our observations is about 0^m.02. The light curves in Fig. 2 were drawn with the light element of Eq. (3) to be discussed in the next section. The light curves of NSVS 1461538 in the upper panel of the figure show that NSVS 1461538 belongs to a typical W UMa contact binary star rather than an Algol/beta Lyr-type binary star classified by Hoffman et al. (2009). Moreover, a total eclipse lasting for about 51 minutes at primary eclipse indicates that NSVS 1461538 belongs to the group of W-subtype. In addition, the O’Connell effect that the luminosity at Max I is brighter than that of Max II is evident for all the BVRI filters. The depths of eclipses for the B filter are 0^m.64 and 0^m.54 for the primary eclipse and the secondary eclipse, respectively, and the difference turns out to be 0^m.10. Also, for the I filter, the depth at the primary minimum is 0^m.55, showing a difference of 0^m.09 from that of the secondary eclipse, 0^m.46. To sum up, as wavelengths become longer, the eclipse depths become shallower; however, the differences in depth between the primary eclipse and the secondary eclipse are almost independent of wavelengths. This indicates that both the primary star and the secondary star of NSVS 1461538 have spectral types of later-type and the temperature difference is small. The differential magnitudes of the light curves at the characteristic phases for NSVS 1461538 are listed in Table 2.

In order to investigate the variation in the orbital period of NSVS 1461538, we first calculated times of minimum light from our filtered observations using the Kwee & van Woerden (1956) method. A total of 20 weighted mean timings were newly obtained. In addition, the same method was applied to super wide angle search for planets (SWASP) data to obtain 15 new times of minimum light and all timings are listed in Table 3. A least squares method was applied to these data to determine the linear light element:

Fig. 3 shows the eclipse timing diagram of NSVS 1461538 plotted using the light elements of Eq. (1). In this figure, filled circles and open circles represent the primary eclipses and the secondary eclipses, respectively. From the figure, we can deduce that the orbital period of NSVS 1461538 could have been undergoing a cyclical change even for a short time period of 4 yr. Thus, we have tried to fit the (O-C₁) residuals in Fig. 3 to a complex ephemeris consisting of linear light elements plus a sine function as follows:

In order to determine 5 unknowns (ΔT₀, ΔP, a, b, c) in Eq. (2), the Levenberg-Marquardt method (LM) was applied and the results are listed in Table 4. The solid curve in Fig. 3 represents the theoretical curve calculated with Eq. (2) and the solution parameters in Table 4. The period of the sine function (P_sine = 2πP/365.24b) was calculated to be about 5.6 yr and the standard deviation turned out to be a relatively large value of ±3.5 yr. This appears to be caused by the relatively large scatter of the timings, about ±0.^d002, as shown in Fig. 3. It is expected that further observations would help improve and verify the period of 5.6 yr. Also, the new light element is obtained as

Eq. (3) was used to generate the light curves shown in Section 2.

BVRI light curves of NSVS 1461538 were analyzed based on the binary star model of Wilson & Devinney (1971, hereafter WD). We first tried to obtain the solution by feeding all the data (9,781 measurements) into the WD code at once but found that it took too long for a single calculation. Hence, as shown in Fig. 4, we have generated normalized BVRI light curves consisting of normal points, which are averaged over a phase interval of 0.005 from the BVRI light curves, and analyzed those light curves with Mode 3 (contact mode) of WD.

In order to obtain a unique solution using the WD binary star model, the initial parameters of the binary system should be known. Among those, the most important are the temperature and the mass ratio of two component stars. However, since those data for NSVS 1461538 are not available, it is necessary to assume these parameters in various ways. As the light curves presented in Figs. 2 and 4 clearly show that the binary system belongs to W-subtype the primary star, which has smaller mass but higher surface temperature, is obscured at the primary eclipse. On the contrary, the secondary star has lower temperature but larger mass and luminosity than the primary star. Thus, we should estimate the initial temperature of the secondary star rather than that of the primary star, which has lower mass. For this purpose, we have collected color indices of NSVS 1461538 from 6 catalogues listed in Table 5 and produced temperatures that correspond to color indices, as listed in Table 5. As can be seen in Table 5, the temperature ranges widely from 4,984 K to 6,366 K. We adopted the average value of 5,340 K as the temperature of the secondary star. The spectral type of the secondary star corresponding to the temperature is G8. Since the component stars of W UMa contact binaries share a common envelope as their outer atmospheres, the temperature difference between the components is small. Hence, since both the primary and the secondary stars were assumed to have a convective gas, gravity darkening exponents (g_h=g_c) and albedos (A_h=A_c) were fixed to theoretical reference values of convective gas, which are 0.32 and 0.5, respectively (Lucy 1967).

In addition, the initial value of the mass ratio was estimated with the so called q-search method. For a series of q the method searches for the global minimum value of the squares of the difference between the observed light curves and the theoretical light curves (Σ(O-C)², abbreviated to Σ). The mass ratio with the global minimum of Σ is chosen as the estimated initial value. Further information about this method was fully described by Jeong & Kim (2013). Fig. 5 indicates that our q-search reaches the minimum at q=3.39. For the other initial conditions, the potential was determined using the relation between the mass ratio and the potential (Hilditch 2001) and the initial value of inclination was deduced to be 88°, because the total eclipse at primary eclipse indicates a high inclination close to 90°.

In Mode 3 calculation, we have adjusted 6 parameters: inclination (i), mass ratio (q), temperature of the primary star (T_h), potential (Ω_h=Ω_c=Ω), relative luminosity of the primary star (l_h), and relative luminosity of the tertiary body (l₃). Among these parameters, the inclination and the mass ratio are strongly correlated and thus weakens the convergence of the solution. Thus, several subsets of parameters that have relatively low correlations between them were used for our calculation (Wilson & Biermann 1976). We repeated the process until the probable errors of the parameters become less than their correction values. The second and the third columns of Table 6 list our results and the dotted lines in Fig. 4 represent the theoretical light curves drawn with our solution. As shown in the figure, the dotted light curves do not meet the observed slight curves between phase intervals of 0.5~1.0, showing the O’Connell effect. In order to explain the O’Connell effect, we have tried to adjust light curves by adding a cool spot on the primary or the secondary star. As a result, the calculated Σ values are almost the same for the two cases. Generally, it is known that the more massive star in W UMa contact binaries is more active magnetically (Mullan 1975; Eaton 1986). Hence, we prefer the solution of the cool spot on the secondary star. Final photometric solutions obtained are listed in the fourth through the seventh columns of Table. 6. The solid lines in Fig. 4 represent the theoretical light curves of the solutions with the cool spot on the more massive secondary star. In addition, the residuals of the observed light curves from the theoretical curves are shown in Fig. 6. The 3-dimensional Roche geometries of NSVS 1461538 at four characteristic phases are drawn in Fig. 7 with the cool spot model on the less-massive primary and the more-massive secondary stars. It is quite interesting that the luminosity of a tertiary body is detected as 5.3 %, 2.2 %, 0.5 %, and 0.1 %, in sequence of BVRI filters, respectively. This corroborates that the variation of the orbital period, 5.6 yr, detected in the previous section might be caused by the light-time effect due to the suggested tertiary body. If there is a tertiary body, the luminosity obtained for each filter suggests that the tertiary body is an early-type star rather than a late-type star.

We have performed the first photometric study of NSVS 1461538, classified as an Algol/beta Lyr-type binary system in the NSVS sky survey. In summary, we have found that NSVS 1461538 is a new W UMa contact binary whose orbital period is 0^d.3913025, and that, based on an analysis of times of minimum light, the orbital period of the binary system shows a cyclic change of about 5.6 yr (±3.5 yr); this period has a significant uncertainty and the periodicity should be verified and improved by further observations, however.

The photometric solutions of NSVS 1461538 reveal that the temperature ratio of the secondary component star to the primary component star is 0.93 with a large difference of 421 K. On the contrary, the mass ratio, radius ratio (volume radius), and the luminosity ratio through the V filter turned out to be 3.46, 1.73, and 2.06, respectively. These results indicate that the primary star has far less mass, smaller radius, and is two times dimmer than the secondary star. The degree of contact of two stars, f, is about 30 % and both component stars are in deep moderate contact with the inner Roche lobe. The inclination of about 89° is so large that there occurs a total eclipse of 51 minutes during the primary eclipse. In addition, the luminosity of a tertiary body was found in all filters. The O’Connell effect shown in the light curves can be explained by a cool spot; however, it is not possible to determine which component has the spot. As any radial velocity curves of NSVS 1461538 are unavailable, it was not possible to obtain accurate absolute dimensions of this binary system. However, for further discussion on the cause of the period change, the evolution of the binary system, and so on, we have estimated the absolute physical quantities of NSVS 1461538 indirectly, as listed in Table 7. The method to estimate these quantities is described in detail in Kim et al. (2014).

As discussed in Section 3, the period variation of 5.6 yr can be explained by a magnetic model of Applegate (1992) or the light-time effect due to a tertiary body. Under the assumption that the variation is caused by the former reason, the Applegate parameters were calculated as listed in Table 8. In the calculation the mass of the convective shell is assumed to be one-tenth of the masses of the primary star and the secondary star, respectively. As can be seen in the table, the Applegate model for both stars show luminosity variation, ΔL/L > ~2, which is too large to be reasonable, and hence this model cannot explain the 5.6 yr cyclic variation of period. Thus, there is a good possibility that this variation is caused by the light-time effect. In this case, the minimum mass of the tertiary body is calculated to be 0.71 M_⊙. Also, as shown in Table 6, since the tertiary body seems to emit more energy in shorter wavelengths, if there exists a tertiary body, it may be an early-type star.

Finally, we have examined the mass-radius (M-R, Fig. 8) and mass-luminosity (M-L, Fig. 9) diagrams, which show the evolutionary status of NSVS 1461538 on a log-scale. The solid and dotted lines in each figure represent zero age main-sequence (ZAMS) and terminal age main-sequence (TAMS) lines, respectively, and both are taken from the PADOVA star evolution model against solar abundance. The W UMa type binaries whose absolute dimensions are well determined from both photometry and spectroscopy are plotted together in each figure. The data are based on 100 contact binaries listed in Yildiz & Doğan (2013). In these figures, the primary stars are more massive than the secondary stars. The cool massive star and hot less-massive star of W-subtype are displayed as an open circle and a plus sign, respectively. And the hot massive star and cool less-massive star of A-subtype are displayed as an open square and a cross sign, respectively. Similarly, for NSVS 1461538, the primary star (hotter but less massive) and the secondary star (cooler more massive) are represented by an open star and a closed star symbols, respectively.

As shown in Figs. 8 and 9, the primary star and the secondary star of NSVS 1461538 are consistent with the distribution of other W UMa component stars. In the M-R diagram, while most of the more massive primary stars are located between the ZAMS and the TAMS, less massive secondary stars are located above the TAMS. This indicates that less massive secondary stars are excessively large for their mass if we assume that primary and secondary stars are main sequence stars. In the M-L diagram, while most of the primary stars are located between the ZAMS and the TAMS or below ZAMS, secondary stars are located above the TAMS. This indicates that less massive secondary stars are excessively bright for their mass and a good number of primary stars are rather dimmer compared to the luminosity corresponding to their mass. There are two different explanations and views for these observations. The first is as follows: since contact binaries are not in thermal equilibrium, they undergo a periodic evolution process called thermal relaxation oscillation (TRO) (Flannery 1976; Lucy 1976; Robertson & Eggleton 1977; Yakut & Eggleton 2005; Li et al. 2008). In this oscillation, contact binaries oscillate between the two states of contact and semi-detached. In order to achieve thermal equilibrium while in a contact state, energy flows from the more massive and hotter primary star into the smaller and cooler secondary star making the secondary star swell to have a larger radius compared to its mass, while the primary star becomes dimmer due to the loss of energy.

The other explanation suggested by Stepien (2007, 2009) comes from the knowledge that the primary star and the secondary star have very different status of evolution, as shown in Figs. 8 and 9. He assumes that the contact binaries are old stellar systems in the point of view of stellar evolution. According to his explanation, a detached binary of a small mass at its early stage had evolved to a semi-detached binary like Algol, which experienced mass ratio reversal, and consequently evolved to form a contact binary. At this stage, while the primary star remains in the main sequence, the secondary star, which has already developed a helium core, belongs to the sub-giant or giant stage. That the secondary star is larger compared to its mass is due to the evolution effect rather than the energy transfer. According to his model, both stars are in a thermal equilibrium state, and thus there is no oscillation like TRO. However, static equilibrium of a common gas envelope is not possible in this model. Since the common envelope has a baroclinic structure to enable a large scale transfer of mass and thermal energy between component stars at the same time, the temperature and luminosity are redistributed in the common envelope, and surface temperatures of both component stars become the same.

At this time, it is not obvious which model better explains the structure and evolution of W UMa binaries. While the TRO evolution model has a flaw that the number of semi-detached binaries predicted in the evolution process are quite small, it is also difficult to explain the frequently observed secular changes of period with the Stepien model. In order to resolve these issues, we may need more observations and a broader theoretical basis.

Although there are many disagreements in details of the evolution process while in a contact state, scientists agree that contact binaries suffer from angular momentum loss (AML) due to the magnetic-braking process via magnetized stellar winds; orbits experience continuous contraction due to AML and tidal friction; rotational angular momentum increases by orbit-spin coupling to make the mass ratio smaller and smaller; when the rotational angular momentum exceeds one-third of the orbital angular momentum, finally it is doomed to coalesce into a single star by suffering from mechanical instability of Darwin (1879) (Huang 1966; Van’t Veer 1979; Vilhu & Rahunen 1979; Mochnacki 1985; Webbink 1985; Guinan & Bradstreet 1988; Eggleton 2012). Therefore, the evolution proceeds toward the direction of a smaller mass ratio (q = M₂/M₁). In order to confirm this, the masses of the primary star and the secondary star of W UMa binaries versus the mass ratio are plotted in Fig. 10. In this figure, two solid lines show the mass variation of the primary star and the secondary star as a function of the mass ratio. The curves were obtained by least-square fits of a quadratic equation to mass variations of the primary and secondary stars, respectively. The equations corresponding to the primary star and the secondary star are as follows:

It is interesting that these two equations have the same coefficients of the second-order term (-0.47) and the constant term (-0.13). While the same coefficient of the 2^nd order term indicates that the shape of the parabola is the same, the constant term indicates that when the mass ratio is equal to 1, the masses of the primary star and the secondary star are M₁ = M₂ = 0.73 M_⊙. Although these two curves do not perfectly represent mass variation for each of the primary stars and the secondary stars relative to the mass ratio, they clearly show the overall trend of mass variation of components as a function of the mass ratio. As the mass ratio decreases, the masses of primary stars show a clear trend of gradual increase, whereas the masses of secondary stars decrease steadily. This is very important observational evidence that shows that the evolution of contact binaries proceeds toward the direction of decreased mass ratio, as indicated by the arrows in Fig. 10, from the right to the left (Stepien 2010). In particular, it should be noted that the distribution of A-subtype and W-subtype show distinct differences at q = 0.25 (log q = -0.6) and q = 0.4 (log q = -0.4). That is, 1) most contact binaries with q ≥ 0.4 are W-subtype (Only 4 out of 29 is A-subtype); 2) contact binaries with 0.25 ≤ q < 0.4 show almost an even distribution of A-subtype and W-subtype (A-subtype: 16, W-subtype: 17); and 3) among 36 contact binaries with q < 0.25, W-subtypes are only 7 and the others belong to A-subtype. This indicates that W-subtype contact binaries prevail in the early stage of evolution with a mass ratio greater than 0.4 and as the mass ratio decreases, W-subtype binaries and A-subtype binaries become equally common, and when the mass ratio is less than 0.25, A-subtype binaries are dominant. The distribution of A-subtype and W-subtype supports the assertion of various scientists that A-subtype is evolved more than W-subtype (Flannery 1976; Lucy 1976; Wilson 1978; Mochnacki 1981; Hilditch 1989), which is still debated (Maceroni & Van’t Veer 1996; Yildiz & Doğan 2013).

In Fig. 10, the mass ratio (q) of NSVS 1461538, marked by a star symbol, is 0.285 and it is at the far left side of the 0.25 ≤ q < 0.4 interval. Similar to other contact binaries of this region, it is expected that NSVS 1461538 will go through the AML process due to magnetic braking causing a gradual decrease of the mass ratio and turn into A-subtype and follow the course of evolution to form a single star consequently (Eggleton 2012).

Finally, since the absolute physical quantities obtained in this paper are based only on photometric observations, further improvement of these quantities is necessary using the radial velocity curves based on spectroscopic observations. Also, more accurate photometric observations are also needed to examine the characteristics of time variable properties of the O’Connell effect and orbital period variation in detail.