Detection of two additional circumbinary planets around Kepler-451 (original) (raw)

Journal Article

,

Department of Astronomy and Space Sciences, Faculty of Science, Ankara University

, 06100 Ankara,

Turkey

Search for other works by this author on:

,

Department of Astronomy and Space Sciences, Faculty of Science, Ankara University

, 06100 Ankara,

Turkey

Search for other works by this author on:

,

Department of Astronomy and Space Sciences, Faculty of Science, Ankara University

, 06100 Ankara,

Turkey

Search for other works by this author on:

Department of Astronomy and Space Sciences, Faculty of Science

, Istanbul University, 34119 Istanbul,

Turkey

Istanbul University Observatory Research and Application Center

, 34119 Istanbul,

Turkey

Search for other works by this author on:

Received:

20 January 2022

Revision received:

03 February 2022

Accepted:

04 February 2022

Published:

10 February 2022

• PDF
• Split View
• • Article contents
  • Figures & tables
  • Video
  • Audio
  • Supplementary Data
• Cite

Cite

Ekrem Murat Esmer, Özgür Baştürk, Selim Osman Selam, Sinan Aliş, Detection of two additional circumbinary planets around Kepler-451, Monthly Notices of the Royal Astronomical Society, Volume 511, Issue 4, April 2022, Pages 5207–5216, https://doi.org/10.1093/mnras/stac357
Close

Navbar Search Filter Mobile Enter search term Search

ABSTRACT

We announce the detection of two new planetary-mass companions around Kepler-451 binary system in addition to the one detected previously based on eclipse timing variation analysis. We found that an inner planet with 43 d period with a minimum mass of 1.76 _M_jup and an outer one with a ∼1800 d orbital period with a minimum mass of 1.61 _M_jup can explain the periodic variations in the residuals of the one-planet fit of the eclipse timings. We updated the orbital period of the middle planet as 406 d, and determined its eccentricity as 0.33. The newly discovered outer planet is also on an eccentric orbit (0.29), while the innermost planet was assumed to have a circular orbit. All three Jovian planets have similar masses, and our dynamical stability test yields that the system is stable.

1 INTRODUCTION

Eclipsing post-common envelope binary systems with a subdwarf B-type primary (sdB) and an M-dwarf secondary (dM) offer a favourable observational opportunity to determine eclipse timings within very good precision. Besides the short-lived eclipse events owing to the small sizes of its components, the primary eclipse, in particular, is also deep, symmetrical, and V-shaped due to the large temperature contrast making the measurements of the profile centre even more efficient. Therefore, the variations in their eclipse timings, also known as ETVs, can be studied in more detail based on more precise and accurate data sets of minima timings. Moreover, additional bodies with substellar masses might have formed in an accretion disc built up by the material lost from the upper layers of the progenitor of the sdB component (Zorotovic & Schreiber 2013) in the short duration following the ejection of the common envelope. In fact, such potentially ‘second-generation’ planets have been suggested in multiple sdB+dM binaries (HW Vir: Lee et al. 2009, NN Ser: Beuermann et al. 2010, NY Vir: Lee et al. 2014). Nevertheless, Bear & Soker (2014) came up with an alternative explanation for the survival of the first-generation circumbinary (CB) planets from the ejection of the common envelope too. An orbital decay can be expected from such systems due to a resonant interaction between the binary and a CB disc (Chen & Podsiadlowski 2017), the evidence of which was observed in the NN Ser system (Hardy et al. 2016).

Kepler-451 (KIC 9472174) is an sdB+dM eclipsing binary displaying a strong reflection effect on the secondary caused by the very high surface temperature of the primary on its light curves. Østensen et al. (2010) was first to analyse the Kepler light curve of the system as well as its radial velocities. They reported the temperature of the primary as 29 564 K from their spectroscopic analyses. They also derived the mass of the sdB component by making use of the surface gravity value (log g) from spectroscopy in a mass-radius calibration as 0.48 M⊙, which is the canonical mass for an sdB pulsator. Based on this value, they were able to compute the secondary mass as 0.12 M⊙ from the light-curve analysis.

The first planet detected in the Kepler-451 system was announced by Baran et al. (2015) based on its ETVs. The suggested Jovian-mass planet was reported to have an orbital period of 416 d about the common centre of mass with the binary pair. In a following study, Krzesinski et al. (2020) were not able to detect the signal and attributed that found by Baran et al. (2015) to a consequence of the data extraction and processing. Therefore, the planet’s existence is disputed. Subsequently, a research report from Baran & Silvotti (2020) claimed the discovery to be valid by recalculating the timings and showing that the ETV signal is indeed periodic.

