2.1. Spectral Observation

We carried out spectroscopic observations of the system with the 1.5-m Russian - Turkish telescope and the Turkish Faint Object Spectrograph Camera (TFOSC¹) at Tubitak National Observatory. Using the échelle mode of TFOSC, we achieve an actual resolution of R=λ/∆λ ≈ 3000 around 5500Å. The whole spectrum covers a wavelength interval between 3900 - 9100Å in 11 échelle orders. All spectra were recorded with a back illuminated 2048 × 2048 pixels CCD camera, which has a pixel size of 15 × 15 µm ².

We recorded eight optical spectra of KIC8043961 between the 2014 and 2017 observing seasons. In addition, we obtained optical spectra of ι Psc (HD 222368, F7V, v_r = 5.656 kms⁻¹) and o Aql (F8V, v_r = 1.40 kms⁻¹) with the same instrumental set-up. The spectra of these two stars were used as templates for radial velocity determinations and preliminary spectral type estimations. Typical signalto-noise ratios (SNR) of observed spectra are between 140 and 190 for a single observation with an hour of exposure time.

We reduced and extracted the observed spectra by following standard steps for reducing and extracting échelle spectra. In the first step, bias correction was applied to all images, and a normalized flatfield image was obtained via bias-corrected halogen lamp frames. Then, all science and Fe-Ar calibration frames were divided by the normalized flat-field frame; scattered light correction was applied to all bias and flat-field corrected images. At that point, reduced science and Fe-Ar images were obtained. Next, spectra were extracted from the reduced images. In the following steps, wavelength calibration was applied to the reduced and extracted science frames via Fe-Ar images, and finally wavelength calibrated science images were normalized to the continuum by using cubic spline functions. The whole procedure was applied under IRAF² environment.

2.2. Photometric Data

In the light curve and O−C analyses, we mainly use long cadence (29.4 minutes of exposure time, hereafter LC) Kepler photometry, which is available in the Kepler eclipsing binary catalogue³ as detrended and normalized fluxes (^{Slawson et al. 2011}; ^{Prša et al. 2011}). We also use the short cadence (hereafter SC) data from the same database in the same format, in order to investigate possible high frequency variations. The data were collected over a broad wavelength interval between 4100Å and 9100Å, which suffers from lack of color information, but provides very high precision. The LC data cover ≈ 4 years of time base without any considerable time gap. The reported light contamination is 0.3% in the MAST database, practically meaning no light contribution of any light sources from the background. Available SC and LC photometry is shown in Figure 1. In the figure, we computed the phase-folded light curve with respect to the light elements in equation (1) given by ^{Slawson et al. (2011)}.

Fig. 1 KIC 8043961 light curve derived with the data taken from the Kepler Database. In the upper two panels, the light variation is plotted versus time and phase depending on available LC Data; it is plotted for the SC data in the bottom two panels. 

3.1. Radial Velocities and Spectroscopic Orbit

We start our analysis with the measurement of the radial velocities of each component, in order to determine the spectroscopic orbit of the system. For this purpose, we cross-correlate each observed spectrum of KIC 8043961 with spectra of ι Psc and o Aql. The cross correlation was applied according to the method proposed by ^{Tonry and Davis (1979)}. The practical application of the method was done with the FXCOR task under IRAF environment. Among the two template stars, ι Psc provided a slightly better cross-correlation and was adopted as the radial velocity template. Cross-correlation was performed for the wavelength range between 4900 - 5400Å, ignoring strongly blended spectral lines. In Figure 2, we show the cross-correlation functions of spectra obtained at orbital quadratures, where both components are clearly resolved.

Fig. 2 Cross-correlation functions of KIC 8043961 obtained at orbital quadratures. The letter φ shows the orbital phase in each panel. P and S denote the primary and secondary component, respectively. The color figure can be viewed online. 

We tabulate the measured radial velocities in Table 2, together with brief observational information. In the table, orbital phases are calculated via ephemeris time and period given by ^{Borkovits et al. (2016)}. These light elements are adopted for phasing the radial velocities and light curve.

Before attempting an orbital solution, we inspected the orbital eccentricity via phased folded Kepler photometry. Primary and secondary eclipses are located precisely at 0.0 and 0.5 phases, respectively, indicating a non-eccentric orbit. Therefore, we model the spectroscopic orbit assuming a circular orbit. The orbital solution was achieved by applying the Levenberg-Marquardt algorithm (^{Levenberg 1944}), and Markov Chain Monte Carlo simulations via EMCEE software developed under a Python environment (^{Foreman-Mackey et al. 2013}). We tabulate the final spectroscopic orbit parameters and their uncertainties in Table 3, and we plot the phasefolded radial velocities, theoretical models and residuals from the best-fitting models in Figure 3. We notice that the less massive, smaller and hotter component is eclipsed at the primary minimum, hence we interpret the numerical values in Table 3 according to this finding.

Fig. 3 (a) Measured radial velocities of the primary and the secondary component (blue and red filled circles, respectively), and corresponding model orbit (blue and red curve). (b) Residuals from the model. The color figure can be viewed online. 

3.2. Spectral Type

