Double eclipsing binary system KIC 3832716

TL;DR

This study reveals that the previously known eclipsing binary KIC 3832716 is actually a quadruple system with two binary pairs orbiting each other, exhibiting complex stellar interactions and activity.

Contribution

The paper demonstrates that KIC 3832716 is a quadruple system with two eclipsing binaries, providing detailed analysis of their orbital parameters and stellar evolution.

Findings

System is a quadruple with two eclipsing binaries.

Binary A shows post main sequence evolution and spots.

Binary B contains a very dim secondary component.

Abstract

In this paper we analyzed the light curve of the object KIC~3832716 observed during Kepler K1 mission. We showed that this previously known eclipsing binary is in fact quadruple system, in double eclipsing binary configuration on a long period orbit with the mass ratio of eclipsing binaries $0.7 \pm 0.3$. The system consists of eclipsing binary A with the orbital period of $P_A \sim 1.1419$\,d. Eclipsing binary A contains larger but less luminous secondary component in the post main sequence stage of its evolution consistent with an "Algol paradox" which can explain parameters of the components by mass transfer from secondary component after leaving the main sequence. Inspection of the light curve and eclipse time variations of the eclipsing binary A also indicates the presence of spots on the surface of secondary component with $57$\,d period of activity probably induced by longitudonal motion of the spots. The second eclipsing binary B with orbital period $P_B \sim 2.1703$\,d contains very dim secondary component with luminosity below 0.09\,$L_{\odot}$.

Tables

Table 1. Table 1  : Basic information about KIC 3832716 (ep=J2000, eq=2000)

Parameter	Value	Source
$RA$	${19}^h{\mathrm{\hspace{0.17em}01}}^m{\mathrm{\hspace{0.17em}34.6}}^s$	Skrutskie et al. (2006)
$DE$	${38}^{\circ }{\mathrm{\hspace{0.17em}54}}^{\prime }{\mathrm{\hspace{0.17em}17.69}}^{\prime \prime }$	Skrutskie et al. (2006)
$Kp/\left[mag\right]$	$13.42$	Conroy et al. (2014)
$T_{eff1}/\left[K\right]\text{<span><sup>a</sup><span><span><sup>a</sup>effective temperature of the primary component</span></span></span>}$	$6500\pm 400$	Armstrong et al. (2014)
$\log \left(g\right)\left(cgs\right)$	4.15	Kepler MAST
$Q_{SC}\text{<span><sup>b</sup><span><span><sup>b</sup>quarters of observational data used for each cadence</span></span></span>}$	2,4,5
$N_{SC}\text{<span><sup>c</sup><span><span><sup>c</sup>total number of data-points analyzed for given cadence</span></span></span>}$	383 532
$Q_{LC}\text{<span><sup>b</sup><span><span><sup>b</sup>quarters of observational data used for each cadence</span></span></span>}$	0-17
$N_{LC}\text{<span><sup>c</sup><span><span><sup>c</sup>total number of data-points analyzed for given cadence</span></span></span>}$	65 307
d/[kpc]dddistance to object (Gaia)	2.1(2)	Lindegren et al. (2018)

Table 2. Table 2  : Photometric solution of EB A binary.

Parameter	Primary (1)	Secondary (2)
	component	component
$T\left(K\right)$	6500	5550(350)
$r_{polar}\text{<span><sup>a</sup><span><span><sup>a</sup>values of radii are listed in semi-major axis (SMA) units</span></span></span>}$	0.20(2)	0.26(1)
$l_i\text{<span><sup>b</sup><span><span><sup>b</sup>relative luminosities of the primary and secondary component <br>of EB A</span></span></span>}$	0.48(8)	0.41(6)
$i{(}^{\circ })$	78(1)
$l_B$cctotal relative luminosity of EB B system, where $l_1+l_2+l_B=1$	0.12(9)

Table 3. Table 3  : Photometric solution of EB B binary.

Parameter	Primary (3)	Secondary (4)
	component	component
