© ESO, 2017

1. Introduction

Stellar astrophysics is experiencing a renaissance thanks to new, highly stabilised spectrographs and very high precision photometry from space-borne observatories like CoRoT, Kepler, or MOST. The fields that, apart from extrasolar planets, benefited the most are probably asteroseismology (e.g. numerous discoveries of solar-type oscillations, or hybrid δ Sct/γ Dor pulsators), and eclipsing binaries (e.g. high-precision light curves of thousands of objects, precise eclipse timing, multi-eclipsing systems, etc.). The former allows us to look into the stellar interiors and study the structure of a star, while the latter brings directly measured, absolute physical parameters of the studied objects. Therefore, targets that combine both are extremely important for the modern astrophysical research.

Another interesting field, which is still relatively poorly studied, are multiple systems, especially those with the order of five or more. This is mainly due to low number of known systems, which decreases rapidly with the number of components. The current version of the Multiple Star Catalog (MSC; Tokovinin 1997)¹ lists >1000 triples, and 220 quadruples, but only 35 quintuples, 13 sextuples, and two septuples (the largest multiplicity order we know). There are many open questions regarding the formation, dynamical interactions, or even the abundance of such systems in the Galaxy.

In this work we present our results of a study of an object observed by the Kepler satellite, that combines all the aforementioned properties – pulsations, eclipses, and high-order multiplicity. It is, therefore, one of the most interesting targets in the field of the original Kepler mission. Since its pulsations have been studied previously by several authors, here we focus on the motion of the components (astrometry, radial velocities, eclipse timing variations) and the absolute parameters, in order to infer the evolutionary status and multiplicity of the system.

Target. Having the observed magnitudes of 7.899 in the Kepler band and 8.00 in Johnson’s V, the KIC 4150611 (a.k.a. KOI 3156, HD 181469, HIP 94924, ADS 12310 AB, WDS J19190+3916AB) is one of the brightest objects listed in the Kepler Eclipsing Binaries Catalog (KEBC; Prša et al. 2011; Slawson et al. 2011; Kirk et al. 2016)². It is usually classified as having a late-A spectral type. Its basic properties are summarised in Table 1.

According to the Washington Double Star catalogue, (WDS; Mason et al. 2001), it is a visual binary, discovered by F. Struve (no particular reference is given). The WDS currently holds 20 position measurements, the first dated 1831, and the last taken in 2015. During this time, the separation between the components A and B changed from ~18 to ~112, and the position angle from ~206° to ~203°.

The eclipses were first noted in the KEBC, where three periods are currently given: 1.5222786(19), 8.6530923(36) and 94.1982(7) d. The second value was first found in Prša et al. (2011), and is incorrectly attributed to the brighter component in the WDS. Additional eclipses were reported later by Slawson et al. (2011), and also by Orosz (2015). Multiple eclipse periods were first given by Rowe et al. (2015), who claim four values: 0.7611224(29), 1.4342003(58), 8.6530988(66), and 94.22454(40) d, which can be found in the Kepler Object of Interest (KOI) data base. The first one is obviously half of the 1.522 d value from KEBC, while the second one (1.434 d) does not appear there.

The source KIC 4150611 was first found to be a spectroscopic binary by Molenda-Żakowicz et al. (2011), but they acquired only three spectra, and did not relate the spectral lines to any particular component, nor gave the orbital solution. They were also the first to mention this object in relation to pulsations. Shortly after it was identified as containing a hybrid δ Sct/γ Dor pulsator by Uytterhoeven et al. (2011). Shibahashi & Kurtz (2012) studied the δ Sct pulsations, and used their new Fourier-transform-based formalism to obtain a radial velocity (RV) curve of the pulsating component, with the orbital period of 94.09 ± 0.11 d. Later, Balona (2014) reconstructed a somewhat different RV curve, with P = 94.3 ± 0.1 d, using a different approach. Finally, Niemczura et al. (2015) performed spectral analysis and obtained atmospheric parameters of the pulsator, showing that it is a rapidly-rotating star (vsin(i) = 128 ± 5 km s^-1) of spectral type F1 V mA9 (T_eff = 7400 ± 100 K).

2. Data and methodology

2.1. Keck II/NIRC2 adaptive optics imaging

In the Keck Observatory Archive (KOA)³ we have found high-angular-resolution images of KIC 4150611, obtained with the NIRC2 camera, which is fed by the natural guide star adaptive optics (NGS AO; Wizinowich et al. 2000) on the Keck II telescope. These observations were taken on 2013 June 15 (filter: Brγ; programme ID: U078N2; PI: Marcy) and 2014 July 07 (filters: Jcont, Hcont, Kcont; programme ID: C191N2; PI: Knutson). We have supplemented these data with our own observations, performed on 2016 September 11 in the Kcont filter, and 2016 October 15 in filters Kcont and K. Unfortunately, on our first night the conditions were bad, causing a poor AO correction. All observations were made in the “narrow” mode, which gives the field of view of ~10′ × 10′. As usual for infra-red observations, dithering was performed in order to remove the flux of the sky in further processing. Table 2 shows the observing log for AO imaging.

