2.1. Kepler Photometry

Photometric data obtained by the Kepler spacecraft cover a broad wavelength range between 4100 ̊A and 9100 ̊A; this has the advantage of collecting many more photons in a single exposure and reaching sub-milli-mag precision, but also has the disadvantage of having no “true” photometric filter, hence no photometric color information. There are two types of photometric data having different exposure times. These are short cadence data (having an exposure time of 58.89 seconds) and long cadence data (having an exposure time of 29.4 minutes). In this study we use long cadence data of KIC 9451096 obtained from the Kepler eclipsing binary catalog. The catalog provides detrended and normalized intensities, which are obtained by application of procedures described by ^{Slawson et al. (2011)} and ^{Prša et al. (2011)}. The whole data covers ≈4 years of time, with 65307 data points in total. The MAST archive reports 0.9% contamination level in the measurements, practically indicating no additional light contribution to the measured fluxes of KIC 9451096.

2.2. Spectroscopy

We obtained optical spectra of KIC 9451096 with the 1.5 m Russian - Turkish telescope equipped with the Turkish Faint Object Spectrograph Camera² (TFOSC) at Tubitak National Observatory. TFOSC enables one to obtain intermediate resolution optical spectra inéchelle mode. In our case, the instrumental setup provides actual resolution of R = λ/∆λ ≈ 2800 around 6500 ̊A, and the observed spectra cover a usable wavelength range between 3900-9100 ̊A in 11 échelle orders. A back illuminated 2048 × 2048 pixels CCD camera, which has pixel size of 15 × 15 μm², was used to record spectra.

We obtained ten optical spectra of KIC9451096 during the 2014 and 2016 observing seasons. In order to obtain enough signal, we used 3600 s of exposure time for each observation. The estimated signal-to-noise ratio (SNR) of observed spectra is mostly between 80-100, except for a few cases, where the SNR is around 50. SNR estimation is based on photon statistic. Together with the target star, we also obtained high SNR optical spectra of HD 225239 (G2V, v_r = 4.80 km s^-1) and ι Psc (HD 222368, F7V, v_r = 5.656 km s^-1), and adopted them as radial velocity and spectroscopic comparison templates.

We reduced all observations using standard IRAF³ packages and tasks. A typical reduction procedure starts with obtaining a master bias frame from several bias frames taken nightly, and subtracting the master bias frame from all object, calibration lamp (Fe-Ar spectra in our case) and halogen lamp frames. Then the bias corrected halogen frames are combined to form an average halogen frame and this average frame is normalized to unity to produce the normalized master flat frame. After that, all target and calibration lamp spectra are divided by the normalized flat field frame. Next, cosmic rays removal and scattered light corrections are applied to the bias and flat corrected frames. At the end of these steps, reduced frames are obtained and these frames are used for the extraction of spectra. In the final steps, Fe-Ar frames are used for wavelength calibration of the extracted spectra and the wavelength calibrated spectra are normalized to unity by using cubic spline functions.

3.1. Radial Velocities and Spectroscopic Orbit

The first step of our analysis is to determine the radial velocities of the components and the spectroscopic orbit of the system. We cross-correlated each observed spectrum of KIC9451096 with spectra of template stars HD 225239 and ι Psc, as described in ^{Tonry & Davis (1979)}. In practice we used the fxcor task in IRAF environment. We achieved better cross-correlation signals (especially for the weak secondary component) when we used HD 225239 as template; thus, we determined all radial velocities with respect to the HD225239 spectrum. We obtained acceptable cross-correlation signals of both components inéchelle orders 5 and 6, which cover a wavelength range between 4900-5700 ̊A. Figure 1 shows the cross-correlation functions of two spectra obtained around orbital quadratures.

Fig. 1 Cross-correlation functions of two spectra obtained around orbital quadratures. The letter ϕ denotes corresponding orbital phase. P and S indicate the primary component and the secondary component, respectively. 

We list the observation log and the measured radial velocities of the components in Table 1. Note that we use the ephemeris and period given by ^{Borkovits et al. (2016)} and listed in their Table 2 to calculate orbital phases and for further analysis.

We achieved a reasonable solution for the spectroscopic orbit assuming zero eccentricity, where an undefined longitude of periastron is taken. We checked this assumption by inspecting the Kepler light curve of the system, where we observe deeper and shallower eclipses at 0.0 and 0.5 orbital phases, respectively, indicating a circular orbit (see § 3.3, Figure 4). In order to reach the final spectroscopic orbital solution, we prepared a simple script written in Python language, which applies Markov chain Monte Carlo simulations to the measured radial velocities, considering their measured errors. We list the final spectroscopic orbital elements in Table 2 and plot the measured radial velocities, their observational errors, the theoretical spectroscopic orbit and residuals from the solution in Figure 2.