$T_{eff4}/T_{eff3}$	${0.70}_{-0.12}^{+0.04}$
$r_{polar}\text{<span><sup>a</sup><span><span><sup>a</sup>values of radii are listed in SMA units</span></span></span>}$	${0.13}_{-0.03}^{+0.02}$	${0.052}_{-0.003}^{+0.016}$
$l_i\text{<span><sup>b</sup><span><span><sup>b</sup>relative luminosities of the primary and secondary component of <br>EB B, where <xmath><xmapp><xmtok>=</xmtok><xmapp><xmtok>+</xmtok><xmapp><xmtok><xmtok>l</xmtok><xmtok>1</xmtok></xmtok></xmapp><xmapp><xmtok><xmtok>l</xmtok><xmtok>2</xmtok></xmtok></xmapp><xmapp><xmtok>⁢</xmtok><xmapp><xmtok><xmtok>l</xmtok><xmtok>B</xmtok></xmtok></xmapp><xmdual><xmref><xmwrap><xmtok>(</xmtok><xmapp><xmtok>+</xmtok><xmapp><xmtok><xmtok>l</xmtok><xmtok>3</xmtok></xmtok></xmapp><xmapp><xmtok><xmtok>l</xmtok><xmtok>4</xmtok></xmtok></xmapp></xmapp><xmtok>)</xmtok></xmwrap></xmref></xmdual></xmapp></xmapp><xmtok>1</xmtok></xmapp></xmath></span></span></span>}$	${0.964}_{-0.014}^{+0.003}$	${0.036}_{-0.003}^{+0.014}$
$i{(}^{\circ })$	${84}_{-2}^{+4}$

Table 4. Table 4  : Luminosity estimation of the components based on total luminosity calculated from Gaia and relative luminosities calculated in tables 2 and 3 .

	EB A		EB B
$l_B$	$L_1\left(L_{\odot }\right)$	$L_2\left(L_{\odot }\right)$	$L_3\left(L_{\odot }\right)$	$L_4\left(L_{\odot }\right)$
0.03	6.8	5.8	0.4	0.01
0.12	6.2	5.3	1.4	0.05
0.20	5.6	4.8	2.5	0.09

Equations

$$P_s=\frac{P_A{P}_B}{2P_A-P_B}\approx 21.8d,$$

$$P_m=\frac{P_B{P}_s}{P_s-10P_B}\approx 350d.$$

$$\begin{array}{rl}MinI=& 2455055.86515(1)+1.14187687(3)\times E\\ & +2.0(7)\times 10^{-10}\times E^2.\end{array}$$

$$\begin{array}{rl}Min{I}_B=& 2455003.9077(2)+2.1702736(3)\times E\\ & -5(1)\times 10^{-9}\times E^2.\end{array}$$

$$Q_{EB}=\frac{M_{E{B}_B}}{M_{E{B}_A}}=\frac{P_B{\dot{P}}_A}{P_A{\dot{P}}_B}=0.7\pm 0.3$$

Full text

00footnotetext: Institute of Physics, Faculty of Science, Pavol Jozef Šafárik University, Park Angelinum 9, 040 01, Košice, Slovakia

Keywords quadruple system; eclipsing binary; eclipse time variations; stellar spots

1 Introduction

The Kepler satellite, launched in 2009, has produced observations with unprecedented photometric precision (Borucki et al., 2010). It has revolutionized the study of extrasolar planets, variable stars and stellar astrophysics; providing photometric data with high-precision, high-cadence continuous light curves. After losing two reaction wheels, the Kepler spacecraft ended its primary mission and started its so-called K2 mission (Howell et al., 2014). The photometric precision of K2 is slightly lower, but still much better than that from ground-based observatories.

The analysis of excellent scientific data from Kepler can reveal not only extrasolar planets and new variable stars (Molnár et al., 2016), but detailed inspection of light curves of known systems can uncover various physical processes, such as pulsations of component(s) in eclipsing binary and/or presence of other bodies in the systems (Gies et al., 2015; Zasche et al., 2015). An example of such a system is also KIC 3832716, the analysis of which is presented in this paper.