Except for the Brγ data, from the KOA we downloaded raw frames, and we reduced and analysed them by ourselves with a combination of IRAF and Python-based routines. Under IRAF we applied corrections for bad pixel, flat field, and sky flux. Because there are no master flats for the -cont filters available from the NIRC2 website, we took our own calibrations for Kcont and Hcont, while for the Jcont band we had to look in the KOA for appropriate calibrations. In case of Brγ frames, we used the calibrated ones, available in KOA directly. All data were then corrected for distortion. This was done with a freely-available, dedicated Python script⁴ which was prepared to follow the prescription described in Yelda et al. (2010) and Service et al. (2016), and which uses distortion maps that are available on line. We note that the maps for observations taken after 2015 April differ from those from before that date. Therefore we used different maps for the KOA (from 2013 and 2014) and our own data (2016).

The images show three point sources: two bright ones, the A and B components that form the visual binary ADS 12310 AB, and a faint one (C), located approximately 3 arcsec NW from the bright pair (Fig. 1). This star was out of the field of view of some 2013 data, due to dithering. We rejected such frames from our analysis.

For all three sources we performed relative astro- and photometric measurements. The positions of stars on the chip (x,y) were measured by fitting a 2D Gaussian to the core of the star’s point spread function (PSF). They were then translated to distances along the chip’s axes in angular scale (, Δŷ), and later, to relative angular separations and position angles (ρ, θ). For this we used the pixel scale and rotation angle values given on the same website as the distortion script, coming from Yelda et al. (2010) and Service et al. (2016). They are: 9.952 ± 0.002 mas/pix and for data taken before, and 9.971 ± 0.004 mas/pix and for data taken after 2015 April. Please note that we did not de-rotate the images, only applied the correction to the measurements of θ.

Finally, we applied corrections for the atmospheric refraction, using the method described in Hełminiak (2009), assuming monochromatic refraction at the effective wavelength of a given filter, and calculating the refractive index n with the method of Mathar (2004, 2007). For the atmospheric conditions (outside temperature, air pressure and humidity) we either used the values listed in the headers of files downloaded from KOA, or (for our own observations) we checked them in the KOA Ancillary Weather Data service⁵. These corrections were, however, smaller than the measurement errors.

The instrumental fluxes were also measured by fitting a 2D Gaussian, but the photometry was done on frames not corrected for distortion. To obtain the magnitude differences we simply translated the flux ratios using the standard formula Δmag = −2.5log (F₂/F₁).

All individual measurements were averaged, and was taken as the uncertainty. The results are shown in Table 2. Very large errors for 2016 September measurements are caused by poor AO correction. That night, the star C was extremely difficult to detect, therefore its measurements are highly uncertain, and we did not use them in the further analysis. Also, due to short integration times and the narrow bandwidth of the filter, this star was barely detectable on the Brγ images. These measurements also have relatively large error bars. The main purpose of our own AO observations was: (1) to determine if the star C is gravitationally bound to the AB pair; and (2) if it shows eclipses. We discuss this below.

2.2. Kepler photometry and light curve analysis

In this study we make use of the Q0-Q17 Kepler mission photometry, publicly available from the KEBC website. We used the de-trended relative flux measurements f_dtr, that were later transformed into magnitude difference Δm = −2.5log (f_dtr), and finally the KEBC value of k_mag was added. We call the resulting light curve (LC) the stage 0 (S0) curve. Due to the amount of data and limited computational resources, only long-cadence data were used in the LC analysis.

Before the light curve (LC) analysis, we filtered out the eclipses of the F1-type star, applying the following ephemerides: 2 544 840.707336 + E × 94.226 d, which are slightly different than those given in the KEBC, but in agreement with Balona (2014) or Rowe et al. (2015). With this procedure we also removed several eclipses that occurred in the component B, and which coincided with the 94-day period. We called this resulting light curve the stage 1 (S1) curve. The main variability in this curve are the eclipses with 8.65-day period.

For all the LC fits we used version 28 (v28) of the code JKTEBOP (Southworth et al. 2004a,b), which is based on the EBOP programme (Popper & Etzel 1981). We fitted for the period P, primary (deeper) eclipse mid-time T₀, eccentricity e, periastron longitude ω, inclination i, ratio of fluxes L₂/L₁, ratio of central surface brightnesses J, sum of the fractional radii r₁ + r₂ (in units of major semi-axis a), their ratio k, and fractional amount of the third light L₃/L_tot. With several other sources of brightness variation (pulsations, additional eclipses), the S1 curve can be treated as affected by a correlated (red) noise. Therefore, for reliable error estimation, we applied the residual-shifts (RS) method. The best fit was done on the complete Q0-Q17 long-cadence light curve, but, due to amount of data, for the error estimation we worked on single-quarter LCs. Our approach is common for all Kepler targets from our programme and is described in details in Hełminiak et al. (2016).

We first made the fit for the 8.65-day pair to the S1 curve, and run the RS to calculate the uncertainties for this system. We then took the residuals and run a Lomb-Scargle (LS) periodogram⁶ in order to identify additional periods.

We found a number of peaks related to pulsations, as well as a period of eclipses of 1.434205 d (see Sect. 3.5). We run the JKTEBOP once again with this period fixed, only to correct for these eclipses and improve the identification of frequencies of pulsation. The residuals of this intermediate JKTEBOP fit we call the stage 2 (S2) curve. We made another periodogram run on it, and we identified the most prominent pulsation periods (Sect. 3.1). We then removed the pulsations from the residuals of the S1 curve, and obtained the stage 3 (S3) curve, in which the main variability are the eclipses of the 1.43-day period. The final JKTEBOP fit for this period, the second one with proper error calculations, was done on the S3 curve. The residuals of this fit to S3 we call the stage 4 (S4) curve, in which we can see eclipses with the 1.52-day period. The last JKTEBOP run was done on the S4 curve.

