The first light curve modeling and orbital period change investigation of nine contact binaries around the short period cut-off

Contact binaries are generally composed of two Roche lobe overfilling F, G, or K type stars. Although they are very common in the universe, their formation, evolution and ultimate fate are still poorly understood (e.g., Bradstreet & Guinan 1994; Qian 2003; Yakut & Eggleton 2005; Stepien 2006; Eggleton2012). Rucinski (1992) firstly noticed a 0.22 day short period cut-off by using the data of contact binaries in the 4th edition of the General Catalogue of Variable Stars (GCVS). More recently, Drake et al. (2014a) presented the period distribution of the potential contact binaries in the Catalina Data Release 1 while Qian et al. (2017) showed the distribution of orbital period of contact binaries based on the Variable Star Index (VSX) and the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) databases. Li et al. (2019) have carried out statistics on the period distribution of contact binaries according to the VSX and GCVS databases and most results of the worldwide photometric surveys. All these period distribution studies show that an obviously sharp decline at about 0.22 days can be identified.

Since the discovery of this short period limit, researchers have tried to solve this issue. Rucinski (1992) claimed that the fully convective limit can possibly explain the short period cut-off. Stepien (2006, 2011) suggested that the timescale of the angular momentum loss (AML) is too long to form extremely short period contact binaries. Jiang et al. (2012) proposed that the low mass limit at about 0.63 $M_{\odot}$ of the primary component causes the short period limit. Very recently, Li et al. (2019) claimed that a third companion is very important during the formation of some ultrashort period contact binaries (USPCBs¹¹1contact binaries with orbital periods shorter than 0.23 days. ) by removing the AML from the central eclipsing star. At present, the short period limit issue is still an open question, more contact binaries around the short period cut-off should be observed and analyzed. In this paper, we present the first photometric investigations and orbital period analyses of nine USPCBs. Table 1 shows the basic information of these nine targets.

The observations of these USPCBs were carried out with five meter-class telescopes from 2016 to 2019. These five telescopes are the Weihai Observatory 1.0-m telescope of Shandong University (WHOT, Hu et al. 2014), the 85 cm telescope at the Xinglong Station of National Astronomical Observatories (NAOs85cm) in China, the 60cm Ningbo Bureau of Education and Xinjiang Observatory Telescope (NEXT), and the 84cm and 2.12m telescopes of the San Pedro Martir Observatory (SPM84cm and SPM2.12m). The PIXIS 2048B CCD was used for WHOT, the Andor DZ936 CCD was used for NAOs85cm, the FLI PL23042 CCD was used for NEXT, the Spectral Instruments 1 CCD was used for SPM84cm, and the Marconi 4 CCD was used for SPM2.12m, resulting in the field of views as $12^{{}^{\prime}}\times 12^{{}^{\prime}}$, $32^{{}^{\prime}}\times 32^{{}^{\prime}}$, $22^{{}^{\prime}}\times 22^{{}^{\prime}}$, $7.6^{{}^{\prime}}\times 7.6^{{}^{\prime}}$, and $6^{{}^{\prime}}\times 6^{{}^{\prime}}$ for the five telescopes, respectively. During the observations, the standard Johnson-Cousins filters were employed. Observing details of all stars including the observing date, exposure time, observational errors, the comparison and check stars, and the telescope are listed in Table 2. The standard aperture photometry reduction technique was applied to process all observing data by using IRAF routines¹¹footnotemark: 1 ²²footnotetext: IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy Inc., under contract to the National Science Foundation.. After the data reduction, the differential magnitudes between the target and a comparison star were derived. In order to confirm that the brightness of the comparison star was constant, a check star was also introduced. Figure 1 shows the light curves of J031700 as an example. As seen in the lower panel of Figure 1, the different magnitudes between the comparison and check stars are flat, thus indicating that the brightness of the comparison star is constant.