Kepler space telescope observed the system in its primary mission continuously for almost 4 yr. This data set forms the backbone of all the studies of the system. Baran et al. (2015) also reported ETVs either due to a secular or a long term cyclic change that could not be covered by the time span of the Kepler data. The variation in that time range resembles an upward parabola, which is not expected considering the lack of mechanisms that would lead to either mass-loss from the secondary or mass-accretion by the primary. On the other hand, another potentially bound, yet unseen CB object can explain the observed part of a variation with a period longer than the time span of the observations. Nevertheless, interactions between the objects and a potential CB disc could be more complicated than theorised.

With the arrival of TESS observations of the system, we investigated this question that long waited to be answered. We modelled the light curves from Kepler and TESS, as well as from our own observations to extend the wavelength coverage and to track the ETVs in the large data gaps between TESS and Kepler observations. We used these models to calculate the mid-eclipse times, which we also derived from the well-established Kwee-van Woerden (KW) technique (Kwee & van Woerden 1956) for comparison, and performed a detailed ETV analysis of the system. We were able to reveal two additional periodic signals in ETVs, which we attribute to two newly found planetary-mass bodies orbiting the binary pair. These discoveries make Kepler-451, the first system with three CB planets in an sdB+dM system, and second to Kepler-47 (Orosz et al. 2012) when all CBs are considered. Finally, we carried out a dynamical stability test for this interesting system, which can open up a new window in the discussion on their formation, whether the numbers, distances, and total masses of the detected bodies can be explained by the ejection of the common envelope material or not, assuming the observed ETVs are due to the suggested unseen second-generation additional bodies. We also investigate if potential magnetic-activity induced orbital period variations can be an alternative or complimentary explanation (Applegate 1992).

2 OBSERVATIONS AND DATA

We gathered Kepler and TESS data through Space Telescope Science Institute’s (STScI) Mikulski Archive for Space Telescopes (MAST)1 portal service. The Kepler data set has a time range from 2009 to 2013. We used only the short cadence pre-search data conditioned simple aperture photometry (PDC-SAP) data that are more suitable for the light-curve analysis compared to the long cadences in which the eclipses are not even detectable due to the short orbital period. We removed the trends on each observational quarters by applying third-degree polynomial fits. Since both ends of almost all the quarters have strong linearity issues, we had to discard corresponding data points as well as random outliers before the detrending procedure. We normalized each cycle to their individual median value (see Fig. 1). For the light curve analysis, we phase-folded these data and binned them to have equally spaced 300 points for an orbital cycle.

Figure 1.

Detrended and normalized light curves of Kepler-451 that were used for producing the phase-folded data. Top panel shows the light curves from Kepler. The bottom left hand shows TESS light curves from sector 14 & 15, and the bottom right panel shows sector 40. Note that only Kepler and TESS sector 14 & 15 were used for the light-curve analysis. See the text for more details.

TESS observed the system in sectors 14, 15, and 40. We applied the same detrending, phase-folding, and normalization procedures that we did for the Kepler data, but this time on each sector separately. The trends in the two parts of each sector separated by a gap due to data downlink seem to be quasi-independent. Therefore, we divided the sectors into two parts from the data gap, where the telescope did not observe its targets. After the normalization, we initially realised that the amplitude of the eclipse light curve is similar in the first two sectors, whereas it changes in the sector 40 (see Fig. 1). Therefore, while we used all TESS data to derive mid-eclipse timings, we only made use of the light curves from sectors 14 and 15 for the light-curve modelling in Section 3 because they form a continuous light curve with more data points and excluded that acquired during the sector-40 from that analysis only.

To provide a wider wavelength coverage for the light curve analysis, as well as additional ground-based eclipse timings especially in the large gaps between space-borne observations, we observed the binary with (i) brand new 80-cm telescope (T80) at Ankara University Kreiken Observatory (AUKR), (ii) with the 1-m telescope (T100) in TUBITAK National Observatory (TUG), and (iii) with Istanbul University’s 60-cm telescope (IUO) at Çanakkale Onsekiz Mart University Ulupınar Observatory. These observations were dark, bias, and flat-field corrected and the differential photometry has been performed by using AstroImageJ software (Collins et al. 2017). The log of the observations and some related statistics are provided in Table 1. We selected the observations from AUKR on the dates of 2021 August 16, 17, 20, and 21 for the light-curve modelling due to their smaller scatter and full orbital phase coverage. The remaining observations were used only for the calculation of eclipse timings due to inferior photometric quality caused by unfortunate weather conditions or insufficient phase coverage.

Table 1.

Log of photometrical data of Kepler-451. The observations marked with asterisks are used in the light curve analysis. For Kepler and TESS, observation dates and photometric scatter (σphot) are for the phase folded and binned data, while the exposure times are the time between two binned data points.

Table 1.

Log of photometrical data of Kepler-451. The observations marked with asterisks are used in the light curve analysis. For Kepler and TESS, observation dates and photometric scatter (σphot) are for the phase folded and binned data, while the exposure times are the time between two binned data points.

Table 2.

Results from the light curve analysis of Kepler-451. The values marked with asterisks are fixed.

Table 2.

Results from the light curve analysis of Kepler-451. The values marked with asterisks are fixed.