Fig. 2 (a) Observed radial velocities of the primary and the secondary (blue and red filled circles, respectively), and their corresponding theoretical representations (blue and red curve). (b) Residuals from theoretical solution. The color figure can be viewed online. 

3.2. Spectral Type

We rely on our intermediate resolution TFOSC optical spectra to determine the spectral type of the components. Most of our spectra correspond to the phases around orbital quadratures, where we observe the signal of the two components separated. However, there are two spectra obtained at phases close to the eclipses, where the two components can not be resolved separately. One of these spectra corresponds to ≈ 0.56 orbital phase (see Table 1), where we cannot observe the radial velocity signal of the secondary component in cross-correlation. Even at the orbital quadratures, the cross-correlation signal of the secondary component is considerably weak compared to the primary component, indicating a very small light contribution from the secondary component to the total light of the system. Our preliminary light curve analysis shows that the contribution of the secondary component to the total light does not exceed ≈ 10%. In this case, the signal from the secondary component becomes almost negligible at the resolution of our observed spectrum at ≈ 0.56 orbital phase. Therefore, we assume that we only observe the spectrum of the primary component and adopt this spectrum as reference spectrum for the primary component. We confirm this assumption by calculating the composite spectrum of the binary via final parameters of the components (see § 3.3), where we observe that the contribution of the secondary component affects the theoretical composite spectrum less than 2% for the wavelength range of 4900-5700 ̊A. We refrain from performing a detailed analysis with spectral disentangling. Future studies could take advantage of this technique and derive atmospheric parameters of the secondary.

We first compare the reference spectrum with the template spectra of HD 225239 and ι Psc. We observe that ιPsc spectrum provides a closer match to the reference spectrum but also indicates earlier spectral type and slightly lower metal abundances for the primary component. At that point, we switch to the spectrum synthesizing method. We use the latest version of python framework iSpec (^{Blanco-Cuaresma et al. 2014}) which enables practical and quick calculation of a synthetic spectrum with a given set of atmospheric parameters via different radiative transfer codes. Among these codes we adopt the SPECTRUM⁴ code (^{Gray & Corbally 1994}), together with ATLAS-9 (^{Castelli & Kurucz 2004}) model atmospheres and the actual line list from the third version of the Vienna atomic line database (V ALD3, ^{Ryabchikova et al. 2015}).

Considering the spectral type of ι Psc, we synthesize spectra for effective temperatures between 6000 K and 7000 K in steps of 250 K, and metallicity values ([Fe/H]) between −1.0 and 0.0 in steps of 0.5. For all synthetic spectra we fix the gravity (log g) to 4.15, which we precisely calculate by light curve modelling (see § 3.3). Since we do not have a high resolution spectrum, we fix the microturbulence velocity to 2 kms^-1. We convolve all calculated spectra with a proper Gaussian line spread function in order to degrade their resolution to the resolution of the TFOSC spectra. Instrumental broadening in TFOSC spectra is 2.2 ̊A, corresponding 119 km s^-1 for wavelengths around 5500 ̊A. The estimated projected rotational velocities of the components are 62 km s^-1 and 36 km s^-1 for the primary and the secondary component respectively (see § 3.3). Since instrumental broadening is the most dominant broadening source in the observed spectra, we do not consider rotational broadening and other line broadening mechanisms.

Among the calculated spectra we find that the model with 6500 K effective temperature and an [Fe/H] value of −0.5 provides the closest match to the reference spectrum. The final effective temperature indicates F5 spectral type (^{Gray 2005}). Considering the effective temperature and metallicity steps in model atmospheres, and the resolution of TFOSC spectra, the final values and their estimated uncertainties are T_eff = 6500±200 K and [Fe/H] = −0.5±0.5 dex, respectively. Note that even if we considered the neglected contribution of the secondary component in the reference spectrum, its effect would be within the estimated uncertainties above. The final T_eff value is ≈ 670 K lower than the 7166 K value estimated in ^{Armstrong et al. (2014)}. In Figure 3 we show portions of the reference spectrum and the model spectrum, calculated with the final parameters above.

Fig. 3 Representation of the observed (black), best matched (red) synthetic spectrum and residuals (blue) for three regions. Note that we shift the residuals upwards by 0.3 for the sake of simplicity. 