Because the pulsations were mainly removed, the phase-folded S3 and S4 curves show a significantly smaller noise than the residuals of S1. We treated S3 and S4 as showing no correlated noise, and used a bootstrap approach to calculate the uncertainties.

We present the results of JKTEBOP fits to curves S1, S3, and S4 in Table 3 and Fig. 2. Apart from the parameters given directly by JKTEBOP, we also calculated the fractional fluxes of each component and used them together with the apparent magnitude of the system to calculate individual apparent magnitudes in the Kepler band (kmag). Then, using the known distance, and assuming it is the same for all three pairs, we estimated the absolute ones (Kmag). The distance modulus is Kmag−kmag = −5.56 ± 0.13 mag. Zero extinction is assumed. The sum of the individual fractional fluxes also allows us to estimate the contribution from the remaining component A to be 84.95 per cent and its absolute magnitude Kmag_A = 2.52(13) mag (assuming that the total light in the Kepler LC comes only from seven stars: the pulsator and components of three shorter-period EBs).

2.3. HIDES spectroscopy and radial velocities

We observed the target as part of our programme of spectroscopic monitoring of bright Kepler eclipsing binaries (Hełminiak et al. 2015, 2016, 2017). We used the 1.88-m telescope of the Okayama Astrophysical Observatory (OAO-1.88), with the HIgh-Dispersion Echelle Spectrograph (HIDES; Izumiura 1999), fed through a 2.̋7 diameter circular fibre (Kambe et al. 2013). Light from both visually separated components was collected. The instrument set-up, observational strategy, and spectroscopic data reduction scheme are described in details in Hełminiak et al. (2016).

Between July 2014 and October 2016, we obtained 17 high resolution (R ~ 50 000) spectra. Despite four eclipsing periods were detected in Kepler data, velocities corresponding to only two of them could be measured – of a pair of G-type stars on a 8.65-day orbit, and of the F1-type pulsator on the 94.2-day orbit. The peculiar cross-correlation function (CCF) is shown in Fig. 3. It is composed of two sharp peaks, that belong to the components of the G-type binary, and a very broad one, coming from the fast rotator of spectral type F1. This contaminator introduces additional systematic variation to the RV measurements of the other two components visible in spectra, and is the probable source of systematic uncertainties that dominate the error budget of their orbital fit.

For calculation of the RVs of the G-type pair we decided to use our own implementation of the TODCOR technique (Zucker & Mazeh 1994), which finds velocities v₁ and v₂ of two stars simultaneously. Individual RV errors were estimated with a bootstrap approach (Hełminiak et al. 2012).

The RVs of the F1-type star (component A) were measured from the position of the H_β line (4861.363 Å). We chose this line due to the fact that the Balmer series are the dominant features in the spectra, and H_β lays in a spectral range that is significantly less affected by lines from the G-type pair and tellurics than, for example, H_α. It is also in the centre of its echelle order, and the signal-to-noise ratio (S/N) around it is significantly higher than around lines at shorter wavelengths, like H_γ or H_δ. We measured its position by fitting a Gaussian to its core under IRAF’s task splot. As the measurement error, we assume a conservative value of 0.1 Å, which corresponds to ~6 km s^-1 at this wavelength. The orbital fit later showed that these errors were in fact overestimated.

2.4. RV orbital fits

The RV orbits were found with the code V2FIT (Konacki et al. 2010), which fits a single or double-keplerian orbit with a Levenberg-Marquardt algorithm. We performed two orbital fits to the measured RVs. First for the 8.65-day pair (component B), treated as a double-lined spectroscopic binary (SB2), and second for the F1-type pulsator on a 94.2-day orbit (component A) treated as a single-lined spectroscopic binary (SB1). In the SB2 fit we fixed the period and eccentricity to the values found with JKTEBOP. We fitted for velocity amplitudes (K₁,K₂), longitude of pericentre (ω), systemic velocity (γ), and moment of the pericentre passage (T_p). With V2FIT it is possible also to fit for the difference between systemic velocities of two components (γ₂−γ₁), but we initially found it undistinguishable from zero.

In the SB1 fit, we set the period to the value that we used to filter out the eclipses, 94.226 d, and also kept it fixed. When set free, the result was nearly the same, but did not reproduce the moments of eclipses very well. Other parameters were essentially the same in both cases. We also kept e fixed, as it was initially found indifferent from zero; that is, we found e < 0.17. This agrees with Shibahashi & Kurtz (2012), who found it smaller than 0.12, and with Balona (2014), who gives the value of 0.043, which is undetectable in our data. Other parameters were essentially the same, no matter if e was fixed, or fitted for. In the final fit, we were only looking for K, γ and T_p, which is defined for circular orbits as the moment of the first quadrature (highest value of RV). The value of γ in this case should be taken with some caution, as it depends on the exact reference wavelength that was used in RV measurement. This, however, has no impact on the velocity amplitude.