Østensen et al. (2010) published radial velocity observations of the system. Their spectra is single lined, therefore, the radial velocity data are only for the primary component. We gathered this data set from CDS’s Vizier archive service2 and used it during the analysis, which we describe in the next section.

3 LIGHT-CURVE MODELLING

The light curves of the binary display a strong reflection effect around the secondary eclipse, due to the large surface temperature contrast and the small orbital separation between the components. While the primary eclipse is deeper than the secondary, the depths of both eclipses are smaller compared to the amplitude of the reflection, most probably due to relatively low orbital inclination. Although it is not easy to notice with the eye, the sdB component’s pulsations are detectable in Kepler light curves and have been reported by Østensen et al. (2010). Since we binned the Kepler and TESS light curves to 300 points each, we (i) reduced the computation time, and (ii) removed the effects of the pulsations of the sdB star. Therefore, the assumption that the light curves represent only the eclipses and phase-dependent modulations is safe.

Because of the radial velocity of the secondary is missing and the eclipses are not total, the mass ratio and related parameters are degenerate. Therefore, we adopted the mass value from Østensen et al. (2010), as 0.48 M⊙, which they derived by using both the spectral logg and the orbital solution, in the range for post-common envelope sdBs. Since, we adopted the mass of the primary, the mass ratio of the system can be calculated easily from the radial velocities as q = 0.25, and the mass of the secondary is 0.12 M⊙. Semi-major axis of the binary was also calculated as 0.891 R⊙.

We adopted and fixed the primary’s temperature to |$29564\, \mathrm{K}$| that Østensen et al. (2010) reported, which is comparable with the value calculated by Lei et al. (2018) from fitting the hydrogen and the helium lines with synthetic spectra. Other than the parameters mentioned above, we fixed the gravity brightening values (g1 = 1, g2 = 0.32) and the albedo of the primary (A1 = 1) to their theoretical values. We set all of the remaining parameters to freely vary. These parameters are, the temperature of the secondary (Teff, 2), radii of the components (R1 and R2), albedo of the secondary (A2), linear limb darkening coefficients, orbital inclination (i), and the passband luminosities. The binary has a very small separation and a short orbital period. Therefore, we assumed synchronous rotation for both of the components.

We used the second version of Phoebe (Prša et al. 2016) for the modelling of the light and radial velocity curves. We first used downhill simplex algorithm (Nelder & Mead 1965) to find an initial global solution that resembles the observations. We then used Levenberg–Marquardt method (LM) (Levenberg 1944; Marquardt 1963) for fine-tuning the parameter values. However, we failed to calculate uncertainties of the parameters, most probably due to the small number of iterations that otherwise would have taken considerably longer computation times. To calculate the uncertainties, as well as to confirm the results of the light curve modelling, we also used the first version of Phoebe (Prsa et al. 2011), by keeping the fixed and free parameters as they were in the modelling with the second version.

The best-fit parameter A2 is considerably different for the analyses with the two different versions of Phoebe. To investigate this difference, we made several fitting trials for different values of A2 with the first version of Phoebe. When we fixed A2 to unity, the value of Teff, 2 becomes 2987 K. On the other hand, the best-fit value of A2 is 0.1520(8) and Teff, 2 is 2829(1) K. Therefore, the difference of the best-fit values of A2 in the modelling attempts with different versions of Phoebe turned out to be due to high parameter degeneracy between Teff, 2 and A2.

Except A2, all the other best-fit parameter values from the two versions of Phoebe are almost the same (see Table 2). Our results have relatively small discrepancies with the observations that can be seen in Fig. 2, on the entire orbital cycle and also around the eclipses. Since our aim is to calculate the eclipse timings from the light curve model, but we could not improve them any better, we scaled the eclipse model to the observations that will be explained in the next section.

Figure 2.

Light curves of Kepler-451 along with synthetic models from Phoebev2 and their residuals. Light curves and the residuals are arbitrarily shifted and colour assigned.

4 ECLIPSE TIMINGS ANALYSIS

In order to determine the mid-eclipse times, we fitted our light-curve model cycle by cycle to each light curve by varying only the conjunction times while keeping all the other fit parameters unchanged. To compensate the small changes on the normalization level of each cycle, we added another parameter that shifts the brightness level of the model. We calculated mid-eclipse timings of our own observations, but we did not include them in the final ETV analysis since we found them not to be sufficiently precise for a meticulous analysis at the end. Therefore, our data set for the final ETV analysis consisted of only Kepler and TESS observations thanks to their superior photometric precision. The secondary eclipses are considerably shallow, which makes timing measurements imprecise, leading to a larger scatter on the ETV diagram. Therefore, we only used primary minima for the ETV analysis.