The light curves of the nine stars were all analyzed by using the 2013 version of Wilson-Devinney (W-D) code (Wilson & Devinney 1971; Wilson 1979, 1990, 1994). First, we have to estimate reliable effective temperatures for these stars. The temperatures can be estimated by using their color indices. All stars lie in distances of several hundred pc, up to more than one kpc. In addition, some of them have small galactic latitudes. Hence, the color indices are very likely reddened. Therefore, the interstellar extinction for each star should be considered when estimating its temperature. The interstellar extinctions of different filters for each target were determined by using the method provided by Schlegel et al. (1998) based on the IRAS database¹¹footnotemark: 1 ³³footnotetext: https://irsa.ipac.caltech.edu/applications/DUST/. For example, the values of E(B-V) for J003033, J014418, J031700, J074350, J104942, J130945, J145224, J161857, and J224015 were determined to be 0.469 mag, 0.055 mag, 0.148 mag, 0.048 mag, 0.033 mag, 0.013 mag, 0.044 mag, 0.055 mag, and 0.056 mag, respectively. The temperatures derived by different dereddened color indices and spectroscopic observations are shown in Table 3. T_BV is the temperature due to the dereddened $B-V$ (APASS DR9, Henden et al. 2015) color index and Table 5 of Pecaut & Mamajek (2013), T_JK is the temperature due to the dereddened $J-K$ (2MASS, Cutri et al. 2003) color index and Table 5 of Pecaut & Mamajek (2013), T_gi is the temperature due to the dereddened $g-i$ (SDSS DR12, Alam et al. 2015) color index and Table 3 of Covey et al. (2007), and T_spec is the temperature due to the spectroscopic observation by LAMOST (Luo et al. 2015) or Drake et al. (2014a). Gaia temperatures (Gaia Collaboration et al. 2018) were not used in our study because of irrespective of interstellar extinction. The last column of Table 3 displays the average values of the temperatures derived by different ways, the corresponding errors represent the standard deviations. Actually, this average value is neither the primary’s nor the secondary’s, it should be a quadrature temperature of the two components. Therefore, when modeling the light curves of these systems, we set this temperature as the primary’s (primary means the hotter component), $T_{1}=T_{m}$, and the final temperatures of the primary and secondary components were derived by using the following equations (Kjurkchieva et al. 2018):

where $\Delta T=T_{1}-T_{2}$ and $c=L_{2}/L_{1}$ can be determined by the adopted final photometric results of each target. Since the temperatures of all the systems are less than 5200 K, the gravity darkening and bolometric albedo coefficients were fixed as $g_{1,2}=0.32$ and $A_{1,2}=0.5$, and bolometric and bandpass limb-darkening coefficients were automatically determined from Van Hamme (1993) by adopting a logarithmic function (ld=-3).

Because all the nine targets were studied for the first time, a $q$-search method was used to determine their mass ratios. Since our light curves show that all the systems are contact binaries, mode 3 (contact mode) of the W-D code was used during the analysis. When searching for the mass ratio, some parameters were adjustable: the orbital inclination $i$; the effective temperature of star 2, $T_{2}$; the monochromatic luminosity of star 1, $L_{1}$; and the dimensionless potential ($\Omega=\Omega_{1}=\Omega_{2}$). The relationship between the resulting sum of weighted square deviations ($\Sigma W_{i}(O-C)^{2}_{i}$) and $q$ for each system is displayed in Figure 2. As seen in this figure, every system has a minimum value of $\Sigma$. This value was set as an initial value and an adjustable parameter, then new solutions were performed until a convergent solution was obtained. The determined parameters for each system are displayed in Table 4, and the corresponding synthetic light curves, solid black lines, are plotted in Figure 3. We found that some of the targets manifest asymmetric light curves, in these cases spot model was applied to fit their light curves (e.g., Qian et al. 2007; Lee et al. 2013; Zhou et al. 2016a). In addition, the search of a third light ($l_{3}$) was carried out for every system. The derived results are also shown in Table 4, and the corresponding theoretical light curves with spot and including third light, respectively, labeled with blue and red lines are displayed in Figure 3. We found that four stars exhibit spot activity and three stars reveal third light. We should state that when we searched for the third light, the spot parameters were fixed to quickly obtain a convergent solution. When we determined the final photometric results of all stars, the temperatures of their two components were calculated using Equation (1) and are listed in Table 5. The errors were derived by error propagate.