In order to estimate atmospheric properties of the components reliably, we use the TFOSC spectrum of KIC 8043961 obtained at HJD 2456842.3645, corresponding to ≈ 0.75 orbital phase, where we can clearly resolve spectral lines of both components. The SNR of that spectrum is around 170, based on photon statistics. Before searching for the best atmospheric parameters, we fixed the micro-turbulence velocity of each component at 2 kms⁻¹, and the logarithm of the surface gravity log(g) of the primary and secondary components at 4.04 and 3.77, respectively. These log(g) values are calculated from combination of spectroscopic orbit parameters and light curve modeling results (see next section). Then we adjusted effective temperatures and metallicities of the components via the spectral synthesis method in a wavelength interval 4750 - 5750Å. In practice, we used the latest version of the Python framework iSpec (^{Blanco-Cuaresma et al. 2014}). Synthetic spectra were calculated for temperatures between 6250 - 7250 K with steps of 50 K, and for metallicities between solar (0.0) and −1.0 with steps of 0.25. Among the different radiative transfer code options provided by iSpec, we adopted SPECTRUM⁴ code (^{Gray and Corbally 1994}) along with ATLAS9 model atmospheres (^{Castelli and Kurucz 2004}) and a line list obtained from the third version of the Vienna atomic line database (VALD3) (^{Ryabchikova et al. 2015}).

All calculated synthetic spectra were convolved with a Gaussian line spread function to match their resolution to the TFOSC spectra. After we obtained the instrumentally broadened synthetic spectra, we started to calculate the composite spectrum of the system. At that point, we selected a trial synthetic spectrum for the primary component. Next, we selected a synthetic spectrum for the secondary component, with a proper temperature as indicated by the preliminary light curve modelling. Then, we applied a proper radial velocity shift to each selected synthetic spectrum to match its spectral lines to the lines of the corresponding component in the observed spectrum. Finally, we calculated the composite spectrum of the system via synthetic spectra mentioned above, and the radius ratio, R₁ /R₂ = r₁ /r₂ = 0.69, where r means fractional radius, found from the light curve analysis. We repeated this process for every possible synthetic spectra configuration and searched the best matched composite spectrum of the system, i.e. the best-fitting synthetic spectrum for each component. We checked the reliability of the parameters suggested by the best matched synthetic spectra, with the best-fitting light curve model parameters. Confirming the light curve and spectroscopic modelling results iteratively, we find the selfconsistent effective temperature of the primary and the secondary component as 6900 K and 6600 K. The temperatures and the surface gravities suggest F3V and F5 III-IV spectral types for the hot and cool components, respectively (^{Gray 2005}). The resulting overall metallicity of the system is found to be [M/H] = 0.25. Estimated uncertainties are 200K in temperature and -0.25 in metallicity, which are predicted by considering the temperature and metallicity steps in the model atmospheres, the SNR of observed spectra and the instrumental resolution of TFOSC spectra.

Furthermore, we observe a perceivable effect of rotational broadening in our spectra. Full width at half-maximum (FWHM) of the instrumental profile of TFOSC spectra is ≈1.8Å around 5500 Å, meaning that any projected rotational velocities (v sini) above ≈ 50 km s⁻¹could be observed. From the combined spectroscopic orbit and light curve modelling, we estimate v sini values of the primary and the secondary component as 81 km s⁻¹and 54 km s⁻¹, respectively, assuming synchronous rotation. Since numerical fitting attempts of v sini velocities did not yield stable and reliable results, we only applied rotational broadening to the corresponding synthetic spectrum of the components, by using the estimated velocities above, and then calculated the composite spectrum of the system. In Figure 4, we plot observed spectrum, best fitting composite synthetic spectrum and residuals.

Fig. 4 Observed spectrum (black line), best-fitting composite synthetic spectrum (red line) and residuals (blue line). Residuals are shifted upwards by 0.4 for a better view. The color figure can be viewed online. 

3.3. Orbital Period Variation

In order to study the orbital period variation of KIC 8043961, we determine mid-times of primary and secondary eclipses, where we obtain 1688 minima times in total, of which 844 are the primary minima and 844 the secondary times of minimum. For each minimum, we only consider the orbital cycle that includes the corresponding minimum, and

then run the differential correction program of the Wilson-Devinney code by only adjusting ephemeris time and keeping all other parameters fixed. This process is fairly straightforward for primary eclipses, while we need to change the role of the components during the iterations, in order to compute secondary eclipses. Ten computed minima times are listed in Table 4, while the rest of them are available upon request. For each observed eclipse time, we compute expected eclipse times by using the following light elements (equation 1) to obtain O − C values:

Inspecting the O − C values, we notice a clear sinusoidal variation, which is presumably due to the light travel time effect (hereafter LTTE) caused by an unseen third body in the system. Taking the light elements given in equation (1) as initial values, we modelled the sinusoidal variation together with the linear trend. Corrected ephemeris time and period are given by equation (2):