Multiplicity of KIC 3832716 was detected during a survey for planetary candidates conducted by Borucki et al. (2011). After automatic detection of transit-like variations, the object was flagged as a false positive due to stellar nature of the eclipses. The object was included in the first release of Kepler eclipsing binary catalog (Prša et al., 2011b) as Algol-type eclipsing binary star with the orbital period of $\sim$1.142 d. In the second data release of Kepler eclipsing binary catalog (Slawson et al., 2011) another ephemeris with the period of $\sim$2.170 d was discovered. According to the catalog, the nature of the second set of eclipses was also stellar but without any further specification. The object was not listed in the list of blended sources and contamination factor of $0.030$ suggests that the light curve suffers minimally from a background light. In this paper we showed that the second ephemeris corresponds to the second eclipsing binary (EB). Additionally, Davenport (2016) detected multiple flares in the light curve using automated procedure which we confirmed by a visual inspection of the observations.

From this point forward, for the sake of clarity, we will refer to the central eclipsing binary with orbital period $\sim$1.142 d as EB A and to the second eclipsing binary with orbital period $\sim$2.170 d as EB B. In this paper we provide times of minima and ephemerides for both sets of eclipses using the method based on fitting of template curves (Sections 3.1 and 3.3) together with a method to separate both sets of eclipses (Section 3.2). Photometric parameters of EB A are derived in Section 4. In Section 7 we examine eclipse time variations and a possible explanation of these changes is given in Section 8. Our results are discussed in Section 9.

2 Observational data

For our analysis we used de-trended (PDCSAP_FLUX) data from the third revision of the Kepler eclipsing binaries catalog (Conroy et al., 2014). We used short cadence (SC) data (sampled every 58.8 s) from quarters Q2, Q4 and Q5 and long cadence (LC) data (sampled every 29.4 min) from quarters Q1 to Q17 (Gilliland et al., 2010). KIC 3832716 was observed approximately 261 d in SC regime and 1334 d in LC regime. Basic information about the object and observational data is listed in Tab. 1. De-trended light curves observed in both regimes are displayed in Fig. 1.

The initial visual inspection of the light curve of KIC 3832716 shows that it is typical for a detached EB with noticeable ellipsoidal variations between eclipses. Eclipses caused by the second eclipsing binary EB B are also visible in the original light curve during coincidences between eclipses caused by both binaries as periodic variations in depth of EB A eclipses. However, much larger variations in flux caused by EB A are preventing us from analyzing the eclipses of EB B directly from the original light curve and thus pre-whitening of the EB A signature becomes necessary (see sec. 3.3).

3 Ephemerides of eclipsing binaries

3.1 Ephemeris of EB A

Preliminary orbital period $\sim$1.14188 d of EB A was calculated using the Phase Dispersion Minimization (PDM) method (Stellingwerf, 1978). This value together with the time of the first observed primary eclipse was used to calculate the positions of all other observed eclipses. The exact times of minima were determined by fitting a template eclipse function to observed eclipses. We used the template eclipse function proposed by Mikulášek (2015) that was originally intended for phenomenological modeling of eclipsing binary light curves. Template eclipse function (eq. 13 in Mikulášek (2015)) contains 6 free parameters for each type of eclipse (primary and secondary) and they were determined by the fitting of phase stacked eclipses. We decided to fit only the light curve around eclipses in phase interval from 0.88 to 1.12 in case of primary eclipses, and from 0.38 to 0.62 in case of secondary eclipses. In order to reduce the number of fitted parameters, we performed fitting separately for both types of eclipses. We also derived different template eclipse functions for eclipses observed in LC and SC regimes.

Preliminary values of fitted parameters were obtained using genetic algorithms. Subsequently, Markov chain Monte-Carlo (MCMC) simulation with $10^{6}$ iterations was used to obtain statistically significant values of parameters.