The period variation study of eclipsing binaries is a powerful tool to investigate their dynamical evolutions and to search for additional companions (e.g., Zhou et al. 2016b; Er-gang et al. 2019; Liao et al. 2019). To analyze the orbital period variations of the nine targets, we tried to collect as many eclipsing times as possible. Some worldwide photometric surveys, such as Northern Sky Variability Survey (NSVS, Wozniak 2004), Wide Angle Search for Planets (SuperWASP, Butters et al. 2010), LINEAR (Palaversa et al. 2013), Catalina Sky Surveys (CSS, Drake et al. 2014b), and All-Sky Automated Survey for SuperNovae (ASAS-SN, Shappee et al. 2014; Jayasinghe et al. 2018) have already observed most of our targets. However, the time resolution of NSVS, LINEAR, CSS, and ASAS-SN is very low. We constructed the complete phase or half of phase light curve using longer time span data. Using these light curves, the eclipsing times were calculated with the K-W method (Kwee & van Woerden 1956) and are displayed in Table 6. Figure 4 displays an example for deriving the complete phase or half of phase light curve. First, we divided the CSS light curves of J073450 into four parts as shown in the upper panel of Figure 4. Second, using the equation, $MJD=MJD_{0}+P\times E$, where $MJD$ is the observing time, $MJD_{0}$ is the reference time for each part, and $P$ is the orbital period of J073450, we moved the long time span data into one period and determined the four lower panels of Figure 4. Then, we obtained seven eclipsing times from these four sets of light curves. Since the observing times of NSVS, LINEAR and CSS are given in $MJD$, we firstly transformed $MJD$ to $JD$ using $JD=MJD+2400000.5$, and then converted to $HJD$. Four targets (J003033, J031700, J104942, and J161857) have been observed by SuperWASP, and Lohr et al.(2013b) have obtained their eclipsing times. We acquired the eclipsing times of the four targets thank to a private communication with Dr. Lohr. Using the following equation,

we computed the $O-C$ values for all the systems. In this equation, $HJD$ is the time of minimum light which is listed in Table 6, $HJD_{0}$ is the initial epoch which is shown in Table 1, and $P$ is the orbital period which is also shown in Table 1. The $O-C$ diagrams of all targets are displayed in Figure 5. For J031700, we found an obvious linear variation from cycle -12372 to -5097 (eight points) which is shown in a dashed line square box. This means that the period of J031700 is incorrect. A linear equation was applied to fit the eight points and a new period was determined to be 0.2256396(3) days. Then, we used this orbital period to reconstruct the $O-C$ diagram of J031700 which is also plotted in Figure 5. Now, we can see that all the $O-C$ diagrams manifest long-term orbital period changes. The following second-order polynomial,

was applied to fit the $O-C$ diagrams. The results are shown in Table 7. We have found that five systems exhibit long-term increasing orbital period and the others manifest secular decreasing orbital period.

Based on our new complete light curves of nine USPCBs, we first determined their photometric elements by using the W-D code. We found that all the targets are W-subtype contact binaries. Eight of the targets are shallow contact systems ($f\leq 25\%$) while one (J130945) is a medium contact system. Four stars show clearly O’Connell effect, which can be modeled by dark spot or hot spot on one of the two components, and three systems have been detected with third light. By using all available eclipse timings, we analyzed the period variation behavior of all targets and found that four of them exhibit long-term orbital period decrease, while five of them display long-term orbital period increase.