3.3. Light Curve Modelling and Physical Properties

Global visual inspection of KIC 9451096 Kepler photometry reflects properties of a typical close eclipsing binary. We start the light curve modelling by phasing the whole long cadence data with respect to the ephemeris and period given by ^{Borkovits et al. (2016)}, and re-binning the phased data with a phase step of 0.002 via the freely the available fortran code lcbin⁵ written by John Southworth. We plot the binned and phased light curves of the system in Figure 4, panels a and aa. The light curve indicates a detached configuration for the system. Mid-eclipse phases are 0.0 and 0.5 phases, indicating a circular orbit. There is no conspicuous asymmetry in the light curve.

Fig. 4 (a) Phase binned light curve of KIC9451096 (black filled circles) together with best-fitting model (red curves). (b) Close up view of the light curve at light maxima. c) Residuals from the best-fitting model. Panels at right (aa, bb and cc) are the same as left panels but for phased long cadence data. The color figure can be viewed online. 

We model the light curve with the 2015 version of the Wilson-Devinney code (^{Wilson & Devinney 1971}; ^{Wilson & Van Hamme 2014}). In the modelling, we first fix the most critical two parameters of the light curve modelling, i.e., the mass ratio (q) of the system and the effective temperature of the primary component (T₁). Since we have reliably derived these parameters in previous sections as q = 0.55 and T₁ = 6500 K, we adopt them as fixed parameters. The calculated atmospheric properties of the primary component reveal that both stars have convective envelopes. Therefore, we set albedo (A₁ , A₂ ) and gravity darkening (g₁, g₂) coefficients of the components to 0.5 and 0.32, respectively, which are typical values for stars with convective outer envelopes. We also consider a slight metal deficiency of the system, and thus adopt the internal stellar atmosphere formulation of the Wilson-Devinney code according to the determined [Fe/H] value of −0.5. We assume that the rotation of the components is synchronous with the orbital motion, and thus fix the rotation parameter of each component (F₁ , F₂ ) to 1.0. We adopt a square root law (^{Klinglesmith & Sobieski 1970}) for limb darkening of each component; this is more appropriate for stars cooler than 9000 K. We take the limb darkening coefficients for the Kepler passband (x₁, x₂, y₁, y₂) and the bolometric coefficients (x_1bol , x_2bol , y_1bol , y_2bol ) from ^{van Hamme (1993)}. In the modelling, we adjust inclination of the orbit (i), temperature of the secondary component (T₂), dimensionless omega potentials of the components (Ω1, Ω2) and luminosity of the primary component (L₁). We also include a phase shift parameter as adjustable in the modelling, since we expect a shift in the ephemeris due to the light-time effect of the third body (^{Borkovits et al. 2016}). The model quickly converged to a steady solution in a few iterations. We list the model output in Table 3 and we plot the best-fitting model in Figure 4, panels a, b, and the residuals from the model in panel c.

In Figure 4, panel b, one can easily see the model inconsistency around 0.25 orbital phase. The inconsistency indicates an additional light variation, which is known as O´Connell effect, i.e. difference between light levels of subsequent maxima in an orbital cycle. Possible sources of the difference may be Doppler beaming, a hot spot or a cool spot on one of the component of the system. KIC9451096 is a detached eclipsing binary, thus we can safely exclude possibility of mass transfer between components, i.e., a hot spot possibility. Doppler beaming was detected observationally among some Kepler binaries (see, e.g. ^{van Kerkwijk et al. 2010}), which becomes important for systems with very low mass ratios, especially for systems with a compact component, such as a white dwarf or a hot sub-dwarf. In addition, if the effect is in progress, then it would change the light levels of each maxima. However, we observe inconsistency only for phase 0.25, while the model fairly represents the light level at phase 0.75. Thus, Doppler beaming should have a negligible effect in the case of KIC 9451096, if any. A remaining possibility is cool spots located preferably on the cooler component.

Here we do not chose to model this inconsistency alone, which would only show the cumulative effect of hundreds of light curves, but instead we subtract the best-fitting model from the whole long cadence data and inspect the residuals in order to investigate further light variations. We will focus on this in § 3.4.

We complete the light curve modelling section with a calculation of the absolute parameters of the system by combining the spectroscopic orbital solution and light curve model results. In Table 4, we give the physical properties of each component. Our analysis reveals that the system is formed by an F5V primary and a K2V secondary component.

3.4. The Out-of-Eclipse Variations