Afterwards, template eclipse functions were used to fit particular eclipses. In this case, we fixed all parameters except parameter $A$ which represents depth of an eclipse and the phase shift parameter $\Delta\vartheta$. $\Delta\vartheta$ was introduced to the model by substituting phase parameter $\vartheta$ with $\vartheta-\Delta\vartheta$. This enabled us to shift the position of template eclipse function by the amount of $\Delta\vartheta$ and finally to determine the time of minimum for a particular eclipse. MCMC simulation with $10^{6}$ iterations was used to find parameters $A$ and $\Delta\vartheta$ for each eclipse and to properly estimate their standard deviations. Especially important was the standard deviation of the parameter $\Delta\vartheta$ from which the error in time of minima can be derived. Examples of primary and secondary eclipses of EB A observed in both regimes with their template eclipse functions are displayed in Fig. 2.

In the case of SC observations, we were able to obtain times of minima with standard deviation of approximately $1.2$ seconds for primary and $2.5$ seconds for secondary eclipses. Surprisingly, almost the same precision was achieved for times of minima of eclipses observed in LC regime. Although LC data have 30 times longer sampling period, much higher precision of LC observations and sufficient width of eclipses ($\sim$3.5 hours) were able to produce times of minima with almost the same precision for both regimes.

All determined times of minima were used to improve the orbital period obtained from the Kepler eclipsing binaries catalog and provide a linear ephemeris of EB A. O-C diagram of primary and secondary minima constructed using this orbital period is displayed in the top panel of Fig. 3.

O-C data of the EB A shows very complex structure which is the result of multiple effects. The most prominent feature is the presence of periodic spikes. The height of the spikes changes with the period of 350 d and the spikes themselves occur with the period of about 22 d.

These spikes are produced only when the eclipse of EB B is coinciding with the primary or secondary eclipse of EB A. Mutual overlap of EB A and EB B eclipses deforms overall shape of the light curve and therefore the time of minimum of EB A eclipse is shifted towards the position of nearby EB B eclipse as illustrated in Fig. 4.

The ratio of orbital periods of both binaries $P_{A}$ and $P_{B}$ (Section 3.3) is close to 1:2, thus the period of spikes can be calculated assuming this resonance of orbital periods:

[TABLE]

and the observed period of amplitude variation of spikes can be explained by 10:1 resonance between the period of spikes $P_{s}$ and orbital period $P_{B}$:

[TABLE]

Knowing the nature of the spikes, we can now remove them and investigate short term variations of the O-C diagram. Analysis of these short term variations is discussed in detail in Sections 7 and 8. We also performed quadratic fit to the O-C data without spikes to examine any long term non-linear effects. The resulting quadratic fit is displayed in the bottom panel of Fig. 3. We detected a slight increase in the observed orbital period of EB A with the rate of $\dot{P}_{A}=dP_{A}/dt\approx 0.011(4)$ s/yr. Therefore, for the final ephemeris of the EB A we can write:

[TABLE]

3.2 Phase curve smoothing and residual curve

Further analysis of the EB A requires us to significantly reduce the number of data points in the phase curve. Also, in order to study the second light curve of EB B we needed to create a residual curve which would be free from the binary signature of EB A. Both problems can be solved by using the smoothed phase curve of EB A. Smoothed phase curve was used to produce residual light curve by subtracting it from the original light curve and it was also used directly for the light curve analysis of EB A (section 4) due to the lower number of data points.

Initially, we divided the phase curve to $N_{b}=P_{A}/T_{SC}$ bins, where $T_{SC}=58.8$ s is the sampling period of SC observations. Afterwards, points of the smoothed phase curve along with their uncertainties were generated as weighted averages of flux for data in each bin. Despite lower sampling frequency, the same number of bins was also used for LC observations because of the sufficient number of evenly distributed data points across the phase curve. Problem emerging for the bins at the edges of the phase curve was solved by introducing cyclic boundary conditions, where the last data bin is immediate neighbor to the first data bin.

