KIC 12268220: A $\delta$ Scuti Pulsating Star and an Active Protohelium White Dwarf in an Eclipsing Binary System

We present a photometric, spectroscopic, asteroseismic, and evolutionary analysis of the Algol-type eclipsing binary KIC 12268220. We find the O'Connell effect and anticorrelated eclipse timing variations in the Kepler light curve, revealing the presence of large starspots. Radial velocities and atmospheric parameters are obtained from ground-based spectroscopic observations. Combined with the radial velocity measurements and Gaia-derived total luminosity, our light-curve modeling yields the solution of the physical parameters for both the primary and secondary components. We find 14 independent frequencies arising from the $\delta$ Scuti primary, and the observed frequencies agree with the frequency range of unstable modes from nonadiabatic calculations. Based on the conclusion from previous literature, we run a grid of models to study the evolution process of our system. The evolutionary tracks of our model suggest that the low-mass ($\sim 0.23\,M_\odot$) evolved secondary shows a similar evolutionary state to the R CMa-type system, which might evolve to an EL CVn system.


INTRODUCTION
Eclipsing binary (EB) systems containing pulsating components can be used for the determination of accurate fundamental parameters. With the photometric and spectroscopic data in the time domain, the mass and radius of the EB system can be derived accurately (e.g., Southworth et al. 2005;Clausen et al. 2008). Asteroseismic modeling can constrain the parameters of pulsators such as δ Scuti and γ Dor stars (e.g., Chen et al. 2016Chen et al. , 2019, which are dwarfs or subgiants located at the lower end of the classical instability strip (Breger 2000). Many δ Scuti pulsating stars in EB systems have been discovered, especially after the Kepler mission (e.g., Guo et al. 2016;Kahraman Aliçavus , et al. 2017;Liakos & Niar-chos 2017;Gaulme & Guzik 2019). Within the group of pulsating binaries, those of Algol-type (oEA) systems (Mkrtichian et al. 2004) might experience mass transfer during their evolution.
Binary star evolution with mass transfer can generate mass-transferring or post-mass-transfer δ Scuti pulsators. One such type is the EL CVn binaries, which contain an A-or F-type primary and a low-mass helium white dwarf (WD) secondary (Maxted et al. 2014;Guo et al. 2017;Zhang et al. 2017). Currently, more than 60 EL CVn binaries are known, with 16 being discovered through the Kepler mission (Lee & Park 2018;Wang et al. 2019;Zhang et al. 2019); however only a few have pulsation signals.
In this paper, we study the activity characteristics of KIC 12268220 from its Kepler light curve and eclipse timing variations (ETV) in Section 2.1. Then, we derive the atmospheric and orbital parameters with the spectra fitting and light-curve modeling (Section 2.2 and Section 3). In Section 4, we study its pulsation properties with the nonadiabatic calculations. In Section 5, comparing to the theoretical evolution models, we suggest KIC 12268220 would evolve to an EL CVn system. Finally, we summarize our results in Section 6.

Kepler Photometry
KIC 12268220 was observed by Kepler from quarters 0 to 17 in the long-cadence mode (29.4 minute sampling) and one-month short-cadence (59 s sampling) mode. The detrended and normalized light curves from the Kepler Eclipsing Binary Catalog (KEBC 1 ; Prša et al. 2011;Slawson et al. 2011) are used in this work. Both long-cadence and short-cadence light curves are shown in Figure 1.
From the light curve we find a strong O'Connell effect (O'Connell 1951, significant flux differences at quadrature phases), which can be clearly found in the shortcadence data (bottom panel of Figure 1). This effect suggests the possible presence of starspots (e.g., Linnell 1986;Kang et al. 2004;Qian & Yang 2005). However, the long-lived O'Connell effect needs long lifetimes of starspots. The lifetimes of starspots vary significantly for different spectral types of stars (Giles et al. 2017). For solar-type stars, the lifetimes of starspots range from 10 days to a year, depending on the area of starspots (Namekata et al. 2019). The lifetimes of cooler stars and active close binary stars (RS CVn-type stars) can be up to several years (e.g., Strassmeier et al. 1994;Henry et al. 1995;Strassmeier et al. 1999;Strassmeier 1999;Giles et al. 2017). Moreover, Hussain (2002) found the spots on tidally locked binary systems live longer than spots on single main-sequence stars. Since the effective temperature ratio of KIC 12268220 can be estimated with the ratio of the eclipse depth, we can infer the secondary star should be a K-type star . Therefore, the O'Connell effect is more likely caused by the starspots on the secondary star.
Additionally, the ETV results of KEBC, KIC 12268220 display an apparent variation for both the primary and secondary eclipses. However, many secondary eclipse time variations have the same maximum or minimum values. These results may be caused by the poor fit for the secondary eclipses in their pipeline (Conroy et al. 2014). We fit the secondary eclipses with a secondorder polynomial function to look for the time of the deepest points, and then calculate new ETVs for those secondary eclipses. Combined with the ETV results of the primary eclipses from KEBC, we find an obvious anticorrelated relation (see Figure 2). This anticorrelation can be successfully explained by the moving of starspots, and the quasi-periodic variation of the amplitude in the ETV curve may imply the long-term evolution of the starspots (Tran et al. 2013;Balaji et al. 2015). However, from the simple model of Tran et al. (2013) and Balaji et al. (2015), the sum of the colatitude and inclination angle should be less than 90 • , which means there could be some polar spots in the KIC 12268220, because its inclination angle is relatively high based on its light-curve shape. Therefore, from the evidence of longlived starspots, polar spots, and short orbital period, we could infer the potentially strong magnetic field of KIC 12268220 (Schuessler & Solanki 1992;Berdyugina 2005).