where ∆t is term for the LTTE of the third body. Using the corrected ephemeris time and the period, we calculate (O−C) _II residuals. Applying the LTTE equation given by ^{Irwin (1959)} to the (O−C) _II residuals, we calculate the orbit of the third body. In Figure 5, we show (O−C) _II residuals together with the best-fitting third body orbit. Calculated third body parameters are listed in Table 5.

Fig. 5 The O-C variations of KIC 8043961. In the figure, the (O-C) II variation, derived after the linear correction, is shown in the upper panel, together with the best-fitting orbit (smooth red curve), while the residuals are shown in the bottom panel. The color figure can be viewed online. 

Assuming that the orbital inclination of the presumed third body is the same as the inclination of the eclipsing binary, and using the calculated total mass of the eclipsing binary (2.892 M⊙) for the absolute parameters of the components, we estimated the mass of the third body as listed in Table 5 for different i values.

_{3.4. Light Curve Analysis}

To determine the physical conditions of the components, we model the average SC light curve by using 2015 version of the ^{Wilson - Devinney Code (Wilson and Devinney 1971}; ^{Wilson and Van Hamme 2014}). The shape of the light curve clearly indicates a detached configuration. In the analysis, we fix the two critical parameters, mass ratio and effective temperature of the primary component (T₁) to 1.09 and 6900 K, respectively, which are found from the spectroscopic orbit modelling and the spectral type estimation. Considering the spectral type corresponding to this temperature, the albedos (A₁ and A₂) and the gravity-darkening coefficients (g₁ and g₂) of the components are adopted for the stars with convective envelopes (^{Lucy 1967}; ^{Rucinski 1969}). The nonlinear limb-darkening coefficients (x₁ and x₂) of the components are taken from ^{Van Hamme (1993)}. In the analyses, the temperature of the secondary component (T₂), the dimensionless potentials of the components (Ω₁ and Ω₂), the fractional luminosity of the primary component (L₁), the inclination of the system (i), the mass ratio of the system (q), the semi-major axis (a) and the third light are adopted as adjustable parameters. Although the synthetic curve derived by these parameters covers the entire observed light curve, a stable sinusoidal variation is still seen at out-of-eclipses. Considering the temperatures of the components, we modelled this variation by a stellar cool spot on the secondary component. We list the best-fitting model results in Table 6, and show the averaged phase-folded light curve together with the model in Figure 6.

Fig. 6 Phase-folded LC light curve (filled circles) taken from the Kepler Database, and the best-fitting model (smooth red curve). The color figure can be viewed online. 

The combination of results obtained from both the spectroscopic orbit and light curve modelling yields the physical parameters of the system listed in Table 7. These parameters show that the smaller, hotter and less massive primary component is still on the main sequence, while the larger, cooler and more massive secondary component appears to have evolved off the main sequence. We show the position of the components on the log Teff-log(L/L⊙) plane in Figure 7. In the bottom panel of the figure, the evolutionary status of the components is shown in the same plane with isochrones for [M/H] = −0.25 (^{Bressan et al. 2012}). The average age of the system is about 2.2 ± 0.3 Gyr.

Fig. 7 Upper panel. Positions of the primary (blue asterisk) and the secondary (red asterisk) components of KIC 8043961 on the logTeff-log⁡L/L⊙ plane. The primary appears to be located among the γ Doradus type stars; the ZAMS and TAMS were taken from ^{Girardi et al. (2000)}. The small filled black circles represent γ Doradus type stars listed in ^{Henry et al. (2005)}. The dash-dotted lines (red) represent the borders of the area, in which γ Doradus stars appear. The blue dashed lines indicate the theoretical cool and hot boundaries of the γ Doradus instability strip, given by ^{Warner et al. (2003)}. In addition, we plot the hot (HB) and cool borders (CB) computed from ^{Rolland et al. (2002)} of the δ Scuti stars for comparison. The small filled purple circles represent detached and semi-detached binaries taken from ^{Soydugan et al. (2006)}. Bottom panel. Individual components in the same plane with three different isochrones (^{Bressan et al. 2012}). The filled circle shows the less massive, smaller but hotter component, while the red filled circle shows the more massive, larger but cooler component. The green, red and purple isochrones are for logage of 9.40, 9.35 and 9.30, respectively. The color figure can be viewed online. 

3.5. The Variation at Out-Of-Eclipses

Inspecting the residuals from the light curve modeling, we find a clear signal of a sinusoidal variation at out-of-eclipses. Considering the place of the primary component of KIC 8043961 in the HR diagram, the source of the variation is possibly a γ Doradus type pulsation. We analyse the residual data to search for possible pulsation frequencies, using the PERIOD04 software (^{Lenz and Breger 2005}) that is based on the discrete Fourier transform (hereafter DFT) method (^{Scargle 1982}). As seen in Figure 8, we found 137 frequencies which exceed the 3 σ level. We tabulate these frequencies in Table 8, and plot the resulting model obtained from these frequencies together with the residual light curve in Figure 9.

Fig. 8 KIC 8043961’s Fourier amplitude spectrum derived from the DFT. 

Fig. 9 The variation at out-of-eclipses and the synthetic model derived by the frequencies obtained from the DFT. The filled circles represent the observations, while the red line represents the model. The color figure can be viewed online.