We also measured mid-eclipse timings from widely used KW method (Kwee & van Woerden 1956) for a comparison with the model-dependent mid-eclipse timings, which provided similar results. However, timings we derived from light-curve models scattered within a narrower range than that calculated with KW, which is valid for both Kepler and TESS data, as well as for primary and secondary minima. We plotted the measurement results of primary eclipses from Kepler with respect to each other (Fig. 3) for an illustration. As a result, we decided to base our timing analysis on the times of minima derived from light-curve model fitting. We provided the list of mid-eclipse timings derived from light-curve model fitting in Table A1.

Figure 3.

Comparison of primary mid-eclipse timings of Kepler data from widely used KW method (_x_-axis) and light curve model fitting method (_y_-axis). The red line shows the best-fit line to data, which has a slope of 0.91, whereas the black line shows the one-to-one relation.

For the Kepler subset of our data, the overall trend is similar to what Baran et al. (2015) showed, hinting a cyclic variation on the order of hundred days and a longer-term variation, that resembles a secular increase on the period. TESS subset, on the other hand, has more evidence for the long-term variation to be cyclic. The trend is a decreasing one during the sectors 14 and 15, but turns and changes to an increase during the sector 40 supported by the positions of the most recent timings from our own observations. This smooth change indicates that the long-term variation in the Kepler data has a periodic nature.

We then performed a frequency analysis to only primary mid-eclipse timings in the Kepler data set. Other than the orbital period frequency and its harmonics, we initially encountered a strong peak at ∼3.7 d and with an amplitude of ∼2.8 s in the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) of the ETV diagram. This variation can easily be noticed with the eye and produces a significant fraction of the scatter in the Kepler data. Assuming the reason for this variation is the light-time effect (LiTE), its source will be unseen third body with a minimum mass was of ∼0.044 M⊙ on an orbit with a minimum semi-major axis of a = 0.038 au. When we phase the ETV diagram with this period, data points clump at regular intervals, indicative of a beat phenomenon between sampling and periodicity. We found out that this clumping emerges in the secondary Kepler minima as well with a phase shift of 0.5, which is not observed in the TESS observations at all. This signal with 3.7 d period is in fact a spurious one caused by the sampling of Kepler data and the pulsations of the sdB component (Barlow, Wade & Liss 2012). Since it has no real physical origin, we ignore this signal totally in our attempt to model the mid-eclipse timing data.

4.1 Two long period circumbinary planets

We first assumed two LiTEs to model (i) the periodic variation that Baran et al. (2015) interpreted as a CB planet and (ii) the longer-term variation. The magnetic stellar winds are expected to remove mass from the sdB primary, hence decrease the orbital period of the binary. Therefore, we also assumed a model with a secular change superimposed on two LiTEs (see Section 4.2 in Esmer et al. 2021 for formulas and definitions of ETV models). We used our own python code for ETV modelling. We used Markov Chain Monte Carlo method for finding the global solution, by making use of the emcee code Foreman-Mackey et al. (2013). We used Gaussian likelihoods and priors and set 100 walkers and 2000 random walks in total. After finding the global minimum of our model by taking the mode of the distributions of each parameter, we set them as initial values for fine-tuning their values with down-hill simplex and LM algorithms.

The best-fit solution for the first LiTE component of our model is similar to what Baran et al. (2015) found (see Table 3). We also found that this CB planet with a minimum mass of 1.86 _M_jup should have an eccentric orbit (e = 0.33) although the period value we found (406 d) is slightly smaller than theirs (416 d). The longer period (P = 1460 d) LiTE, on the other hand, corresponds to another circumbinary planet with a similar minimum mass (1.61 _M_jup), on a wide (a sin i = 2.1 au), and eccentric orbit (e = 0.29) (see Fig. 4). The uncertainty on the long-term LiTE is large, likely due to insufficient data coverage in one cycle, and heteroscedasticity of the data set, which is formed by a combination of Kepler and TESS observations.

Figure 4.

Mid-eclipse timings of Kepler-451 along with quadratic and two planet model (solid red). Each model component shown with different colours and line styles. Our observations from TUG, AUKR, and IUO are shown in magenta. Note that the residuals divided in two regions and shown with different scales.

Theroot mean square (rms) of the residuals of the two planet solution is 4.85 s. The 3.7 d periodic signal is ignored as explained above, however, it affects the goodness of fit statistics of our best-fit solution. The |$\chi _\nu ^2$| of the best-fit is 12.625 for a degrees of freedom of 8970, which is usually interpreted as either a poor fit or a result of underestimated observational errors. If we remove 3.7 d signal by fitting a sine function, the |$\chi _\nu ^2$| decreases to 9.332, which is still not statistically satisfactory. This is, again, likely due to either poor representation of the sine function or underestimated timing uncertainties. We leave this discussion here and adopt our best-fit two planet solution.

4.2 Other potential explanations of the ETV