In this section, we subtract the best-fitting light curve model from the whole long cadence data and obtain residuals. Here, we first divide the whole long cadence data into subsets, where each subset covers only a single orbital cycle, resulting in 1026 individual light curves. Then we apply the differential corrections routine of the Wilson-Devinney code and fix all parameters, except the ephemeris reference time. In this way, we find a precise ephemeris reference time for each individual subset, therefore eliminating any shift in the ephemeris time due to the third body reported by ^{Borkovits et al. (2016)}, and obtain precise residuals. In Figure 5, we plot three different parts of the residuals. Note that we remove data points that correspond to the eclipse phases due to the insufficient representation of the model at those phases. This mainly arises from the inadequacy of radiative physics used in light curve modelling for a very high photometric precision and can clearly be seen in Figure 4 panel c.

Fig. 5 (a) Residuals from whole long cadence data. Remaining panels show different time ranges of residuals, where we observe different light curve shapes, and flares. 

Inspecting residual brightness, we immediately see a variation pattern which changes its shape from time to time. Furthermore, we observe a sudden increase and gradual decrease in the residual brightness which occasionally occurs over four years of time span and has short time scale of a few hours. These patterns are traces of magnetic spot activity, which is very possible for the K2V secondary component. Observational confirmation of this possibility can be done by inspecting magnetic activity sensitive spectral lines, such as the Hα and Ca II H & K lines. We inspected these lines in our TFOSC spectra and did not notice any emission features, which could be considered as the sign of the activity. However, one should consider that the contribution of the secondary component to the total light does not exceed 10% at optical wavelengths and will steeply decrease towards the ultraviolet region of the spectrum. Furthermore, the variation patterns observed in Figure 5 exhibit very small amplitudes. Therefore, the existence of magnetic spot activity cannot be confirmed or excluded via spectral line inspection in the case of KIC 9451096. Nevertheless, variation patterns and flares observed in the residuals indicate weak magnetic spot activity in the secondary component, which can still be detected with the very high precision of the Kepler photometry.

We analyze rotational modulation and flares of the secondary component via residuals by assuming that the source of all variation patterns is only the secondary component.

3.4.1. Photometric Period and Differential Rotation

Conventional periodogram methods for determining rotational period do not perfectly work in our case because the observed variation patterns exhibit quick changes in amplitude and mean brightness level over short time scales of a few days, which is comparable to the orbital period. Moreover, since we remove data points at eclipse phases, this causes regular gaps in the data which repeat each ≈ 0.625 day (i.e. half of the orbital period); thus, it causes an alias period and its harmonics, and disturbs the real periods. Furthermore, one can clearly see that the rotational modulation of residuals has an asymmetric shape. Considering an individual light curve with an asymmetric shape, it is not possible to find a single period to represent the whole light curve perfectly, and additional periods (i.e. harmonics) are required. Therefore we apply an alternative method based on tracing the time of a minimum light observed in an orbital cycle, which was previously applied to RS CVn system HD 208472 (^{özdarcan et al. 2010}). For each orbital cycle, we find the time of the deepest minimum in the cycle by fitting a second or third order polynomial to the data points around the expected minimum time. The order of the polyno Subset mial depends on the light curve shape. After obtaining all minimum times, we construct an O − C diagram by adopting the first minimum time observed in the residuals as initial ephemeris reference time, and the orbital period as the initial period, and obtain O−CI values. Then we apply a linear fit to the O−CI values and calculate an average ephemeris reference time and period given in Equation 1, together 6 with statistical uncertainties given in parentheses for the last digits.

In the equation, T₀(BJD) and E denote ephemeris reference time and integer cycle number, respectively. We plot O − CI values and linear fit in Figure 6, panel a. After obtaining an average ephemeris and period, we subtract the linear fit from O−CI data and obtain O−CII data, which in principle shows the real period variation for a given time range. Figure 6, panel b shows O − CII data. We divide O − CII data into 30 subsets by grouping data points that appear with a linear slope. The linear trend of a subset gives the difference between the best-fitting photometric period of the subset and the grand average photometric period given in Equation 1. Therefore we can calculate a mean photometric period for each subset. We plot the calculated mean photometric periods versus time in Figure 6, panel c, together with the statistical uncertainties. We list photometric periods for 30 subsets in Table 5, and tabulate our O − C analysis results in Table 8.

Fig. 6 (a) O−CI diagram of observed minimum times (blue filled circles) and linear fit (red line). 