Due to different slopes of the phase curve during eclipses and outside of eclipses, we used different bin widths for each part of the phase curve as shown in Fig. 5 (middle). Phase width of the bins was set to 0.001 during eclipses and to 0.003 outside of eclipses. Linear transitions between two bin widths were made at the edges of eclipses to prevent the creation of artifacts. The resulting smoothed phase curve of EB A is shown in Fig. 5 (top). We tried several different widths of bins for both SC and LC regimes in order to find the optimal values for them. We found, however, that the above mentioned values of bin widths were able to smooth out all of the fluctuations of flux in the resulting smoothed phase curve below the precision of observations and simultaneously introduced the least amount of additional smearing into the data (see Fig. 5, bottom).

3.3 Ephemeris of the EB B

Smoothed light curve of EB A was subtracted from the original Kepler light curve (for LC and SC regimes separately) to obtain the residual curve, which is actually light curve of the second binary EB B. These residual curves for SC and LC observations are almost free from flux variations caused by EB A as one can see in the phase curve of EB B in Fig. 6. Using the approach of smoothed light curve enabled us to recover even eclipses which coincide with EB A eclipses and use them for the determination of the times of minima.

Minima times of EB B were calculated by a similar procedure as those of EB A (Section 3.1). One additional step has been taken before template fitting. Due to variable spot activity on the surface of EB A component (see Section 8), smoothed phase curve was not able to completely remove the signature of EB A from the residual curve. Therefore, an additional detrending of each eclipse was performed by a quadratic fit. We reached the precision of 20 s for primary eclipses and 90 s for secondary eclipses. The resulting quadratic ephemeris of the second EB B was determined:

[TABLE]

O-C diagram of EB B is shown in Fig. 7, where shortening of the orbital period $P_{B}$ with the rate of $\dot{P}_{B}\approx-0.15(4)$ s/yr is visible.

4 Parameters of EB A binary system

Smoothed phase curve calculated in section 3.2 is now suitable for the light curve analysis of EB A binary system. For this purpose we used software PHOEBE 0.31a (Prša & Zwitter, 2005). The orbital period and shape of the phase curve strongly suggests a system with circular orbit, therefore we fixed the value of eccentricity to $e=0$. Effective temperature of the primary component was fixed to value $T_{1}=6500\pm 400\,K$ (Armstrong et al., 2014). We decided to adopt this value instead of value 5926 K listed in Kepler Input Catalog (KIC) mainly due to the fact that KIC value is based on a single star model, whereas value published by Armstrong et al. (2014) was calculated considering both stars in the eclipsing binary using color information from multiple photometric surveys. We also fixed coefficients of gravity darkening $g_{1}$ = $g_{2}$ = 0.32 and the bolometric albedo coefficients $A_{1}$ = $A_{2}$ = 0.6, which are appropriate values for the convective envelopes (Prša, 2011a). The Kurucz (1993) model of stellar atmospheres was applied to the stars assuming solar composition. The limb darkening coefficients were interpolated from the van Hammes tables (van Hamme, 1993) using linear cosine law.

In our solutions, we were looking for the temperature of the secondary component $T_{2}$, surface potentials of both of the components $\Omega_{1}$, $\Omega_{2}$, orbital inclination $i$ and the photometric mass ratio $Q$. We also searched for the contribution of the second binary EB B in the system to total flux as the third light parameter $l_{B}$. The resulting parameters are listed in Table 2 and the best fit is shown in Fig. 8.

Due to possible correlation between mass ratio $Q$ and inclination $i$ we decided to perform an extensive search for other suitable solutions in $Qi$ plane around initial solution found via PHOEBE’s differential corrections method. During this search we fixed values of $Q$ and $i$ and the rest of the parameters were fitted in similar fashion as in the case of the initial solution. We searched the $Q$ parameter space in range from 0.4 to 2.5 and $i$ parameter space in range from 75.5 ∘ to 81.3 ∘. The resulting $\chi^{2}$ distribution is displayed in Fig. 9. As depicted in the Fig. 9 we were unable to determine photometric mass ratio of the system. We suspect that the main reason behind this uncertainty is the fact that the amount of third light contributed by EB B was unknown prior fit and it varied widely in our solutions across $Qi$ plane. This creates large uncertainty in absolute amplitude of ellipsoidal variations which is crucial in deriving photometric mass ratio.