To explain whether the ETVs of Kepler-451 is caused by the magnetic activity of the cooler component we used the online tool3 from Völschow et al. (2016) that calculates the energy required to produce the quadrupole moment based on three different models, (i) finite-shell two zone model, (ii) finite-shell constant density model (Völschow et al. 2016), and (iii) thin shell model (Tian, Xiang & Tao 2009). For all three models, the ratio of the required energy (Δ_E_) to the available energy for the magnetically active star (E_⋆) has to satisfy the condition Δ_E/_E_⋆ < <1. By using the value of the coefficients _k_1 and k_2 (see the formulation in Völschow et al. 2016) as 0.133 and 3.42, we calculated the energy required for the cooler component to derive Applegate (1992) mechanism based on the three models mentioned above. Considering the 1460 d periodic modulation, we calculated the minimum value of Δ_E/_E_⋆ for finite-shell two zone model as ∼12, for finite-shell constant density model as 103 and for the thin shell model as ∼0.7. We also calculated the same parameter for the case of 406 d modulation as 143 for finite-shell two zone model, 103 for finite-shell constant density model and 8 for the thin shell model. Since the magnetic activity of the dM component should have weakened significantly over the time it takes for the progenitor of the sdB component to have evolved from the main-sequence, Applegate mechanism is not expected to be effective in causing quadruple moment changes that lead to ETVs. In addition, magnetic-activity induced ETVs have more of a cyclic nature rather than a periodic one, with changing amplitudes from one cycle to another. Color, brightness, temperature, and maxima-level differences expected from strong magnetic activity variations are not observed in the light curves of the system either. Therefore, we concluded that the cyclic ETVs of Kepler-451 can not solely be explained by magnetic activity modulation. On the other hand, the orbit of the binary is circular, and the secondary minima are not out-of-phase with the primaries, hence apsidal motion cannot be a possible explanation to the ETV of the system.

4.3 The third planet on a smaller orbit

We then continued and searched for signs of periodicity in the residuals of the two planet solution in order to check if we are under-fitting since the |$\chi _\nu ^2$| is still large. For this purpose, we used our own script that makes use of the astropy (Astropy Collaboration 2013, 2018) module’s LS periodogram tool (VanderPlas 2018). We ignored the binary’s orbital frequency and the 3.7 d signal by limiting the range of frequency spectrum to a maximum of 0.2 d−1.

We first analysed the primary mid-eclipse timings of Kepler data. The peak frequency in the LS power spectrum (Fig. 5) corresponds to a periodicity at 43.0 d. This peak has a False Alarm Probability (FAP) of 4.91 × 10−15, and an amplitude of 0.554 s. The same peak exists in the Kepler secondary timings with almost the same amplitude, however, FAP of the peak for the secondaries is unity. Next, we searched for the same frequency in the TESS subset. The total time span of sector 14 and 15 is 54 d. Sector 40 is separated from the other two sectors by around 1.8 yr. So, we searched for potential frequencies only in sectors 14 and 15. The frequency with the highest power corresponds to ∼39 d then, which is close to what we have found from the Kepler subset. However, due to relatively larger scatter and timing uncertainties, the FAP value is almost unity. The amplitude, on the other hand is around 20 times larger, again, compared to the Kepler subset. We think that the highest peak in TESS spectrum is likely to be not related to 43 d signal that we found in frequency spectrum of the Kepler primary minima.

Figure 5.

Lomb–Scargle spectrum of the residuals of mid-eclipse timing data from the fitting of two planet solution. The highest peak corresponds to 43-d periodic signal with almost zero FAP.

We investigated if there is an improvement in the ETV fit when we add the 43-d frequency to our model as an additional LiTE. We set the initial values of the period to 43 d and the semi-amplitude (ALiTE, 3) to 0.5 s. Since the amplitude is too low to model eccentricity, considering timing uncertainties, we assumed a circular orbit for this putative third planet. We also kept every other parameter related to the secular change and the previous two planet solution to their best-fit values. The only free parameters of this model were the period, the amplitude and the mean anomaly for the third LiTE, the best-fit values of which, corresponds to a minimum mass of 1.76 M|$\rm _{jup}$| orbiting the binary at a minimum separation of 0.2 au. The |$\chi _\nu ^2$| value reduced from 12.625 to 12.502 with the addition of the third LiTE

We phase folded the Kepler light curves to 5, 10, 20, and 50 cycles and recalculated the mid-eclipse times with light curve fitting as we did for the original timings. We confirm that the 43-d signal is present on all of the binned data sets. We plotted only the 50-cycle binned mid-eclipse timings in Fig. 6 as well as binned original Kepler timings with the same number of data points to make the third planet signal for better illustration. Finally, we plotted the binned data sets phased to the third planet period in Fig. 7.

Figure 6.

Three planet solution of the mid-eclipse modelling of Kepler-451 in a close-up view to the Kepler data, excluding Quarter 0. Black circles represent the mid-eclipse timings from light-curve model fitting to 50-times binned light curves. Blue circles represent the binned O-C data. Note that both binned data follows the trend of the 43-d modulation.

Figure 7.