Results of both fits are summarised in Table 4. Observed and modelled RV curves are shown in Fig. 4. Uncertainties were estimated with a bootstrap procedure, which properly accounts for possible systematics. The individual RV measurement errors were re-scaled, so the final reduced χ² is close to one.

2.5. Eclipse timing variations

Two of our LCs – the S3 and S4 curves – were also checked for the eclipse timing variations (ETVs). The 8.65-day eccentric binary (our S1 curve) has been analysed by Borkovits et al. (2016), who did not detect any significant signal.

We used the radio-pulsar-style approach, presented in Kozłowski et al. (2011). In this method, a template LC is created by fitting a trigonometric (harmonic) series to a complete set of photometric data. Then, the whole set of photometric data is divided to a number of subsets. Their number is arbitrary, but for this study we set it to 200. For each subset, the time shift is found by fitting the template curve with a least-squares method. This approach is well suited for large photometric data sets, especially those obtained in a regular cadence, such as from the Kepler satellite.

3. Results

3.1. The hybrid pulsator and its motion

In the residuals of the first fit we see several periodicities, which are even clearer in the S2 curve. The most prominent are the δ Scuti type pulsations of the F1-type star (highest peak at 20.243 d^-1), also reported by Shibahashi & Kurtz (2012) and Balona (2014), and the γ Dor type pulsations (highest peak at 2.6064 d^-1). The δ Sct and γ Dor pulsations have similar amplitudes, and the star meets the criteria of being a hybrid pulsator (Bradley et al. 2015). The corresponding periodogram is presented in Fig. 5. An enlargement of the γ Dor area reveals a peak that coincides with the 23-rd harmonic of the 8.65-days orbital period. This mode is clearly seen in the phase-folded S1 curve (Fig. 2; top), but it is rather a coincidence than a pulsation mode induced in one of the G-type stars during close periastron passages, as it takes place in “heartbeat” (HB) stars (Beck et al. 2014). No such oscillations have been observed in main sequence stars of this type, and in the HB stars they usually have a declining amplitude. The exact frequency f₂₃ = 2.6578921 1/d correspondst to the period of 0.376238 d. Multiplied by 23, it gives 8.653474 d, which is close to but significantly different than P_B = 8.6530941(16) d. The S1 curve phase folded with the period of 8.653474 d looks much worse than with the best-fitting period.

Another small but clearly seen group of peaks can be found at frequencies around 0.185 d^-1, shown in an enlargement in Fig. 5. They most likely come from rotation of one of the components (probably the pair B), and their complicated structure reveals that the rotation is differential, and/or that they may be coming from two stars.

Both Shibahashi & Kurtz (2012) and Balona (2014) used the pulsations to detect the orbital motion of the F1-type component with the 94.2-day period. They used different approaches, which resulted in different values of asin(i) and the predicted velocity amplitude K. Shibahashi & Kurtz (2012) used the split of periodogram peaks of pulsation in the frequency domain, caused by the orbital motion, but Balona (2014) argues that such side-lobes of the main peaks may be obscured or distorted by other, independent pulsation modes of similar frequencies. He also notes that amplitude variations may also generate side lobes in the frequency domain, and decided to work with all frequencies, and look for time delay.

Closer inspection of the two results is confusing. Shibahashi & Kurtz (2012) clearly miscalculated (or made a typo in) the values of a₁sin(i) from their Table 4, giving them close to 1.2 AU, while the value expected from their RV amplitudes (K ≃ 23.5 km s^-1) would be rather 0.2 AU. Similar value can be deduced from Fig. 15 of Balona (2014), but in his Table 3, he gives a₁sin(i) = 0.140, and a (properly) corresponding K = 16.2 km s^-1. This solution is in agreement with his measurements of pulsation time delay, presented later in Fig. 16.

Thanks to our direct RV measurements, and the orbital solution (Table 4), we can clarify the situation. Our value of K₁ = 22.6 ± 2.2 km s^-1 agrees better with the one of Shibahashi & Kurtz (2012), and the corresponding a₁sin(i) is almost exactly 0.2 AU (0.196 ± 0.019). The value that can be deduced from Balona’s Fig. 15 is therefore correct.

3.2. Astrometry of the AB pair

From the existence of eclipses with the 94.2-day period we know that the inclination of this orbit is close to 90° (sin(i) ≃ 1), so the true a₁ is close to 0.2 AU. At the distance to the system, this corresponds to the projected angular value of the major semi-axis â₁ ≃ 1.55 mas. This means that we can expect a peak-to-peak astrometric displacement of ~3 mas of star A relatively to B, which should be measurable (Neuhäuser et al. 2007; Röll et al. 2008; Hełminiak et al. 2009, and the errors in Table 2). Comparison of two images, taken in 2013 and 2016, presented in Fig. 6, shows that we see a real relative motion between the two stars. The difference of ~20 mas (~2.6 AU) can not be explained by the improper distortion correction (this would increase individual errors), nor the atmospheric refraction (measurements from observations in different filters from the same night are in agreement, the scale of refraction is much smaller).