Despite the previously mentioned issues, statistical analysis of models with quality of the fit $\chi^{2}<0.010$ provided us with very consistent values of effective temperature for the secondary component, polar radii and relative luminosities across the searched parameter space $Qi$. Final parameters of the binary system are listed in Table 2.

Detailed inspection revealed that out of eclipse parts of the phase curve are slightly deformed with different maxima in flux. The most probable cause of this effect is the presence of spot(s) on the surface of the secondary component. Such deformation also decreases the accuracy of photometric mass ratio, determined by our solution.

Our model indicates that EB A is currently a detached binary with tidally deformed secondary component. The radius ratio of 1.30 indicates a larger secondary component. This accompanied by lower effective temperature than the primary component suggests that we are probably seeing an evolved, post main sequence (MS) secondary component.

We decided to test if our set of effective temperatures for primary and secondary component is consistent with effective temperature provided by KIC by calculating approximate aggregate effective temperature of both components. Resulting value of 5960 K is in agreement with KIC value of 5926 K and value from Gaia archive 5821 K (Sartoretti et al., 2018).

5 Parameters of EB B

Despite large uncertainty in the light contribution of the EB B to the overall flux of the system we attempted to derive its parameters. It was assumed that all of the third light detected in Section 4 originated from the EB B and there are no other sources of light present in the system. We decided to investigate cases with different light contributions of the EB B based on the third light $l_{B}$ (Tab. 2) within its uncertainty to assess the reliability of EB B parameters. The residual light curve after prewhitening of the EB A signature (see fig. 6) was renormalized using corresponding value of the third light and the output was analyzed in PHOEBE with no other third light. Resulting parameters of the EB B are listed in Tab. 3. Parameters such as polar radii and relative luminosities suggest a system with much larger and luminous primary component accompanied by darker secondary component that contributes only 4 ‰  to the total flux of the quadruple system. Due to the same reason as in the case of EB A, uncertainty in $l_{B}$ precludes the determination of photometric mass ratio and it is presumably responsible for the lack of precision of inclination as well. Additionally, because of the lack of any color information about EB B, only ratio of effective temperatures is listed in Tab. 3.

6 Luminosity estimation based on Gaia data

We estimated total luminosity of the system using distance obtained from Gaia (see Tab.  1) and coordinate dependent band extinction parameters $A_{\lambda}$ from Schlafly & Finkbeiner (2011). Estimates were done in 2MASS J,H,K bands. Resulting luminosities ranged from 12.8 to 13.4 L⊙ with average uncertainty 2.4 L⊙.

Using 13 L⊙ as an approximate total luminosity of the system and the relative luminosities calculated during light curve analysis we calculated estimated luminosities for all components of the quadruple system for 3 different values of the third light based on value $l_{B}$ from Table 2 and its uncertainty in order to investigate the effect of uncertainty in $l_{B}$ on luminosities of the quadruple system members (see Tab. 4).

Luminosity and effective temperature of the primary component of EB A is consistent with MS star with mass around 1.5 $M_{\odot}$ based on simulations carried out in Modules for Experiments in Stellar Astrophysics - MESA (Paxton et al., 2018). However, luminosity of the secondary component of EB A is much too high for given effective temperature which puts it among evolved stars above MS as illustrated in Hertzprung-Russell (HR) diagram in Fig. 11. Uncertainty in the light contribution of the EB B prevents us from drawing any meaningful conclusions about its primary component. However, in case of the secondary component, we can conclude that obtained luminosity range is consistent with the expected luminosity range of a red dwarf.

7 Eclipse time variations