Phased O-C signal corresponding to the third planet model with the period of 43 d and its residuals. Only the binned data set plotted for clarity.

The final search for a periodicity in the residuals of this three-planet model; however, resulted in no-detection with a significant FAP value. Therefore, our final ETV solution consists of three similar-mass Jovian planets orbiting the binary, with at least two of them on eccentric orbits (Table 3). We investigate if such a system will be dynamically stable in the next section.

Table 3.

Results from the eclipse timing variation analysis of Kepler-451. The values in parentheses are uncertainties on the last digits.

Notes. *Eccentricity of the inner planet is assumed to be zero. **The time and the argument of the periastron of the inner planet are correlated and arbitrary. n denotes the planet number from inside out while the planet names are given in parentheses.

Table 3.

Results from the eclipse timing variation analysis of Kepler-451. The values in parentheses are uncertainties on the last digits.

Notes. *Eccentricity of the inner planet is assumed to be zero. **The time and the argument of the periastron of the inner planet are correlated and arbitrary. n denotes the planet number from inside out while the planet names are given in parentheses.

5 DYNAMICAL STABILITY OF THREE PLANET MODEL

We made an orbital simulation that corresponds to the best-fit parameters to check if the system keeps the initial orbital configuration for a long duration. We assumed all the objects in the system to be co-planar. We also assumed the binary as a single object with a total mass of 0.6 M⊙, which is a very good assumption considering the significantly shorter orbital period compared to that of even the shortest orbital-period planet. We used the Rebound code (Rein & Liu 2012) with WHFAST integrator (Rein & Tamayo 2015) for the simulation. Total time of the simulation was set as 10 Myr with a time-step of 0.5 d. In this test, the orbital coordinates of the objects chaotically oscillate, therefore, the system is unstable.

In order to investigate if a stable solution exists in a parameter space close to the best-fit solution, we performed a Frequency Map Analysis (Laskar 1990, 1993), similarly to our earlier work in Esmer et al. (2021). The timing data cover fewer orbital cycles of the outermost object, hence the uncertainty on its parameters are larger compared to the others. Therefore, we varied the outermost object’s semi-major axis (a) and eccentricity (e) around the best-fit values, while keeping the parameters of the remaining objects to their best-fit values.

We set the total integration time to 106 d, which is around ∼680 orbital revolution of the outermost object, and we set the time-step as 0.5 d. The semi-major axis was changed in the range of 1.6–2.6 au, while the eccentricity was changed between 0 and 0.5. For each integration, we calculated the change of the frequency of the mean motion in two halves of the total integration time, and calculated a normalized stability index, D. In the case of stability the value of D has to satisfy the condition, |$\rm \log D \lt -6$| (Correia et al. 2005). For the calculation of the frequencies of the mean motion, we used the TRIP code (Gastineau & Laskar 2011).

We initially searched the a–e ranges mentioned above in 60 × 60 resolution. Then, we identified potentially stable sub-ranges and performed the same search for them to reduce the computational time needed for integrations. The three black rectangles shown in Fig. 8 correspond to the potentially stable regions. From left to right, the first region around a = 1.95 au and e = 0.18 has very narrow and scattered a–e coordinates that have normalized stability indices which are closer to the stability condition, none of which actually satisfies it. We think that outermost object is likely to be unstable because it is found to be closer than the best-fit orbital separation to the remaining objects. Therefore, this first candidate region does not represent a stable solution.

Figure 8.

Stability map of the three planet solution of Kepler-451. The central black circle corresponds to the best-fit ETV solution, while the error bars corresponds to 1σ uncertainty ranges. Colour bar indicates the logarithm of normalized stability index, and the stability occurs when its value is < −6 (black regions). Three rectangles are potentially stable regions. See the text for more details.

The second region around a = 2.2 au and e = 0.15, corresponds to an orbital separation for the outermost object closest to the best-fit solution. This darker and narrow region has very few a–e coordinates that actually satisfies the stability condition.

The third region that corresponds to a = 2.5 au has a wide range of stable configurations. The allowed eccentricity values are in a range of 0 – 0.4, for which the upper limit even exceeds the best-fit value for this parameter. The separations of the stable region, on the other hand, are larger around almost 0.3 au (1.5 σ) from the best-fit value.

We would like to remind that, our eclipse timing solution does not take into account the mutual interactions. One should also keep in mind that the fixed parameters of the outermost object as well as the parameters of the other objects have uncertainties that we did not take into account for the sake of computation time. Our assumption of co-planar orbits is likely to be a good approximation, but there may be slight differences between the mutual orbital inclinations of the objects. Despite all of these, the dynamical stability test shows that stable orbital configurations close to our three planet best-fit eclipse timing solution exist.

6 CONCLUSIONS