The average period given in Equation 1 represents the average rotation period for magnetic activity features on the surface of the secondary component, which are typically cool and dark regions, i.e., star spots, and indicates a slightly (∼0.5% day) shorter period compared to the orbital period. This is clearly observed in Figure 6 panel c, where the mean photometric periods of subsets are mostly shorter than the orbital period. Assuming a solar type differential rotation, this means that the orbital period is slightly longer than the equatorial rotation period of the secondary component. Under the same assumption, the differential rotation coefficient can be estimated from (P_max −P_min)/P_equ = kf, where P_max , P_min , k and f denote observed maximum and minimum period, differential rotation coefficient and a constant that depends on the range of spot forming latitudes, respectively (^{Hall & Busby 1990}). Considering the small amplitude of rotational modulation of residuals, we assume that the secondary component is not largely spotted and that the total latitudinal range of the spot distribution is 45 degrees, which causes the f constant to take values between 0.5 and 0.7 (^{Hall & Busby 1990}). Using maximum and minimum photometric periods from the O − C analysis, and assuming that the shortest period corresponds to the equatorial rotation period of the star, we find k = 0.081 ± 0.011 and k = 0.058 ± 0.006 for f = 0.5 and f = 0.7, respectively. Since these k values are calculated via boundary values of f, the real differential rotation coefficient must lie in the range of k values calculated above. An average k is found as 0.069±0.008.

3.4.2. Flares

We detect 13 flares in the residuals from long cadence data. In the flare analysis, it is critical to determine the quiescent level, which denotes the brightness level in the absence of a flare. In our case, we determine the quiescent level by applying Fourier analysis to the single orbital cycle where the flare occurs. The Fourier analysis represents the rotational modulation of residuals in the cycle, and then we remove the Fourier representation from the data. The remaining residuals show only the quiescent level and the flare itself. We show such a flare light curve in Figure 7.

Fig. 7 An example of a flare light curve. The filled black circles represent the observations, while the red line represents the quiescent level derived from the data out-of-flare. The color figure can be viewed online. 

The energy (E) is a very important parameter for a flare. However, the energy parameter has the luminosity L of the star as a factor in equation E = P × L described by ^{Gershberg (1972)}. Due to the disadvantages described in ^{Dal & Evren (2010)}, we use the flare equivalent duration instead of the flare energy, which is more proper. We compute the equivalent durations of flares via the equation P = [(I_flare − I₀)/I₀]dt (^{Gershberg 1972}), where P is the flare equivalent duration in seconds, I₀ is the quiescent level intensity, and Iflare is the intensity observed at the moment of the flare. Considering the quiescent level, the times of flare beginning, flare maximum and flare end are determined, together with flare rise duration, flare decay duration and flare amplitude. We list all computed values in Table 6 for each of the 13 flares.

^{Dal & Evren (2010}, ²⁰¹¹⁾ suggest that the best function to represent the relation between flare equivalent duration and flare total duration is the OPEA, where the flare equivalent duration is considered on a logarithmic scale. The OPEA function is defined as y = y0 +(Plateau−y0)×(1−e−kx), where y is the flare equivalent duration on a logarithmic scale, x is the flare total duration, and y0 is the flare equivalent duration in the logarithmic scale for the least total duration, according to the definition of ^{Dal & Evren (2010)}. It should be noted that the y0 does not depend only on the flare mechanism, but also depends on the sensitivity of the optical system used in the mission. The most important parameter in the model is the Plateau value, which defines the upper limit for the flare equivalent duration on a logarithmic scale and is defined as the saturation level for a star (^{Dal & Evren 2011}). Using the least squares method, the OPEA model leads to the results in Table 7. We plot the resulting model in Figure 8 with its 95% statistical sensitivity limit.

Fig. 8 The OPEA model obtained over 13 flares. The blue filled circles show each flare while the continuous red line shows the OPEA model and the dotted red lines show the sensitivity range of the model. The color figure can be viewed online.  

We tested the derived model by using method proposed by ^{D’Agostino & Stephens (1986)} to understand whether there are any other functions to model the distribution of flare equivalent durations on this plane. In this method, the probability value (P value), is found to be ≈ 0.10, which means that there is no other function to model the distributions (^{Motulsky 2007}; ^{Spanier & Oldham 1987}).

^{Ishida et al. (1991)} described a frequency for the stellar flare activity as N₁ = Σn_f /ΣT_t, where Σn_f is the total flare number detected in the observations, while ΣT_t is the total observing duration from the beginning of the observing season to the end. In the case of KIC9451096 we find the N₁ frequency as 0.000368411 h^-1 adopting the total long cadence observing duration as 1470.2786 days from the times of the first and last long cadence data points.