In Fig. 7 we show the measurements of the position of A relatively to B on the Δα/ Δδ plane. The left panel shows all data available from the WDS (starting from year 1831). Unfortunately the uncertainties are not given or can not be estimated in most cases. Nevertheless, one can easily note that in 185 yr the two stars approached each other, moving with the average speed of 5.26(56) × 10^-3 mas/d or 1.9(2) mas/yr in α, and 9.55(95) × 10^-3 mas/d or 3.5(3) mas/yr in δ. This seems to be the orbital motion, as the measured proper motion is over 2 times larger. One can deduce that in ~200 yr the two stars will pass very close to each other. The orbital period is most likely of the order of single thousands of years. Assuming P_AB ~ 1000 yr, and taking the estimates of total masses of A and B (m_AB,tot ≃ 4.32 M_⊙; see next sections), we can estimate the physical major semi-axis to be ~165 AU, or ~1.28 arcsec in angular units, at the distance to the system (129.4 ± 7.7 pc). This suggests a high-inclination orbit and/or significant eccentricity.

The middle panel of Fig. 7 depicts only our Keck II/NIRC2 astrometry of star A relatively to B. The orbital displacement, predicted for the time span of observations, is overplotted. One can see that the measurements do not lay on the orbital motion path, and their spread is much larger than 3 mas (predicted for the 94.2-day orbit). The orbital motion seems to be dominant here, therefore, in the right panel of the same figure, we show our measurements corrected for the orbital motion, and shifted to have the first point in (0, 0). The spread of the data is still larger than 3 mas, and the points still do not lay along a single line, which would be expected for i ≃ 90°. Moreover, when orbital phases are calculated (from the ephemeris given in Sect. 2.2), it turns out that all Keck observations were done in phases 0.1−0.3 (as labelled). If only the 94.2-day orbit was responsible for the observed positions, one would expect all the measurements to be well within ~1.5 mas. A possible explanation is another body orbiting one of the components, but it is impossible to say which one, from the relative astrometry alone. It is difficult to estimate the parameters of the putative orbit, because the orientation of the one with P = 94.2 d is unknown, and we only have precise measurements from three epochs. More high-precision astrometric data, from AO or interferometric, are needed to properly model all possible orbits. One should also keep in mind, that the long-term, gradual motion of A relatively to B (or vice versa) that we found from archival data, is quite uncertain. The archival astrometry from the WDS is not precise enough, and in most cases the uncertainties are not available. Nevertheless, it shows that the relative position of B vs. A has not changed much for nearly two centuries, so our results can not be explained by different proper motion and/or parallax of the two components. We also see no significant variations in the RV residuals.

3.3. Absolute parameters of KIC 4150611 B

Absolute values of stellar parameters of the components of the 8.65-day pair were calculated with the JKTABSDIM procedure, available together with the JKTEBOP. This simple code combines several spectroscopic and LC parameters (i.e.: P,K_1,2,i,e,r_1,2) to derive a set of stellar absolute dimensions (M_1,2,R_1,2,a), and related quantities (v_syn1,2,log (g_1,2)). Using formalism of the theory of tidal interactions, it also predicts the time scales of spin-orbit synchronisation (τ_syn), and circularisation of the orbit (τ_cir). If desired, the JKTABSDIM also calculates radiative properties (L_1,2/L_⊙,M_bol1,2) and distance, but requires both effective temperatures, multi-colour photometry, and E(B−V) as input. Due to lack of such data, we did not attempt to estimate the distance, but instead rely on the Gaia parallax (Gaia Collaboration 2016).

The parameters are presented in Table 5. Indices “Ba” and “Bb” are used instead of “1” and “2”. We reached ~1.1 per cent precision in masses, and 4.6–5.5 per cent precision in radii. The latter is hampered by additional photometric variability in the system, presence of the (dominant) third light, and the fact that in JKTEBOP analysis we did not use spectroscopic flux ratios, which help to constrain the ratio of the radii. From the parameters from Tables 4, 3, and 5, one can see that the components of the studied pair are nearly identical, both being smaller, and slightly less massive than the Sun, and the more massive component (here, the primary Ba) seems to be smaller. All the important ratios – mass, radii, and fractional fluxes – agree with unity within errors. The contribution to the system’s total flux from this binary is 14.12(45) per cent, with 6.86(6) and 7.25(6) per cent individual contributions from the primary and secondary, respectively.

3.4. The 94.2 and 1.52-day periods

In Sect. 3.1 we discussed the orbital motion of the component A (F1-type pulsator) with the 94.2-day period. As it was mentioned before, the related eclipses have quite peculiar shapes, which is due to the binary character of the object revolving around this star. Those eclipses come in two shapes: “triple” with three small dips, or “deep” with one long drop in brightness, which deepens in its central part. Examples of a triple and a deep eclipse are shown in Fig. 8.

The triple eclipse occurs when one star passes quickly in front of the pulsator, then is followed by a transit of the other star, and again by a passage of the first one, which this time is moving in the opposite direction. In the deep event, the first star starts to transit the pulsator, but the change of the observed direction (due to its orbital motion) occurs when it is still in front of the pulsator, therefore this transit is relatively long. During this eclipse the second star also transits the pulsator, which is the cause of the deeper part of the event, when larger area of the eclipsed component is obscured. Note also that during this deeper part there is a small brightening. It happens when the two transiting stars eclipse each other, and the obscured area of the star behind them is smaller. It is therefore easy to associate these brightening events with eclipses that occur with one of the other periods. They coincide with the 1.52-day period (see curve S4), meaning that this is the eclipsing pair that revolves around the F1-type star on the 94.2-day orbit. Figure 9 shows three deep eclipses (short-cadence data) phase-folded with the ephemeris for the S4 curve from Table 3. The brightenings occur exactly at phases 0.0 and 0.5.