Accuracy of the minima times calculated in Sections 3.1 and 3.3 enabled us to study eclipse time variations in more detail. After removing spikes caused by coincidence of EB A and EB B eclipses (see Section 3.1), we analyzed the remaining data using generalised Lomb-Scargle periodogram (Zechmeister & Kürster, 2009) which is a very effective method for analyzing unevenly spaced data. Horne Baliunas normalization was used and we searched for significant periods in range from Nyquist period up to 1400 d. The same approach was used also for O-C data from EB B eclipses. Resulting period spectra are displayed in Fig. 12.

Periodic variations in O-C data can be a sign of the light-time effect (LITE) caused by the orbital motion of bodies already detected by eclipses or yet more bodies present in the system. However, filtered O-C data (Fig. 13) and their corresponding period spectra in Fig. 12 shows semi-regular behavior of the O-C data with no dominant period that would indicate presence of the LITE especially in EB A O-C data.

Period of 57 d was detected in the O-C data of the primary and secondary minima of EB A (indicated by gray vertical dashed lines in Fig. 12). However, no such period was detected in the data from EB B. This fact combined with different amplitudes and phase shifts of the 57 d O-C variations of primary and secondary minima are ruling out LITE as a possible explanation of this periodicity. Possible explanation for behavior of the O-C data is discussed in Section 8.

8 Presence of spots

As we showed in the previous Section 7, detected 57 d periodicity of O-C variations cannot be explained by the gravitational interaction of both binaries nor by the presence of another body orbiting central eclipsing binary EB A. Semi-regular nature of O-C variations of EB A and its irregular amplitude changes, mainly for secondary minima, strongly suggest that variations in the O-C data can be caused by the presence of spot(s) on the surface of one or both components in the central EB A.

We tested this hypothesis by subtracting template eclipse curve from observed EB A’s primary and secondary eclipses. For this purpose, template eclipse curve was obtained using the same approach as smoothed phase curve in section 3.2, but during its creation, eclipses that produced spikes were omitted. Template eclipse function from section 3.1 was not suitable in this case because of artifacts present in the residual curve after subtracting the template eclipse function. Due to the insufficient sampling period of LC data, only eclipses covered by SC data were suitable for this analysis.

Residuals of the EB A’s primary eclipses showed no significant variations. On the other hand, residuals of secondary eclipses show noticeable variations as can be seen in Fig. 14. We detected two types of variations in the residuals. The first type of variations (Fig. 14 a,c,d,f) is believed to be caused by the presence of a bright or dark spot on the surface of secondary component of EB A. The second type of variations (Fig. 14 b,e) is of a similar nature but in this case flux does not return to its previous level. It is most likely due to the fact that the spot is no longer in the line of sight because of the rotation of the secondary component. Movement of these spots across the stellar surface or periodic occurrence of spots could be responsible for detected periods in O-C data of both primary and secondary eclipses. Presence of such spots on the surface of the EB A component periodically affects the overall shape of the light curve especially during eclipses and therefore the positions of primary or secondary minima can be shifted as it was observed.

9 Discussion

Although both sets of eclipses are clearly visible in the light curve of KIC 3832716, an overall hierarchy of the system is not obvious without spectroscopic observations.

Initially we analyzed the possibility of a hierarchic system with central eclipsing binary and circumbinary third component. We tried to determine the stability of such an orbit by numerical integration using Gauss-Radau integrator RA15 (Everhart, 1985). We were unable to achieve stable prograde orbit with orbital period of only 2.17 $d$. Periodic close encounters with components of central binary caused instability in such an arrangement of bodies. Although retrograde circular orbit was stable in the time range of our numerical integrations (150 d), from the perspective of system’s evolution however, such a tight retrograde orbit is very unlikely. More importantly, the dispersion of O-C data from second set of eclipses is too low to support the possibility of such an arrangement of components. The third body would in this situation eclipse primary or secondary component which is also moving around the system’s center of mass. This would produce variations in O-C data with amplitude in the order of at least tens of minutes which was not observed (see fig. 13 bottom panel).