Planets around evolved stellar systems are interesting because they may have formed in the accretion discs made up of the ejected material from the stars at some stage of their evolution. Kepler-451 is one such system, with an evolved sdB primary and a dM secondary. We detect two more additional periodic signals in the ETV that we interpreted as planets in the system in addition to the confirmation of the Kepler-451 b (Baran et al. 2015) based on our analysis of the ETVs of the central binary. These two planets have similar minimum masses (m|$\rm _c \sin i$| = 1.61 ± 0.14 _M_jup and m|$\rm _d \sin i$| = 1.76 ± 0.18 _M_jup) to that of the already found planet-b (m|$\rm _b \sin i$| = 1.86 ± 0.05 _M_jup). Since the sdB primary of Kepler-451 system must have ejected a significant fraction of its mass, these three planetary-mass bodies might be second-generation planets formed in this material.

The projected orbital separations are a|$\rm _d \sin i$| = 0.20 ± 0.03, a|$\rm _b \sin i$| = 0.90 ± 0.04, and |$a\rm _c \sin i = 2.10 \pm 0.02$| au, respectively. These three similar-mass, Jovian planets of Kepler-451 should be stable on their orbits even they were eccentric, if the outermost planet’s orbital separation is ∼0.3 au larger than our best-fit Keplerian solution, which requires the outermost planet’s orbital period to be at least ∼1800 d. Nevertheless, slight differences in their mutual orbital inclinations and/or eccentricities may lead to other stable configurations too.

The total radial velocity contribution of the three planets then should have a semi-amplitude of ∼250 m/s depending on their instantaneous mean anomaly, which can be traced with radial velocity observations having achievable precision but on long time-scales for verification with an independent method. On the other hand, the orbital inclination of the binary is 70°. It is very likely for the circumbinary planets to have similar orbital inclinations. Therefore, they should not show any transits at all, and any information of their radii will be absent until another method is suggested. However, if they have formed from the ejected outer layers of the sdB primary, they should be expected to have very low densities. Therefore, finding planetary-mass bodies in transit of their sdB hosts would be very interesting.

The light curve modelling of Kepler-451 showed that the reflection effect on such sdB+dM systems is not trivial to handle. The Kepler light curve and their correlated residuals are not symmetrical when we fold the light curve with respect to the second eclipse, implying that the brightness distribution due to reflection is not homogeneous on the secondaries reflected hemisphere. This can be due to poor-handling of heat re-distribution over the stellar surfaces by the models. However, the effect is small and should have no influence on the determination of the mid-eclipse times.

Calculations of the mid-eclipse timings from the light-curve models provide much more reliable data sets for an ETV analysis, as is the usual practice in TTV studies based on model-dependent transit timings. Model-based timings are especially handy in the presence of phase-dependent brightness variations such as spot or pulsation modulations, which pose major problems in traditional methods such as KW due to the light curve asymmetries they cause.

We would like to mention that the mid-eclipse timings calculated by Krzesinski et al. (2020) from the Kepler light curves are significantly different from the values we derived. We calculated mid-eclipse timings with KW and light curve model fitting from Kepler PDC-SAP light curves and found similar results that only differ in their scatter. However, our mid-eclipse timings do agree with that Baran et al. (2015) derived from the same Kepler light curves, the analysis of which ended up in sinusoidals with very similar parameters that we also attributed to the first CB planet (Kepler-451 b) in the system suggested by them.

Together with the suggested two additional planets in Kepler-451 system, the number of CB planets discovered with the timing technique increases up to 22 when the NASA Exoplanet Archive and exoplanet.eu data bases are considered together. Analysing the eclipse timing variations will continue to play a major role not only in revealing additional CB bodies on larger orbits which allow us to probe a different parameter space, but also study current and past astrophysical phenomena in the evolutions of these systems.

SUPPORTING INFORMATION

lightcurves.dat