Since we have recorded the orbital motion of the pulsator, we also attempted to do it for the 1.52-day pair, using the ETVs. This is however very difficult, as the eclipses are shallower than the rms of the curve. The expected signal has the amplitude A_Ab of the order of single minutes, but the errors of measurements themselves and their spread have larger values. Therefore we could not detect the signal securely, even though the period is known. We can only determine the detection limit of A_Ab < 283 s (rms of the ETVs), which translates into the limit of RV amplitude of the centre of mass K_Ab < 66 km s^-1 (from Eq. (8) in Hełminiak et al. 2016). Together with the RV amplitude of the pulsator itself (K_Aa = 22.6 ± 2.2 km s^-1), we can determine the limit for ratio of masses: q_A ≡ (M_Ab1 + M_Ab2) /M_Aa = K_Aa/K_Ab ≳ 0.34.

3.5. The 1.43-day period and the star C

The last period, not discussed so far, is the 1.4342 d, which is the period of the most prominent brightness variation in the S3 curve (Fig. 2). Because it was not listed in the KEBC, we were not aware of its existence, until we run the periodogram on the residuals of the JKTEBOP fit to the S1 curve. We were looking for the 1.52-day period, instead we noted the peak at 1.43420486 d (Fig. 10). The LC phase-folded with this period (and cleared from pulsations) is the S3 curve (middle panel of Fig. 2). It is difficult to associate this peak with any of the two bright visual components. The Kepler pixel mask of the target is elongated and relatively large (due to the brightness), therefore covers several nearby faint stars, but changes from quarter to quarter and the common area is small. Because the eclipses are seen in all quarters, the source must be close to the bright visual pair. Therefore we suspect that the third star seen in the AO images is also an eclipsing binary, and 1.4342 d is its orbital period.

In order to confirm that, we performed photometric measurements on the archival and our own AO observations. We discovered that the data from 2014 July 7 were taken during the beginning of the secondary eclipse, which helped us choose the observing set-up for our own observations. Before our run, we performed preliminary photometric measurements, and found that the lowest rms is reached in Hcont and Kcont filters. We decided to observe in Kcont, using exactly the same settings as in 2014, because the expected brightness variation was larger than for Hcont (for eclipsing binaries composed of stars of significantly different T_eff, the longer the wavelength the deeper is the secondary eclipse). Additionally, in October we decided to also observe the target in the K (clear) filter. This set-up is more convenient and time efficient that the previous one, that includes a neutral density filter, and requires longer integrations and calibrations.

Unfortunately, the conditions were too poor on the night of 2016-09-11, when the model predicted maximum brightness of the S3 curve, and quality photometric measurements of the star C were impossible. To make matters worse, the orbital phase during the other night (2016-10-15) was almost exactly the same as in 2014, so no significant brightness variation is observed (Fig. 11). Therefore, with our current data, we can not confirm that the star C is the eclipsing pair with the period of 1.4342 d, but we can not exclude it either. Additional observations are required, and they should optimally be taken with the K+clear set-up.

We also attempt to verify if the star C is bound to the AB system. In Fig. 12 we show how the measured positions of stars A and B relatively to C change over time, assuming that C is a distant background object that has no measurable proper motion. Our measurements do not lay on the path predicted by the parallax and proper motion of KIC 4150611 from Gaia DR1, but they also do not agree with C being gravitationally bound to AB. We already discussed the possible systematics and the difference in observed paths of stars A and B in Sect. 3.1. We conclude that C is most likely a background object, but has a measurable proper motion of few mas/year.

Finally, we have checked the ETVs of S3, obtained as described in Sect. 2.5. After removing three outliers, we run an LS periodogram on our 197 measurements. Their mean separation is ~7.35 d, therefore we looked for periods longer than 14.7 days. We detected five candidate peaks, with the strongest one at P ~ 85 d. We immediately noticed that frequencies of at least two others are multiples of the strongest one. To confirm the significance of the set of frequencies we found, we run a multi-frequency periodogram with the FREquency DEComposer algorithm (FREDEC; Baluev 2013a,b). We found that this combination of peaks is statistically significant, with all having false alarm probability below 0.05%. These results are presented in Table 6 and Fig. 13.

In Fig. 13 we also show our ETVs phase-folded with the 85-day period, and the best-fitting sine function. Despite the large scatter, the modulation with the amplitude of 804 ± 85 s is clearly seen, even if it is lower than the rms of the fit (1281 s). The amplitudes of other periods are no larger than 500 s, and when fitted for, the rms drops only to 1200 s. It is difficult to assess if these four frequencies are physically real, or just some sort of artefacts in the data. With current sampling, each subset of ~7.35 d covers about 5 orbital periods of the EB. Given the observed depth (or shallowness) of the eclipses, this seems to be a feasibility limit of the method we used – individual ETV measurement errors are quite large already. With longer sampling (e.g. 100 subsets, 14.7 d mean cadence), the ETV errors are obviously smaller. The 85 d period is still very well visible, but the other four become shorter than the approximate Nyquist cut-off (~29.4 d) and their secure detection is impossible. We then treat only the longest period as realistic.

However, assuming that this is the case, and that the modulation is caused by another body orbiting the EB, the observed amplitude can be translated into an unrealistic value of the mass function f(M) ≃ 77.5 M_⊙. Such a configuration would be difficult to explain. Therefore, the ETVs we observe for the S3 curve are rather not caused by another body. They might possibly reflect a different phenomenon, occurring on one of the stars of the whole KIC 4150611 system, such as evolution of spots, which we know exist on the Ba+Bb pair, because we see the rotational period in the periodogram in Fig. 5. They may also be a result of a systematic factor we did not take into account.

3.6. Comparison with isochrones

We compare our results – masses, radii and absolute magnitude – with the theoretical PARSEC isochrones (Bressan et al. 2012), that include calculation of absolute magnitudes in the Kepler photometric band. The estimate of iron abundance from Niemczura et al. (2015) is log ϵ(Fe) = 7.33 ± 0.10. Assuming the customary logarithmic abundances scale, with log ϵ(H) = 12, and solar iron abundance of log ϵ(Fe)_⊙ = 7.50 ± 0.04 (Asplund et al. 2009), we get the value of [Fe/H] = −0.17 ± −0.11. We assume that this represents also the metallicity [M/H] of KIC 4150611, and use PARSEC isochrones for this [M/H] value, which translates into Z ≃ 0.0103 for this set.

Comparison of our results for Ba+Bb with the models on the M/R and M/Kmag planes shows a good agreement with the 35 Myr isochrone (Fig. 14). For the estimated absolute magnitude of the pulsating star Aa (2.52 ± 0.13 mag), we can find its theoretical mass, radius, log (g), and temperature – 1.64(6) M_⊙, 1.376(13) R_⊙, 4.38(1) dex, and 8440(280) K, respectively. This is in disagreement with Niemczura et al. (2015), who give 3.8(2) dex and 7400(100) K, but it may be at least partially explained by uncertainties in the metallicity and age determinations, or the influence of additional flux on the spectral analysis. This may also be a hint for the existence of additional, relatively bright source in the system, meaning that the pulsator constitutes less than ~85 per cent of the total flux. In such situation the isochrone-predicted T_eff would be lower.

Another relatively good fit is found for a 3.16 Gyr isochrone, but in such case the component Aa would be an F5-F6 type star (5000 < T_eff < 6600 K), and this is in strong disagreement with any spectral type estimation for this star. It would also have a mass of 1.32–1.36 M_⊙. We find this scenario to be unlikely for a δ Sct pulsator, and adopt 35 Myr as the age of the system.

Taking our isochrone-based estimation of the mass of the star Aa, and the mass function value derived from the RVs (0.113 ± 0.033 M_⊙), we can estimate that the total mass of the Ab1+Ab2 pair is 0.90(13) M_⊙ (assuming sin³(i) ≃ 1). It is in a reasonable agreement with masses predicted by the isochrone from our Kmag estimates, which are 0.44 ± 0.13 and 0.32 for the primary and secondary, respectively (total of 0.77).

The agreement would be slightly better if the component Aa was fainter (cooler and less massive), which may be another suggestion of an additional flux in the Kepler LC.

At the assumed age of 35 Myr, the PARSEC isochrone predicts that the Ab1+Ab2 pair would be still in pre-main-sequence stage, with radii of 0.61 ± 0.09 and 0.52. It is therefore plausible that this pair is still very active, and responsible for the relatively weak X-ray emission detected by ROSAT (Voges et al. 2000). Taking the isochrone-predicted masses and the orbital period 1.5222468(25) d, we can estimate the major semi-axis of the orbit: . This leads to fractional radii of and . They are in agreement with the results of the JKTEBOP fit to the S4 curve (Table 3). The ETV signal of this pair, produced by the orbital motion around the component Aa with P = 94.226 d, would be ~180–210 s, which is below our detection limit.

From the isochrones, we can also estimate effective temperatures of Ba and Bb, to be 5680 and 5640 K, respectively. This is consistent with the observed G spectral type. The age of 35 Myr suggests that this pair has not reached the state of pseudo-synchronisation yet. The theoretical time scale for its masses and period is τ_syn ≃ 57 Myr (as calculated by JKTABSDIM). Therefore, the peaks in the periodogram at frequencies 0.18–0.19 d^-1 (Fig. 5) can be explained by a super-synchronous rotation of one or both components. In a pseudo-synchronous case, the frequency would be ~0.11 d^-1 (P_rot ≃ 8.7 d).

3.7. Galactic kinematics