4-body configuration of the system is therefore the most likely explanation for these observations. In this configuration the system consists of two eclipsing binaries EB A and EB B orbiting around a common centre of mass. This arrangement of the components is much more consistent with observed data than the 3-body variant. The remaining question concerns the orbital period of the two EBs around a common center of mass. We analyzed the possibility that the detected period of 57 d could be caused by LITE effect produced by the EB B or another body present in the system. However, different amplitudes of detected periods (by a factor of 4) suggest that those variations in O-C data are not caused by the LITE. Due to the absence of any other significant common periods in EB A period spectra we can conclude that EB B or other presumed body on high inclination orbit is not located in orbit around EB A with period lower than roughly 1400 d (time span of Kepler’s observations) because such an orbit would produce observable LITE effect in O-C data and the precision of detected minima times of the EB A is sufficient to discover such objects in its vicinity due to the LITE effect.

It leaves the possibility of EB B on long period orbit around EB A and also the possibility that these two EBs are not gravitationally bound at all. The first case is however more probable due to detected non-linear effects in O-C data from both of the EBs that we described by quadratic fit in ephemerides in Eq. 3 and Eq. 4. Their opposite effect on observed orbital periods of both EBs may be explained by the orbital motion of both EBs around common center of mass. Due to very short duration of observations compared to presumed long orbital period of the EBs we did not attempt to find parameters of such an orbit. Rough estimate of the mass ratio of EBs can however in this case be calculated by using relative rate of change of orbital periods:

[TABLE]

We can conclude that this estimation of mass ratio of EBs is consistent with luminosity and mass estimations of the components in Sec. 6 assuming that both components of EB B are still on MS due to lower luminosity and insufficient time for them to evolve beyond MS.

10 Conclusion

We have identified that object KIC 3832716 consists of two eclipsing binaries with orbital periods $P_{A}\sim 1.1419$ d and $P_{B}\sim 2.1703$ d. that are in double EB configuration, orbiting a common center of mass with period much longer than the time span of Kepler photometric observations (1400 d). From derived luminosities and effective temperatures we conclude that EB A system consists of MS primary component and a more evolved post MS secondary component. However, positions of the EB A components in the HR diagram (see Fig. 11) suggest more evolved but lighter secondary component that is impossible to explain by standard stellar evolution given the sensible assumption that those two components were formed at the same time from the same cloud of interstellar matter. This situation is very reminiscent of well known case of eclipsing binary Algol (Baron et al., 2012) which is a prototype for such behavior. The situation of ”Algol paradox” was explained by mass transfer from secondary component after secondary component evolved beyond MS. We propose that the same mechanism is in action here as it can very well explain present state of the EB A, especially the relative positions of components in HR diagram and the fact that we consistently achieved very good quality of the light curve fit in case of the mass ratio between 0.5 and 0.6 for very wide range of inclinations which can be explained by unknown light contribution of the EB B system (see Fig. 9).

We also identified spot activity on the surface of secondary component, causing semi-regular variations in the observed times of minima of EB A with period about 57 d. We presume that this period is connected with some sort of cycle in spot activity or it may be the period of longitudinal travel of long term spots across the surface of the secondary component of EB A system.

EB B contains a red dwarf with luminosity below 0.09 $L_{\odot}$ orbiting around primary component on circular orbit. Due to the large uncertainty in the light contribution of the EB B we were unable to derive mass ratio of both EBs but uncertainties of other parameters such as relative radii and effective temperatures or effective temperature ratio were affected very little by this issue.

Despite the fact that analyzed data from Kepler photometry are characterized with exceptional precision and duration, additional photometric and especially spectroscopic observations are required to confirm the findings of this paper. Further photometric observations of eclipses will help to confirm the gravitational bound between the EBs, their orbital period and hopefully other parameters of the orbit. Data from spectroscopic observations in combination with photometric observations could help to determine the absolute parameters of both EBs.

Acknowledgments

This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This work was supported by the Slovak Research and Development Agency under the contract No. APVV-15-0458. The research of M.F. was supported by the internal grant No. VVGS-PF-2018-758 of the Faculty of Science, P. J. Šafárik University in Košice. The author would like to thank E. Paunzen and M. Lach for constructive criticism of the manuscript.