timings.dat

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We gratefully acknowledge the support by The Scientific and Technological Research Council of Turkey (TÜBİTAK) with the project 118F042. We thank TÜBİTAK for the partial support in using T100 telescope with the project number 19AT100-1471 and 17BT100-1208. This work is also supported by the research fund of Ankara University (BAP) through the project 18A0759001. EME acknowledges support from TÜBİTAK (2214-A, number 1059B141800521). IST60 telescope and its equipment are supported by the Republic of Turkey Ministry of Development (2016K12137) and Istanbul University with the project numbers BAP-3685, FBG-2017-23943. Some/all of the data presented in this paper were obtained from the Multimission Archive at the Space Telescope Science Institute (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia. We thank the developers of the phoebe code, especially Dr. Kyle Conroy from Villanova University for his help in using the code.

DATA AVAILABILITY

Eclipse timing measurements and the light curves appearing for the first time in this work are presented as online material in the machine readable format and can be found in the CDS linked to this study. Although we provide the light curves within the online data; we can also provide raw data sets of our own observations on request from the corresponding author.

Footnotes

REFERENCES

Applegate

J. H.

,

1992

,

ApJ

,

385

,

621

Astropy Collaboration

,

2013

,

A&A

,

558

,

A33

Astropy Collaboration

,

2018

,

AJ

,

156

,

123

Baran

A. S.

,

Silvotti

R.

,

2020

,

Research Notes of the AAS

,

4

,

237

Baran

A. S.

,

Zola

S.

,

Blokesz

A.

,

Østensen

R. H.

,

Silvotti

R.

,

2015

,

A&A

,

577

,

A146

Barlow

B. N.

,

Wade

R. A.

,

Liss

S. E.

,

2012

,

ApJ

,

753

,

101

Bear

E.

,

Soker

N.

,

2014

,

MNRAS

,

444

,

1698

Beuermann

K.

et al. ,

2010

,

A&A

,

521

,

L60

Chen

W.-C.

,

Podsiadlowski

P.

,

2017

,

ApJ

,

837

,

L19

Collins

K. A.

,

Kielkopf

J. F.

,

Stassun

K. G.

,

Hessman

F. V.

,

2017

,

AJ

,

153

,

77

Correia

A. C. M.

,

Udry

S.

,

Mayor

M.

,

Laskar

J.

,

Naef

D.

,

Pepe

F.

,

Queloz

D.

,

Santos

N. C.

,

2005

,

A&A

,

440

,

751

Esmer

E. M.

,

Baştürk

Ö.

,

Hinse

T. C.

,

Selam

S. O.

,

Correia

A. C. M.

,

2021

,

A&A

,

648

,

A85

Foreman-Mackey

D.

,

Hogg

D. W.

,

Lang

D.

,

Goodman

J.

,

2013

,

PASP

,

125

,

306

Gastineau

M.

,

Laskar

J.

,

2011

,

ACM Commun. Comput. Algebra

,

44

,

194

Hardy

A.

et al. ,

2016

,

MNRAS

,

459

,

4518

Krzesinski

J.

,

Blokesz

A.

,

Siwak

M.

,

Stachowski

G.

,

2020

,

A&A

,

642

,

A105

Kwee

K. K.

,

van Woerden

H.

,

1956

,

Bull. Astron. Inst. Neth.

,

12

,

327

Laskar

J.

,

1990

,

Icarus

,

88

,

266

Laskar

J.

,

1993

,

Physica D: Nonlinear Phenomena

,

67

,

257

Lee

J. W.

,

Kim

S.-L.

,

Kim

C.-H.

,

Koch

R. H.

,

Lee

C.-U.

,

Kim

H.-I.

,

Park

J.-H.

,

2009

,

AJ

,

137

,

3181

Lee

J. W.

,

Hinse

T. C.

,

Youn

J.-H.

,

Han

W.

,

2014

,

MNRAS

,

445

,

2331

Lei

Z.

,

Zhao

J.

,

Németh

P.

,

Zhao

G.

,

2018

,

ApJ

,

868

,

70

Levenberg

K.

,

1944

,

Quarterly of Applied Mathematics

,

II

,

164

Lomb

N. R.

,

1976

,

Ap&SS

,

39

,

447

Marquardt

D. W.

,

1963

,

SIAM J. Appl. Math.

,

11

,

431

Nelder

J. A.

,

Mead

R.

,

1965

,

Comput. J.

,

7

,

308

Orosz

J. A.

et al. ,

2012

,

Science

,

337

,

1511

Østensen

R. H.

et al. ,

2010

,

MNRAS

,

408

,

L51

Prsa

A.

,

Matijevic

G.

,

Latkovic

O.

,

Vilardell

F.

,

Wils

P.

,

2011

,

Astrophysics Source Code Library, record ascl:1106.002

.

Prša

A.

et al. ,

2016

,

ApJS

,

227

,

29

Rein

H.

,

Liu

S. F.

,

2012

,

A&A

,

537

,

A128

Rein

H.

,

Tamayo

D.

,

2015

,

MNRAS

,

452

,

376

Scargle

J. D.

,

1982

,

ApJ

,

263

,

835

Tian

Y. P.

,

Xiang

F. Y.

,

Tao

X.

,

2009

,

Ap&SS

,

319

,

119

VanderPlas

J. T.

,

2018

,

ApJS

,

236

,

16

Völschow

M.

,

Schleicher

D. R. G.

,

Perdelwitz

V.

,

Banerjee

R.

,

2016

,

A&A

,

587

,

A34

Zorotovic

M.

,

Schreiber

M. R.

,

2013

,

A&A

,

549

,

A95

APPENDIX A: MID-ECLIPSE TIMES OF KEPLER-451

Table A1.

Mid-eclipse times of Kepler-451 from our observations.

Table A1.

Mid-eclipse times of Kepler-451 from our observations.

© 2022 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society

Supplementary data

Citations

Views

Altmetric

Metrics

Total Views 924

530 Pageviews

394 PDF Downloads

Since 2/1/2022

Citations

9 Web of Science

×

Email alerts

Citing articles via

• Latest

• Most Read

• Most Cited

More from Oxford Academic