To verify the young age of the system, we checked its Galactic velocity. Using the known parallax and proper motion (Table 1) and our value of the systemic velocity v_γ for the component B (the more reliable one, we have calculated the spatial motion components: U = 2.52 ± 0.28, V = −23.29 ± 0.16 and W = −2.49 ± 0.13 km s^-1 (no correction for the solar movement has been done). These values put KIC 4150611 well within the thin disc, probably in a moving group called Coma Berenices or “local” (Nordström et al. 2004; Famaey et al. 2005; Seabroke & Gilmore 2007). Famaey et al. (2005) have shown that ages of stars from this group vary from several to few hundreds of Myr. This confirms our “young” isochrone age of 35 Myr, and makes the “old” one (3.16 Gyr) even less probable.

4. Summary and future prospects

In Fig. 15 we present the configuration of the multiple system KIC 4150611. The orbital periods and spectral types (at the “ends” of each branch) are given to the best of our knowledge. The uncertain character of the source C is taken into account. The additional body that could explain the AB pair astrometric measurements (Sect. 3.2) is not shown.

By analysing the eclipses seen in the Kepler light curve, HIDES radial velocities and AO observations from the Keck II telescope, we obtained a consistent image of a bright, interesting multiple system KIC 4150611. We managed to directly measure physical parameters of two of the components (Ba and Bb), which allowed us to find that this is a relatively young system. From its age we inferred properties of three other component (Aa, Ab1 and Ab2). Our results are still incomplete though.

Detailed analysis of the currently available data still leaves some open questions:

• 1.

  What is the 1.43-d eclipsing binary? The faint star C seems to be a good candidate, but we can not confirm it. We lack photometry taken during the primary (deep) eclipse. Sufficient observations do not have to be made with AO facilities, but good seeing conditions and a relatively large mirror is necessary to resolve the star from the AB pair, and obtain sufficient S/N of the target. If the star C is not the 1.43-d eclipsing binary, then another source must exist in the vicinity of the AB pair.

• 2.

  What is the proper motion of the system? The star C seems to be a distant background object, so a good reference point for astrometric measurements of the AB pair, but our measurements do not match the Gaia DR1 results. The faintness of C suggests a large distance, so a very slow proper motion would be expected. The discrepancy may be caused by an incorrect determination of μ_α and μ_δ. We hope that this will be clarified with future data releases from Gaia.

• 3.

  Where does the 85-day period in ETVs come from? The signal in our ETVs of the 1.43-d binary is clear and statistically significant, but it leads to an improbable physical configuration (very high mass function). It is thus not clear if this modulation is an artefact, or has a physical origin, like evolution of spots. Confirmation would come from further timing measurements, but the source has to be resolved from the bright components A and B. This will only be possible if the star C is the 1.43-d binary. Other periods we found in the ETVs also require confirmation.

• 4.

  What is the true number of components? We presented several hints suggesting that another body may exist in the system. First, the isochrone-predicted properties (like temperature of Aa or masses of Ab1+Ab2) would be more consistent with observational constraints if the pulsator was fainter than what we found from our JKTEBOP fits. Second, as mentioned above, the 1.43-day period may not originate from the (seemingly unrelated) star C. Finally, our relative astrometry of A and B show a motion that can not be explained with the 94.2-day orbit. To verify this, more quality AO observations are required, including coronagraphic images of the surroundings of stars A and B. Precise (~1 mas level) astrometry from speckle observations, and optical or infra-red interferometers, is also welcome. Confirmation of the existence of another body would make KIC 4150611 a sextuple, or, still possibly, a septuple, which would be only the third known case.

The fact that in the system we see four periods of eclipses, with at least three coming from the system itself is intriguing. The triple sub-system A may have two co-planar or nearly co-planar orbits (94.2 and 1.52 days), as it is observed in other objects with similar architecture, for example, KOI-126 (Carter et al. 2011), HD 181068 (Derekas et al. 2011; Borkovits et al. 2012), or KIC 2856960 (Lee et al. 2013). The pairs with other two periodicities (8.65 and 1.43 days) do not, however, need to share the same orientation. Their inclinations are only calculated relatively to the plane of the sky, so the true orientations of their orbital angular momenta can still be very different. Astrometric detection of the 94.2-d motion, and complete modelling of the orbits and eclipses in the component A, will give valuable information about the distribution of momenta in this sub-system, giving important insight into the history of formation of this (young) object.

KIC 4150611 is one of the most interesting astrophysical discoveries of the Kepler mission, and we believe it deserves further attention. Additional insight into the evolutionary status and structure of the Aa component may come from detailed asteroseismic studies of its pulsations. We note that Shibahashi & Kurtz (2012) and Balona (2014) focus only on its orbital motion, and use only the δ Sct pulsations, while Uytterhoeven et al. (2011) only give the classification as a hybrid. A proper asteroseismic study is still missing. We would also like to encourage the community to perform new AO observations, in order to answer the questions described above.

---

Acknowledgments

We would like to thank Prof. Andrzej Pigulski from the Astronomical Institute of the Wrocław University, and Prof. Krzysztof Goździewski from the Toruń Centre for Astronomy of the Nicolaus Copernicus University for fruitful discussions and valuable suggestions, and Dr. Andrei Tokovinin from the Cerro Tololo Inter-American Observatory for valuable comments and corrections. This research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. The authors recognise and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. K.G.H. acknowledges support provided by the Polish National Science Center through grant 2016/21/B/ST9/01613, and by the National Astronomical Observatory of Japan as Subaru Astronomical Research Fellow. This work is supported by the Polish National Science Center grant 2011/03/N/ST9/03192, by the European Research Council through a Starting Grant, by the Foundation for Polish Science through “Idee dla Polski” funding scheme, and by the Polish Ministry of Science and Higher Education through grant W103/ERC/2011. C.B. acknowledges support from the Alfred P. Sloan Foundation.

References

Appendix A: RV measurements

In Table A.1 we present single measurements of RVs of components Aa, Ba, and Bb. For Aa, for every observation, we initially assumed equal uncertainties of 6 km s^-1. For Ba and Bb they were obtained from our TODCOR runs. Final measurement uncertainties ϵ, given in the Table, are scaled to have the reduced χ² of the fit close to one.

All Tables

All Figures