Spectral Analysis
We secured seven nights from 2018 to 2019 on the 2.16 m telescope at Xinglong Station, which is administered by National Astronomical Observatories, Chinese Academy of Sciences. We obtained 12 echelle spectra with the Beijing Faint Object Spectrograph and Camera (BFOSC) E9+G10 (R ∼ 2500), G11 and G12. They are three combinations of echelle and grisms with some different wavelength ranges and resolutions (details in Fan et al. 2016). After removing three spectra with signalto-noise ratios (S/Ns) less than 30, all the spectroscopic data are reduced using the Image Reduction and Analysis Facility (IRAF) package (Tody 1986(Tody , 1993, following the standard procedures introduced by the 11th Xinglong Observational Astrophysics Training Workshop 2 . To measure the radial velocities (RVs), we apply the template fitting method with iSpec (Blanco-Cuaresma et al. 2014;Blanco-Cuaresma 2019). Before the measurements, we convert the spectrum from the air wavelength to the vacuum wavelength based on the method of Birch & Downs (1994). According to the parameters in Table 1, we choose T eff = 7800 K, log g = 3.6, and [Fe/H] = −0.3 as the initial values to generate a tem-  plate spectrum. Then, we apply the cross correlation algorithm to find the best fitted RV for each observed spectrum. iSpec can be used to analyze the double-lined spectroscopic binary. However, limited by the resolution and the luminosity difference between the primary and secondary stars, we can only measure the RVs from the primary star. Also, although most of the spectra were observed consecutively in 2018, two spectra were obtained in 2019. Thus, we add the barycentric velocity correction for each RV. The RV results are listed in Table 2. Our measurement uncertainties correlate well with the wavelength calibration results of the BFOSC (J. Zhang et al. 2020, in preparation). In order to obtain the atmospheric parameters, we combine three best-quality spectra of E9+G10 to increase the S/N. The combined spectrum is also resampled to keep the wavelength steps consistent with the raw spectrum. Then, the atmospheric parameters are derived by using iSpec with the synthetic spectral fitting technique. iSpec implements some commonly used  synthetic models and atomic line lists. In this work, considering a higher T eff (> 7000 K) and for the wavelength of the spectra, we choose the Vienna Atomic Line Data Base (VALD) line lists (Ryabchikova et al. 2015), the SPECTRUM code (Gray & Corbally 1994), the Grevesse 2007 solar abundances (Grevesse et al. 2007), and the ATLAS9 Castelli atmosphere library (Kurucz 2005). The initial parameters are set to the same values as the radial velocity calculation, and the resolution is fixed at 2500. Because of the relatively short orbital period, synchronous rotation could be expected. Thus, we specify the v rot sin i at 37 km s −1 with the estimated radius R ≈ 3.5 R and inclination angle i ≈ 70 • . iSpec applies the Levenberg-Marquardt algorithm to fit the spectrum, and the iteration stops when the ftol (relative error desired in the sum of squares) or xtol (relative error desired in the approximate solution) less than 10 −10 . The corresponding errors are calculated from the covariance matrix. However, our internal errors are probably underestimated. Thus a more robust error determination is needed. Because of the slightly higher resolution,

Parameters (Units)
Fitted Results Note-The microturbulence velocity and macroturbulence velocity are adopted by an empirical relation considering the effective temperature, surface gravity, and metallicity. The relation was constructed by the Gaia-ESO Survey.
the BFOSC was often used as the follow-up confirmation and calibration of the Large Sky Area Multi-Object Fiber Spectroscopic Telescope 3 (LAMOST; Fan et al. 2016). Therefore, we choose the mean errors of the late A-type subgiants from the LAMOST spectrum as the estimated errors. After we consider the statistic errors, the best-fitting results and their errors are listed in Table 3. The observed composite spectrum and the model spectrum are shown in Figure 3. The observed spectrum is dominated by the strong hydrogen Balmer series, without obvious peculiar metallic-line weakness or enhancement (e.g., Ca II K, Si II, Cr II, and Sr II), which matches the model spectrum well. Since we have the atmospheric parameters and RV measurements of the primary star, we can combine with the light curve to constrain the parameters of the invisible secondary star. To do so, we use the PHOEBE (Prša & Zwitter 2005), which is based on the Wilson-Devinney (Wilson & Devinney 1971;Wilson 1979Wilson , 1990Wilson & Van Hamme 2014) code, to fit the phase folded light curve in the semidetached mode. However, because the variation of starspots dominates the dispersion of the folded light curve, we use the Savitzky-Golay filter (Savitzky & Golay 1964) to obtain a smoothed version of the folded light curve, and our fitting result is based on these light curves.

LIGHT-CURVE MODELING
The effective temperature of the primary star is fixed to the value from our spectra fitting results (7843 K). We also set e = 0 with the assumption of circular orbit based on the short orbital period and the phase difference between the two eclipses of the folded light curve (Zhang et al. 2018). The albedos are set to 1.0 and 0.5 for primary and secondary stars, respectively (Lucy 1967;Ruciński 1969); the gravity brightening coefficients are adopted to 1.0 and 0.32, respectively.
In addition, as we discussed in Section 2.1, the O'Connell effect, and the anticorrelated ETV curves evidently show the existence of starspots in KIC 12268220, and the starspots could seriously affect the results of our fitting. The location, size, and temperature of a starspot are usually strongly correlated. To solve this problem, we apply some prior knowledge to the parameters of starspots. Firstly, although there are some A-type stars showing the spot-like features in the light curves (e.g., Balona 2013Balona , 2017, because of the deeper convective envelopes of the late-type stars, starspots are more likely to appear on these types of stars (McQuillan et al. 2014). We also exclude the possibility that the spots are caused by chemically abundant inhomogeneities (classified as ACV variables; Bernhard et al. 2015) through our spectrum observation. Therefore, the starspots should be on the secondary star. Second, the ratio of starspot temperature to stellar effective temperature could be derived from Berdyugina (2005) and Maehara et al. (2017) as Equation 1 T spot /T star = 1−3.58×10 −5 T star −0.249+808/T star . (1) Thirdly, although there could be multiple starspots at different locations, we only add one minimal starspot to prevent overfitting. To achieve this, we choose the colatitude of the starspot equal to the inclination angle, which means the center of the starspot is facing the line of sight. During our fitting process, the temperature ratio of the starspot and the colatitude are calculated after each iteration. Finally, with applying the differential corrections routine, we obtain a local minimum of the parameters as an initial guess.
However, without a secondary RV curve, there is still a lot of degeneracy. Since many eclipsing binaries can be used to estimate distance accurately (e.g., Guinan et al. 1998;Pietrzyński et al. 2013), to remove this degeneracy, we use the parallax from Gaia DR2 as a piece of extra information, although the errors would be significantly large. Because the flux of Kepler light curves are not calibrated, the model cannot derive the distance directly. However, the absolute luminosity of the binary system can be estimated from the Gaia DR2 parallax. To do so, the absolute magnitude can be cal- where the m V = 11.471 ± 0.02 is adopted from Everett et al. (2012) and the A V = 0.155 ± 0.062 is calculated based on the E(B − V ) = 0.05 ± 0.02 from the 3D dust map (Green et al. 2018(Green et al. , 2019 majek et al. 2015), and the BC V is the bolometric correction of V band. Since the luminosity of the invisible secondary star is significantly smaller than the primary star, the BC V = 0.05 ± 0.02 is interpolated from the MIST bolometric correction grids (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015Choi et al. 2016;Dotter 2016) with our fitted atmospheric parameters. Finally, the total luminosity is log L/L = 1.60 ± 0.04.
Then, we apply a Markov Chain Monte Carlo (MCMC) sampler to explore the posterior distribution of the binary parameters. Using MCMC to sample the posterior probability distribution is quite common to obtain more robust results of many binary systems (e.g. Schmid et al. 2015;Hambleton et al. 2018;Iglesias-Marzoa et al. 2019;Mahadevan et al. 2019). Our MCMC sampler is based on the emcee package (Foreman-Mackey et al. 2013), which is an affine invariant version of the MCMC method (Goodman & Weare 2010). Because the computational load is heavy for the MCMC method, we choose the 2015 version of the Wilson-Devinney LC code to generate model light curves and RV curves. Nine parameters are free in our sampling: the effective temperature of the secondary star, T eff,2 ; the semi-major axis, sma; the mass ratio, q; the inclination angle, i; the primary star surface potential, Ω 1 ; the center-of-mass velocity, vga; the passband luminosity of the primary star, HLA; the longitude of the starspot, xlong; and the radius of the starspot, radsp. The colatitude of the starspot is set to the inclination angle, and according to Equation 1, the temperature ratio of the starspot is restricted by the T eff,2 .
Our likelihood function is written as where the χ 2 Lum is calculated from the Gaia-derived total luminosity and the sum of model luminosity of the primary and secondary stars. The prior distribution is a uniform distribution for each parameter, and the ranges of the prior distributions are large enough to cover reasonable models. The initial parameters are obtained from the PHOEBE results. The number of walkers is 128, and to ensure convergence, we choose 30 times of the integrated autocorrelation time as the "burn-in" steps. After "burnt-in," more than 50,000 steps are restarted, and we also thin the chains with the autocorrelation time to reduce autocorrelation. Then, we can derive the final parameters and uncertainties from their marginalized posterior probability distributions. As shown in Figure 4, for each parameter, we adopt the median value as the best-fitting value, and the 16th and 84th percentiles as the upper and lower uncertainties. Moreover, the final values and errors of masses, radius, loggs, passband luminosity ratio, colatitude of starspot, and temperature ratio of the starspot are also derived from their corresponding distributions (calculated after each iteration).
The light-curve modeling results are given in Table 4, and the values without errors are the fixed parameters. The observed Kepler light curve and RV curve with the best-fitting models are shown in Figure 5. The randomly distributed residual shows our model is a good solution of the light curve and RV curve.

PULSATION ANALYSIS
From the short-cadence light curve, we can clearly see the effect of pulsation signals on the shoulders between the two eclipses. As shown in the Fourier ampli- tude spectrum (see Figure 6), the low-frequency region (f 5 day −1 ) is dominated by the orbital frequency (f orb = 0.22616 day −1 ) and its harmonics. The high frequency region shows typical δ Scuti pulsations with frequencies ranging from 20 to ∼ 24 day −1 . Indeed, our spectroscopic parameters of the primary pin point the star in the δ Scuti instability strip (see Fig. 8).
To investigate the pulsation properties, we apply the SigSpec (Reegen 2007) to the 4 yr long-cadence light curve after removing the modeled EB light curve. Because the short-cadence data shows no significant frequencies higher than the Nyquist frequency (24.510 day −1 ), we calculate the significant frequencies from 0 to the Nyquist frequency. SigSpec performs a prewhitening procedure for a given light curve. The prewhitening method calculates the Discrete Fourier Transform and fits the signal with a sinusoidal of variable amplitude and phase, then iteratively subtracts the fitted light curve from the previous light curve. In each iteration, SigSpec calculates a sig (spectral significance), which is defined as the logarithm of the inverse false alarm probability, and the false alarm probability shows that the probability of a peak is caused by pure noise in a non-equidistantly spaced data set. With the Equation  (31) of Reegen (2007), sig could be conveniently converted to the S/N. In our case, the procedure stops when sig < 5.46, which is approximately equivalent to the empirical criterion: S/N < 4. After the prewhitening, besides the low-frequencies caused by the starspots and imperfect EB fitting, 19 δ Scuti frequencies are extracted and listed in Table 5. Following the method introduced by Kallinger et al. (2008), the errors of the frequencies, amplitudes, and phases are calculated based on their sig. Because of the nonlinear effect, δ Scuti stars often show combination frequencies, we search for those combinations by computing where n and m are integers (1,2,3), and ε is the Rayleigh resolution (ε = 1/∆T ≈ 0.00068 day −1 ; ∆T ≈ 1470.46 days). If the difference between two frequencies is less than the Rayleigh resolution, the two frequencies are indistinguishable (Pápics 2012); and if a frequency could be combined by two parent frequencies with larger amplitudes, it is marked in Table 5. Finally, we obtain 14 independent frequencies.
In addition, although the pulsation properties of the primary may differ from the results of single star evolution due to the mass transfer in a binary system, we also check the frequency range of unstable modes with the non-adiabic calculation. To do so, we calculated the non-adiabaic eigen-functions and eigen-frequencies with the Dziembowski's oscillation code (Dziembowski 1971(Dziembowski , 1977) for a stellar model representing the observed parameters. The stellar model from the MESA evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018 has the following parameters: M = 2.17 M , R = 3.11 R , Z = 0.01. The mass and radius are essentially within 1σ of the observed parameters. We find that p-modes at radial orders of n p ≈ 5 − 7 are excited, having positive stability parameters (see Figure 7) that agree with the observed frequency range.

EVOLUTION AND DISCUSSION
With an A-type subgiant primary star and a low-mass (≈ 0.2 M ) evolved secondary star, KIC 12268220 might have a similar evolutionary history to the EL CVn system such as KIC 8262223 (Guo et al. 2017) and KIC 7368103 . They are formed through the case-B evolution (Paczyński 1971), which leads to a mass exchange; the initial high-mass star evolves quickly to fill its Roche lobe and transfer its mass to the lowmass secondary. Recently, Chen et al. (2017) introduced the nonconservative stable mass transfer channel, which successfully explained the formation of the EL CVn and they also showed a grid of possible parameter space.
To study the evolution of the KIC 12268220, we follow the method of Chen et al. (2017) and run a grid of models. However, solving the initial parameters from the current evolution stage is the inverse problem. It is not only time-consuming but also highly sensitive to the initial parameters. As discussed in Chen et al. (2017), the mass transfer rate, angular momentum loss, and metallicity could have significant effects on the parameter space. Therefore, we run some theoretical models with different initial parameters and only display similar evolutionary tracks for some typical final WD masses.
We use the MESA evolution code with the Ritter (Ritter 1988) mass transfer scheme and a 50% mass transfer rate; the initial metallicity is set to Z = 0.02 and the initial helium abundance is Y = 0.28. The gravitational wave radiation and magnetic braking are also switched on. In our models, we find some evolution tracks are similar to the KIC 12268220. For example, we choose the initial primary mass M 10 = 1.55 M , initial secondary mass M 20 = 1.23 M , and initial orbital period P orb0 = 2.75 days. As shown in Figure 8, the two lines in color indicate the evolutionary tracks for the primary and secondary stars respectively. This system ends up with the M 1 = 1.89 M , M 2 = 0.227 M , and orbital period P orb = 4.87 days. Although the model parameters and observed parameters of KIC 12268220 are slightly different, they are likely to experience a similar evolutionary process. Thus, we could use this model to study the evolutionary details.
From the evolutionary tracks shown in Figure 8, because the mass transfer leads to the mass and radius decrease of the primary star, the most luminous point of the primary star after first evolving from the zero age main-sequence (ZAMS) indicates the onset of mass transfer; moreover, with the end of stable mass transfer, the orbital period keeps stable, so the ending point of the mass transfer is shown as the start of the yellow range. After the typical L shape phase in the evolutionary track of the primary star, it evolves to a low-mass helium WD, with the secondary leaving its main sequence.
Therefore, for KIC 12268220, as a semidetached system, according to the position of its secondary star on the evolutionary tracks of the model, we could infer it is about to or has just finished its mass transfer. Several studies also investigate some similar cool EL CVn candidates to KIC 12268220, including KIC 10661783 (Southworth et al. 2011;Lehmann et al. 2013), KIC 8262223 (Guo et al. 2017), KIC 7368103 (Wang et al. 2019, AS Eri (Mkrtichian et al. 2004), and R CMa (Lehmann et al. 2018). The parameters of their cool secondaries are listed in Table 6.
They are classified as the R CMa-type system, which is an Algol-type system with a low-mass ratio and short orbital period (Budding & Butland 2011). Previous works believe they would most likely evolve into an EL CVn system (Lee & Park 2018;Wang et al. 2019). To compare with those cool progenitors of EL CVn, we also plot them on the HR diagram in Figure 8. Some typical evolutionary tracks of different WD masses are also plotted as references. Since the M 2 of those R CMa-type systems are around 0.2-0.21 M , we could estimate their evolutionary status from the relative position of those objects to the 0.205 M track. The objects with relatively higher temperatures are more likely to start their L shape evolution; however, the lower temperature objects are still losing their masses during the mass trans-   fer. After the mass transfer is finished, the final WD mass would be lower than the current M 2 and the orbital period would increase. Chen et al. (2017) found a tight relation between the WD mass and the orbital period (M WD -P ), which follows the formula of Lin et al. (2011). Since this relationship is determined by the degenerated core massluminosity relation, it would be nearly stable no matter what the initial parameters of a binary are. We put KIC 12268220 on Figure 9 to check its M WD -P relation. We also add the known WDs in the EL CVn systems found by Kepler (collect from Zhang et al. 2017) and the five cool R CMa objects to Figure 9. It is clear that most of the WDs follow this relation and KIC 12268220 locates near the predicted line, although the error is relatively large. From the M WD -P results in Figure 11 of Chen et al. (2017), most of the final products of their models locate on the left side of the predicted line except for a larger dispersion for M WD ≈ 0.23. However, unlike WDs, all the R CMa-type objects locate on the right side of the predicted line (red triangle), which could be confirmed by Figure 7 of Wang et al. (2019). Here we explain this with the ongoing evolution for some of the R CMa-type systems, which means they would move to the upper left on the M WD -P relation slightly after further evolution. Thus, R CMa-type system might be a transition to the EL CVn system for low-mass and short orbital period (M WD 0.22 M , P orb 3 days) objects; similar mass but longer orbital period (P orb ≈ 4 days) like KIC 12268220 would have slightly higher WD mass after evolved to be an EL CVn system. Another possible effect on the evolution is the magnetic field (e.g., Mestel 1968;Rappaport et al. 1983): a strong magnetic field would lead to more angular momentum loss. KIC 12268220 potentially possesses a strong magnetic field, which makes it an ideal target to study the role of the magnetic field during the evolution as an EL CVn precursor.

SUMMARY
In this paper, we investigate the eclipsing binary KIC 12268220 with the photometric and spectroscopic data. The spot features in the light curve suggests that KIC 12268220 has a strong magnetic activity. Combined with the atmospheric parameters, RV data, and Gaiaderived luminosity, the modeling of the light curve yields the fundamental parameters for both primary and secondary stars. The A-type primary star shows δ Scuti pulsations. We confirm the frequency range of unstable modes with the nonadiabatic theory. After running a grid of models following Chen et al. (2017), similar to the R CMa, we suggest the low-mass secondary star is a precursor of helium WD.
We appreciate the referee for several suggestive comments that significantly improved our article. We thank the help of Xinghao Chen, Gang Li, Namekata, Daniel Hey, Liang Wang, Yaguang Li, Changqing Luo and all the people who are interested in my poster at TESS Science Conference I. We also thank the Xinglong Observational Astrophysics Training Workshop for the spectrum reduction and the Simulating Stars Summer School, 2018 in University of Chinese Academy of Sciences (UCAS) for the MESA technique. We acknowledge support from the Chinese Academy of Sciences (grant XDB09000000), and from the National Science Foun-dation of China (NSFC, Nos. 11988101, 11425313, 11933004, 11903047). New data presented here were observed with the 2.16 m telescope at Xinglong Observation base in China. We thank the support of the staff of the Xinglong 2.16 m telescope. This work was partially supported by the Open Project Program of the Key Laboratory of Optical Astronomy, National Astronomical Observatories, Chinese Academy of Sciences. This work was also partially supported by funding from the Center for Exoplanets and Habitable Worlds. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. The paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission Directorate. All of the Kepler data presented in this paper were obtained from the MAST. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts.