© R. Heller et al. 2020

1 Introduction

The star Kepler-160 (KIC 7269974, KOI-456), at a distance of light years (Gaia Collaboration 2016, 2018; Andrae et al. 2018), was almost continuously observed during the Kepler primary mission from 2009 to 2013. It is a Sun-like star with a radius of solar radii (R_⊙), an effective temperature of K, and a resulting Sun-like luminosity of times the solar value¹. The star has been identified as a host of two transiting-planet candidates in very close orbits (Borucki et al. 2011) that were later statistically validated (Rowe et al. 2014) as planets Kepler-160 b and c. Kepler-160 b has a radius of 1.7 Earth radii (R_⊕) and is in a ~4.3 d orbit, while Kepler-160 c, with a radius of about 3.1 R_⊕, orbits the star with a period of ~13.7 d. Kepler-160 c has been shown to exhibit transit-timing variations (TTVs) with a period of about 879 d (Holczer et al. 2016; Ofir et al. 2018), whereas Kepler-160 b does not seem to exhibit TTVs.

As part of our transit least-squares survey (Heller et al. 2019a,b), we analyzed the light curve of Kepler-160 to search for additional transiting planets that have previously been missed. We also examined the TTVs of Kepler-160 cwith orbital simulations of the Kepler-160 planetary system in combination with a Markov chain Monte Carlo (MCMC) fitting procedure to the Kepler transit photometry to search for additional planets. Here we report the detection of a new periodic transit feature in the light curve that might originate from a 1.9 R_⊕-sized planet receiving virtually Earth-like insolation from its Sun-like host star. We also show that the TTVs of Kepler-160 c are neither caused by this new transit candidate nor by the previously known planet Kepler-160 b and that there must be yet another nontransiting planet in this system.

In Sect. 2 we explain the methods we used for the light-curve detrending and transit search, for the MCMC-based system characterization, for the automated false-positive vetting, and for the TTVs study through N-body simulations. In Sect. 3 we present our results, Sect. 4 is dedicated to a discussion of our findings, and in Sect. 5 we summarize and conclude.

2 Methods

2.1 Light-curve detrending and transit search

We used the publicly available pre-search data conditioning simple aperture photometry (PDCSAP) flux of the star Kepler-160 of the Kepler Data Release 25 (DR25) for Kepler Quarters 1–17 (Thompson et al. 2016a; Twicken et al. 2016) that is accessible through the Mikulski Archive for Space Telescopes (MAST)². This light curve contains about 60 000 brightness measurements of the star between 13 May 2009 and 11 May 2013, mostly with a cadence of 30 min. The time frame of the data is the Barycentric Kepler Julian Day (BKJD), which is equal to the Barycentric Julian Date (BJD) subtracted by 2 454 833.0 d. BKJD = 0 corresponds to January 1, 2009, 12:00:00 UTC.

In our search for new transiting planets around Kepler-160, we first masked any data from the light curve that were affected by transits of the previously known planets Kepler-160 b and c. We removed any variability on timescales longer than the expected transit duration using the Wōtan software package³ (Hippke et al. 2019). Wōtan offers more than a dozen detrending filters, and we decided to use Tukey’s biweight filter (Mosteller & Tukey 1977), which has been shown to yield the highest recovery rates of injected transits from simulated data (Hippke et al. 2019). We chose to use a running window of three times the duration of a central transit on a circular orbit at any given trial orbital period to detrend the light curve from long-term stellar variability and from any systematic trends on timescales longer than this value. This width has been shown to yield the best compromise between the trend removal and a high recovery rate of synthetic transits from simulated noisy light curves (Hippke et al. 2019).

Then we passed the resulting detrended light curve to the transit least-squares (TLS) algorithm⁴ (Hippke & Heller 2019). We used an oversampling factor of five to increase our sensitivity at the cost of computational effort, and we used estimates of the stellar limb-darkening coefficients, mass, and radius from the revised Kepler Input Catalog (Brown et al. 2011; Mathur et al. 2017) to adapt the TLS transit search function to Kepler-160. The key detection criteria for our search were the signal-to-noise ratio (S/N) and the signal detection efficiency of TLS (SDE_TLS; Hippke & Heller 2019). The S/N was calculated as (δ∕σ_o) n^1∕2 with δ as the mean depth of the trial transit, σ_o as the standard deviation of the out-of-transit points, and n as the number of in-transit data points (Pont et al. 2006). For a detection to be significant, we chose a minimum S/N of 7.1, which has been shown⁵ to yield one false positive for a threshold-crossing event over the entire Kepler mission (Jenkins et al. 2002). For the SDE_TLS value we applied a minimum threshold of 9, which results in a false-positive rate < 10⁻⁴ in the limiting case of white noise (Hippke & Heller 2019).

We also inspected the light curve for evidence of an eclipsing binary. First, our automatic search with TLS only accepted candidate transit sequences for which the average depths of the odd and even transits agreed within 3 σ. In the case of the candidate transit signal that we report here, we had two odd transits and a single even transit. Second, we ensured that there was no evidence of a secondary eclipse at the > 3 σ level compared to the local noise at half an orbital phase after the candidate transit. Third, to relax the circular-orbit assumption, we performed a TLS search for secondary eclipses with orbital phases ranging between 0.1 and 0.9 (the transit being at 0). Although any secondary eclipse should have the same period (P) as the transit signal that we found, we allowed for a small variation of P in our TLS search. As a result, the strongest signal always had both SDE_TLS < 3 and S∕N < 3, and so we did not find any signal of a secondary eclipse with statistical significance.

We also verified that the period of the transit signal differs sufficiently from any periodic modulation of the light curve due to stellar and possibly systematic variability. As a tool for inspection, we chose the Lomb-Scargle power spectrum (Lomb 1976; Scargle 1982) algorithm for unevenly spaced data (Press & Rybicki 1989) as implemented in the astropy software package⁶.

All of these automated search metrics, together with the Kepler PDCSAP flux, the Wōtan detrending curve, the phase-folded light curve, secondary eclipse search, etc. were displayed on a vetting sheet (for an example, see Heller et al. 2019a). Our vetting sheet for Kepler-160 is one of hundreds that we inspected in our TLS survey. Some of the aspects that typically prompted us to reject a candidate based on the vetting sheet after it had already passed the automated detection criteria stated above were the following: (1) the transit period was very close to a period of activity in the Lomb-Scargle power spectrum, (2) there was substantial periodic activity in the light curve, (3) there were long gaps in the data, (4) the transit depths were highly variable, (5) there were features in the phase-folded light curve with amplitudes comparable to the primary signal, and (6) the transit shape was not convincingly “planet-like”. Many of these aspects are highly subjective, but we nevertheless considered them very helpful in rejecting a large part of the automated candidate detections. Kepler-160 passed all of these human-based interventions in addition to our automated search criteria.

2.2 MCMC transit fitting of the light curve

Our post-detection characterization of this putative three-planet transiting system involved the MCMC sampler emcee (Foreman-Mackey et al. 2013). As inputs, we provided the mid-times of the first transits and P for all three planets. The free parameters of our three-planet transit model were P, the transit epoch (T₀), the planet-to-star radius ratio (R_p∕R_s), and the transit impact parameter (b) for each planet as well as the stellar density (ρ_s) and two limb-darkening coefficients for a quadratic limb-darkening law (Kipping 2013) that we used as global parameters for all transits. Our MCMC analysis was initialized with 100 walkers, and each walker performed 300 000 steps. The first half of each walk was discarded to ensure that we preserved only burned-in MCMC chains. We visually checked that the Markov chains were sufficiently burnt-in. We selected the burn-in time by inspecting the chains for trends and changes in the spread, and we selected the acceptance fraction using a trace plot of the chains (Roy 2020).

We computed the planetary radii from the posterior distribution of R_p ∕R_s and from the estimate of R_s from the Gaia DR2 (Gaia Collaboration 2016, 2018), which was derived from a combination of Gaia photometry and parallax measurements (Andrae et al. 2018). For the sake of a self-consistent data set, we used K from the Gaia DR2 to estimate the stellar luminosity (see Sect. 4.3 for a discussion of T_eff). From the posterior distributions of ρ_s and the orbital parameters and using the Gaia DR2 values for R_s and T_eff as priors, we calculated the top-of-the atmosphere stellar flux received by the new transiting planet candidate.

Fig. 1 Panel a: Kepler CCD photometry of Kepler-160 from the Data Validation Report. The dotted line outlines the aperture mask used to extract the light curve. The cross denotes the target position, and the stars (with annotated KIC numbers) indicate the positions of stars that might act as contaminants of the target light curve. The color scale illustrates the difference between mean flux out-of-transit and in-transit of Kepler-160 b (KOI-456.01) after normalizing by the uncertainty in the difference for each pixel. North and east are indicated with yellow lines. The Kepler CCD pixels are squares and appear stretched in this figure from the Data Validation Report. Panel b: Pan-STARRS image that shows Kepler-160 in the center. The imageis properly oriented so that north coincides with the Kepler image in panel a. The distance of 6.3′′ to one nearby star is indicated with an arrow.

2.3 False-positive vetting

We inspected the Kepler CCD pixel maps from the Data Validation (DV) Report of Kepler-160 (available at the NASA Exoplanet Archive) to search for possible stellar contaminants of the light curve. Stars that are physically unrelated but apparently close to the target on the plane of the sky can contribute to the total flux received within the aperture selected on the Kepler CCDs. If these stars are, for example, members of an eclipsing background binary or triple system, or if the target star is accompanied by a grazing eclipsing stellar binary, then the resulting flux measurement in the Kepler photometry can lead to a false-positive detection of a planet-like transit signal.

Figure 1a shows the Kepler CCD image of Kepler-160 from Kepler Quarter 5 (Q5), during which the first transit of our new exoplanet candidate was observed. The color scale illustrates the S/N for the previously known planet Kepler-160 b. The position of the target star is denoted with a cross and additional stellar sources are marked with stars. Two stars north of the target (i.e., to the lower right in that image) are within the aperture mask used for the extraction of the light curve, which is outlined by a dotted line. Three more stars south-east of the target (left on the image) are about 1.2 pixels away from the edge of the aperture mask that was used to extract the light curve of Kepler-160. Figure 1b shows a high-resolution image of the region around Kepler-160 from Pan-STARRS, which we used as an additional tool to assess the possibility of stellar contamination.

We used the lightkurve software (Lightkurve Collaboration 2018)⁷ to create per-pixel S/N charts for the new transit signal from differential CCD images and to validate that the newly discovered transit signal comes from the Kepler-160 target pixel and not from any of these nearby sources (see Sect. 3.1.1). Briefly, the per-pixel S/N of the difference images is the difference between in-transit and out-of-transit flux in each pixel divided by the propagated uncertainty in each pixel (Bryson et al. 2013).

Because the edge length of a Kepler CCD pixel corresponds to 3.98′′ on the sky, we chose a value of ρ = 3.98′′ for the maximum aperture radius interior to which the signal must have been produced. This is a conservative estimate compared to the value of 0.5′′ used by Morton et al. (2016) to evaluate the false-positive probabilities (FPPs) of Kepler-160 b and c. We used the vespa software⁸ (Morton 2012, 2015) to evaluate the FPP of the signal, for example, the probability that it was caused by a blend or stellar eclipsing binary or triple system. For the vespa maxrad value we adopted a value of 3.98′′, and we supplied vespa with a phase-folded light curve of our new transit candidate from which all the transits of the previously known planets Kepler-160 b and c had been removed.

In addition to the Pan-STARRS optical images of Kepler-160 and its celestial neighborhood, we also inspected high-resolution adaptive optics images⁹ of Kepler-160 obtained by the Robotic Laser-Adaptive-Optics imaging survey (Law et al. 2014). These images show no evidence of contaminants with a magnitude difference <3 within about 4′′ around Kepler-160. Although we cannot use this information for our analysis, which requires the input of the maximum possible radius for contamination, this information further increases our confidence that the transit signature that we discovered is not a false positive caused by contamination from an eclipsing or grazing binary or an unrelated transiting planet in the background.

To assess the possibility that our new candidate might be caused by instrumental effects, such as the rolling-bandeffect or a combination of rolling-band and statistical artifacts, we vetted it using the Model-Shift software¹⁰ (Coughlin 2017).

2.4 Transit-timing analyses of Kepler-160 c

We analyzed the TTVs of Kepler-160 c with a dynamical model to determine the origin of the TTVs and to constrain the physical parameters of the perturber. Our model is similar to the photodynamical model applied by Freudenthal et al. (2018, 2019), except for the fact that it calculates the transit mid-times instead of the light curve. We used the N-body code Mercury (Chambers 1999) to simulate the three-body (two planets plus one star) orbital dynamics of the Kepler-160 system over the duration of the observations. From the resulting orbital positions we calculated the transit times. The out-of-transit dynamics were numerically integrated with a typical time step of a twentieth of the period of the innermost planet and a time step of 0.05 d near the transits. In a first attempt, we investigated whether Kepler-160 b might be the cause of the observed TTVs of planet c. The results of this initial study were negative, and so we focused on additional planets in the system.

Then we explored the parameter space of our model to identify the regions that maximize the likelihood of the set of parameters given the data. For this purpose, we again used emcee (Foreman-Mackey et al. 2013) to deduce the posterior distributions of the planetary parameters. Our focus was on the orbital period and mass of the unknown perturber of Kepler-160 c. In order to generate the simulated TTVs of Kepler-160 c from our Mercury simulations, we used the following jump parameters for the MCMC fitting: the semimajor axis of planet c as well as the eccentricities, arguments of periastron, and the true anomalies of planets c and the test planet. The inclinations of both planets were fixed to 90° (see Sect. 4.4 for a discussion).

Because the TTV signal of Kepler-160 c is periodic (see Sect. 3.2), we studied the TTVs caused by test planets nearvarious orbital mean motion resonances (MMRs) with low integer ratios. All jump parameters had uniform priors with broad boundaries in order to avoid unphysical results. Starting values for the period of the perturber were calculated from the super-period given by the periodicity in the TTVs of Kepler-160 c and its own period (Eqs. (5)–(7) in Lithwick et al. 2012). The starting values of the remaining jump parameters were chosen manually sothat the model approximately matched the data and were then optimized according to the MCMC procedure.

3 Results

3.1 Transit detection

Figure 2 shows the results of our TLS search after masking out the known transits of Kepler-160 b and c. Panel a illustrates the Lomb-Scargle power spectrum of the PDCSAP Kepler light curve prior to our own detrending, and panel b shows the spectrum of the detrended light curve. A comparison between panels a and b shows the effect of the removal of stellar activity by Wōtan at frequencies ≲ 0.1 d⁻¹, or equivalent at periods ≳ 10 d. Figure 2c shows the SDE_TLS spectrum with a notable peak at a period of 378 d. Any periodic activity in the Lomb-Scargle periodograms of panels a and b is separated by at least one order of magnitude in period from the 378 d period of the suspicious transit signal. As a consequence, we do not find any clear signs of a connection between stellar activity and this candidate signal.

We also note a spurious peak at a period of 6.54206 d, or a frequency of 0.15286 d⁻¹ in the SDE_TLS spectrum of Fig. 2c. This value is not near any integer multiple of the period of our candidate signal. We have also verified thatthe spurious signal is not an alias of the new transit candidate signal (and vice versa). Moreover, the spurious signal remained in the SDE spectrum after the three transits of our new candidate were cut out from the light curve. We also split the Kepler light curve into two halves, which showed that the spurious signal is significantly more enhanced in the first half. We conclude that the spurious signal and the new transit candidate signal have independent origins. Instead, the spurious SDE peak at 0.15286 d⁻¹ in Fig. 2c coincides with a moderate peak in the LS power spectrum in Fig. 2a, which disappears after detrending in Fig. 2b. We interpret this as evidence that the 0.15286 d⁻¹ signal is caused by stellar activity. We also inspected the phase-folded light curve of this spurious signal, and it looks only marginally consistent with a transit shape. Finally, its S/N of 5.6 (compared to 9.5 for our new candidate) also suggests that this signal might be caused by alternative noise sources rather than by a transit.

Our detection of the candidate signal with TLS gave us first estimates of P, T₀, the transit depth (δ), and transit duration (d) of the candidate. The signal is related to a sequence of three periodic transit-like events distributed over Q5, Q9, and Q13. Our initial measurements of the individual transits depths were obtained with TLS from the mean in-transit flux (Heller 2019; Hippke & Heller 2019). The measured depths for the triple transit sequence are 207 (± 54) ppm, 237 (± 46) ppm, and 370 (± 41) ppm, respectively. The mean value is 271 ppm, and the individual deviations are − 1.2 σ, − 0.8 σ, and + 2.4 σ, respectively. We also searched for further periodic transit sequences beyond those of Kepler-160 b, c, and the new candidate with TLS by masking the transits of the two planets and of the new candidate. This final search, however, did not yield any further planetary candidates with both S/N and SDE values higher than 6.

In Fig. 3 we show (a) the entire Kepler photometry of Kepler-160, (b–d) close-up inspections of the three transits of the new candidate prior to detrending and (e–g) after detrending, and (h–j) the phase-folded transits of Kepler-160 b and c and of the new candidate. The resulting system parameterization from emcee is shown in Table 1.

Fig. 2 Lomb-Scargle power spectrum and signal detection efficiency spectrum. Panel a: Lomb-Scargle power spectrum of the raw PDCSAP Kepler light curve that shows stellar (and possibly systematic) activity mostly at periods of between about 10 d and 50 d. The peak at roughly 22 d might correspond to the stellar rotation period. Panel b: the activity has been removed from the light curve, which was detrended with a running biweight filter. Panel c: signal detection efficiency spectrum produced with the TLS algorithm reveals a new candidate transit signal at a period of about 378 d.

Fig. 3 Kepler photometry of Kepler-160. Panel a: entire light curve of Kepler-160 with four years of data. Black dots illustrate the PDCSAP flux, and the red line shows our detrending as obtained with the biweight filter. All values have been normalized by the arithmetic mean of the respective Kepler quarter. Green data points highlight the in-transit flux measurements detected with TLS after detrending. Panels b–d: three new candidate transits highlighted in green and labeled with an x. Transits of the previously known planets Kepler-160 b and c are labeled in violet and orange and labeled b and c, respectively. Black dots refer to the Kepler PDCSAP flux without transits, and the red line shows the detrending function obtained with a running biweight window of 2.2 d length. Panels e–g: detrended light curves and best-fit model photometry of the system centered around the candidate transits. Panels h–j: phase-folded light curves of Kepler-160 b, Kepler-160 c, and new candidate, respectively. The red line shows the respective best-fit transit model and the 17-bin running mean is plotted as the blue line.

3.1.1 Contamination from nearby sources

As a first estimate for the possible contamination of the light curve from the three sources visible in the southeast direction on the Kepler CCD pixel map (to the left of the target in Fig. 1), we explored the Gaia DR2 data and found that these three sources have Gaia magnitudes m_G,con ≳ 20.5 compared to m_G,K160 = 14.6 for Kepler-160. The Gaia bandpass is a good proxy for the Kepler bandpass. The magnitude modulus then gives the flux ratio (F_K160∕F_con) between the sources as m_G,con − m_G,K160 ≈ 2.5log₁₀(F_K160∕F_con). Equivalently, F_K160 = F_con × 10^{(m_G,con−m_G,K160)∕2.5} = 229. In other words, the flux from the nearby stars is lower than one percent of that from Kepler-160, but it cannot be fully excluded that an eclipsing binary might mimic our transit signal with its depth of ~ 300 ppm.

Hence, we used the lightkurve software to generate differential CCD images of the pixel region around Kepler-160 in order to constructper-pixel light curves and measure the S/N of the transit signal for each pixel (Bryson et al. 2013). In Fig. 4 we show the S/N landscape of the new transit candidate in the CCD region around Kepler-160. Visual inspection of the different images in panels b, d, and f shows that the offset between the source of the new transit candidate signal and Kepler-160 is less than a pixel. The different images also rule out the possibility that the transit source is due to one of the known stars in the aperture used to extract the light curve (dotted line in Fig. 1). Moreover, Fig. 4 rules out the possibility that the source of the transit is due to one of the three stars just outside the aperture. Figure 4 does not positively identify the transit signal as being on the target star, but it also does not rule it out. It does confirm, however, that the transits are certainly not caused by any of the nearby sources. We adhere to the “innocent until proven guilty” philosophy of the Kepler classification of planet candidates (Morton et al. 2016; Thompson et al. 2018) and identify this transit signal as belonging to a new planet candidate, KOI-456.04.

We conclude that the three stars just outside the Kepler aperture did not induce the transit candidate signal. The signal rather comes from within a confusion radius of no more than 3.98′′ around Kepler-160, which is the edge length of one Kepler pixel corresponds to 3.98′′ on the sky.

Fig. 4 Mean flux and difference S/N of the new transit candidate on the Kepler pixel map around Kepler-160 (center pixel) for Q5 (a, b), Q9 (c, d), and Q13 (e, f). The positions and Kepler magnitudes of Kepler-160 and two nearby sources are indicated with dots and labels. The difference S/N maps show that the transits are not caused by the nearby sources.

3.1.2 Planet vetting: astrophysical false positives

We set the maxrad parameter of vespa equal to ρ  =  3.98′′ and analyzed the phase-folded transit light curves of Kepler-160 b and c and of the transit candidate signal. Planets b and c had been statistically validated before (Lissauer et al. 2014; Rowe et al. 2014), but no FPPs were published in these studies. Morton et al. (2016), using vespa, derived FPPs of 3.1 × 10⁻¹ for planet b and 1.1 × 10⁻² for planet c,both of which would not qualify them as statistically validated. We therefore decided to consistently derive the FPPs for all three planets again. For planets b and c, we supplied vespa with the data within ±1 d around the transit mid-point of the phase-folded detrended light curve, which we constructed from the median values for T₀ and P as obtained with our MCMC fitting method. For both planets, we searched for secondary eclipses with TLS, and we used our constraints from the nondetections to set the maximum depth for a possible secondary eclipse (the secthresh parameter in vespa) to 10 parts per million (ppm). For the new transiting planet candidate, we used ± 2 d worth of data around the transit mid-point in the phase-folded light curve and a maximum depth for a possible secondary eclipse of 38 ppm, using constraints obtained with TLS.

For planets b and c and for the candidate, vespa predicts FPP₁ values of 1.95 × 10⁻⁵, 6.76 × 10⁻², and 8.33 × 10⁻², respectively, where the subscript 1 refers to the assumption of each planet being the only planet around the star. While our FFP₁ value for planet c is about a factor of six lower than the value obtained by Morton et al. (2016), our value for planet b is about four orders of magnitude lower than their value. Part of this discrepancy might lie in the fact that Morton et al. (2016) used DR24 data, whereas we used data from DR25. Although there have been substantial improvements in the Kepler Data Processing Pipeline between DR24 to DR25 (Twicken et al. 2016), our detailed detrending of any nontransit activity combined with our stringent upper limits on the presence of secondary eclipses (see Sect. 2.1) probably also contributed to these improvements of the FPP₁ estimates.

For each planet, the corresponding probability of planethood, assuming it were the only transit candidate around Kepler-160, is given as P₁ = 1 −FPP₁. Taking into account the presence of two additional transit sequences in the light curve of Kepler-160, we can use the empirical evidence of the nonrandomness of transiting exoplanets among Kepler for each of the three transit sequences. Stars that are known to have one transiting planet are more likely to have more transiting planets than a randomly chosen star. In particular, the probability of a planet candidate that is detected around a Kepler star known to have two additional transiting planet candidates can be estimated as (Lissauer et al. 2012) (1)

In the cases of Kepler-160 b, c, and the new candidate, this “multiplicity boost” yields statistically corrected FPP₃ values of 3.90 × 10⁻⁷, 1.45 × 10⁻³, and 1.81 × 10⁻³, respectively. The presence of the new nontransiting planet Kepler-160 d (see Sect. 3.2), however, is not taken into account in these calculations.Strictly speaking, the FPPs of Kepler-160 b and c and of the new transiting candidate signal are therefore even smaller than our FPP₃ values.

All these FPP₃ values are significantly below the value of 1% that is commonly used as a threshold for statistical validation (Morton et al. 2016). This means that each of the three transit sequences could be formally considered as belonging to a statistically validated planet with very low astrophysical FPPs due to a blending, grazing, or background eclipsing binary or due to a blend of a transiting planet around an unrelated star near the Kepler CCD location of Kepler-160. The entire picture, however, is more complicated because instrumental or statistical false alarms are possible.

3.1.3 Planet vetting: instrumental and statistical false alarms

Although the new transit candidate is very unlikely to be an astrophysical false positive, its 378 d orbital period is suspiciously close to the Kepler heliocentric orbital period of 372 d, a regime that has been known as a source of the rolling-band artifact (Mullally et al. 2018; Thompson et al. 2018). During the primary Kepler mission, the telescope was rolled every 90 d to keep the solar arrays oriented toward the Sun and the radiator pointed toward deep space (Thompson et al. 2016b). This maneuver meant that the Kepler target stars moved on and off a given set of Kepler detectors with some periodicity that can result in a transit-like feature in the light curve with a period close to 372 d after four such rolls. The technical source of the rolling-band effect is in an aliased high-frequency (GHz) noise in the amplifier circuit of the local detector electronics (Kolodziejczak et al. 2010; Sect. 6.7 in Van Cleve & Caldwell 2016). The exact amplitude and frequency of the resulting artifact depend on the CCD detector channel and the board temperature.

Kepler has 21 science CCD modules and four CCD modules that served as fine-guidance sensors. Each of the 21 science modules has 4 output channels, which yields a total of 84 science channels. We infer from the DV Report of Kepler-160 that this target moved repeatedly over modules/outputs 8/2 (Kepler quarters Q1, Q5, Q9, Q13, Q17), 12/2 (Q2, Q6, Q10, Q14), 18/2 (Q3, Q7, Q11, Q15) and 14/2 (Q4, Q8, Q12, Q16). Most alarmingly, module/output 8/2 has been flagged in the Kepler Instrument Handbook as having a high Moiré amplitude that may affect the interpretation of the resulting time series (Table 13 in Van Cleve & Caldwell 2016). The three transits that we found occurred in Q5, Q9, and Q13, when Kepler-160 was located on this known rolling-band channel (module/output 8/2).

We used Model-Shift to evaluate the possibility that this new transit sequence might be caused by rolling-band artifacts alone or in combination with other statistical or instrumental outliers. We fed Model-Shift with our detrended light curve from which we had removed all the transits of planets Kepler-160 b and c. The resulting light curve thus only contained the three transits of our new candidate. Model-Shift then performed a range of metric tests with the light curve. These metrics have been empirically tested with detections from the Kepler exoplanet search pipeline and simulated data (Thompson et al. 2018). Model-Shift computes the multiple event statistic (MES; Jenkins et al. 2002) for the primary transit sequence (Pri), the second most significant sequence (Sec), the third most significant sequence (Ter), and for the most significant positive (or inverted) transit-like feature (Pos).

Figure 5 is an automated output file from Model-Shift. The table at the top of Fig. 5 contains an overview of the Pri (10.7), Sec (5.03), Ter (4.71), and Pos (5.40) metrics. These values need to be compared to the FA₁ value of 5.39, which is the 3 σ threshold atwhich any event is considered to be significant when the measured systematic noise level is taken into account. Although the primary signal is highly significant, the multitude of near-3 σ features suggests that the rolling-band phenomenon is in fact present in the data. This indicates weaker features in the light curve that have a shape similar to the planet candidate transit signal with almost half the strength (see the bottom panels in Fig. 5). While our planet candidate is significantly stronger, it is still possible, though unlikely, that some combinationof rolling bands and other noise is causing our candidate signal. As for the Shape metric in the table at the top of Fig. 5, a value near 0 would suggest a transit-like signal and a value near 0.5 a sinusoidal or heartbeat-like signal (see Coughlin 2017). Our value of 0.336 is just slightly above an arbitrary threshold of 0.3, beyond which the entry is highlighted in red. Because our TLS transit-search algorithm searches for transit-like events with the well-known signal detection efficiency metric (SDE; Kovács et al. 2002; Hippke & Heller 2019), we do not consider this flag as critical.

We used the Pri MES value of 10.7 and the orbital period of 378 d to infer a ~ 85% reliability of the candidate against instrumental artifacts using Fig. 22 of Bryson et al. (2019). We conclude that this newly discovered transit signal passes the astrophysical vetting (Sect. 3.1.2) and that it has a reasonable reliability against instrumental artifacts. Formally speaking, it remains an (unvalidated) planet candidate, however, because validation would require both FPP₃ < 1% and reliability against instrumental artifacts >99%.

Fig. 5 Automated candidate vetting with Model-Shift. The MES values for the Pri, Sec, Ter, and Pos events along with other metrics are listed in the table at the top. The plots show (from top to bottom) the phase-folded light curve without the transits of planets b and c (red: data; blue: binned data; black: model), a zoom into the binned data light curve, and six subpanels illustrating the phase-folded primary, odd, even, secondary, tertiary, and positive events.

3.1.4 Characterization of the transiting habitable zone candidate

We find that the stellar radiative energy flux of the new transit candidate would be times the insolation atthe top of the Earth’s atmosphere. An Earth-like Bond albedo of 0.3 would result in a globally averaged surface temperature of about K (see Fig. 6). This value isjust 10 °C lower than the nominal value of 255 K (− 18 °C) for Earth computed from its net absorption of the top-of-the atmosphere solar flux in thermal equilibrium and in absence of an atmosphere.These calculations neglect the additional heating of the greenhouse effect, for instance, from water vapor (H₂ O), carbon dioxide (CO₂), methane (CH₄), nitrous oxide (N₂O), ozone (O₃), and other greenhouse gases. On Earth, the greenhouse effect amounts to about + 33 °C. If the transit candidate signal belongs to a genuine planet and if this planet has an atmosphere that provides an Earth-like greenhouse effect, then the globally averaged surface temperature is around + 5 °C.

In Fig. 7a we show the distribution of all known transiting exoplanets in the context of the conservativestellar habitable zone (HZ, green shaded area; Kasting et al. 1993) around their individual stars, 3038 planets in total. Most transiting exoplanets are likely gas-dominated giant planets, their median (mean) radius is 2.3 R_⊕ (4.3 R_⊕) or 0.6 (1.1) Neptune radii. The hot limit of the conservative HZ is given by the runaway greenhouse effect derived from one-dimensional radiative-convective cloud-free climate models of water-rich, rocky planets (Kopparapu et al. 2013). The cold limit is imposed by the maximum greenhouse effect that could be provided by a CO₂-rich atmosphere at the limit of condensation. Only a few transiting planets are known to orbit their star in the habitable zone.

Figure 7b shows a zoom into the region of planets with Earth-like insolation around Sun-like stars. In this region of the conservative HZ, the new transiting candidate (denoted with an arrow) would be the smallest planet by far and the only super-Earth with a radius smaller than two Earth radii.

Fig. 6 Posterior distribution of the surface temperature of the new transiting candidate as derived from 150 000 MCMC samplings of the Kepler light curve of Kepler-160. We assumed an Earth-like Bond albedo of 0.3. The median value of 244.8 (+2.1, −2.9) K is indicated with a vertical marker, as is the Earth’s effective global mean surface temperature as calculated from its top-of-the atmosphere solar insolation, that is, neglecting the atmospheric greenhouse effect.

3.2 Nontransiting planet Kepler-160 d

In Fig. 8 we show the TTVs of Kepler-160 b and c as derived with our fittings of a one-planet transit model to each individual transit¹¹. The decrease in the width of the error bars after BKJD = 1000 d is due to the availability of short-cadence data from that point on, whereas earlier measurements were taken in long-cadence mode. The TTVs of planet b do not show any significant periodicity, but the periodicity of the TTVs of planet c is very pronounced. The TTV amplitude of approximately 20 min and the TTV period of about 879 d is in good agreement with earlier studies by Holczer et al. (2016) and Ofir et al. (2018).

We also inspected the Lomb-Scargle periodograms of the TTV series of both planets (not shown). For Kepler-160 b we find no clear signal, and for planet c we find that the 879 d period is the strongest periodic signal.

In Fig. 9 we visualize the results of our three-body TTV modelling including the star Kepler-160, planet c, and one perturber near a low number period commensurability to planet c. Panel a shows the posterior distribution of the mass of the planet (M_d) that would act as a gravitational perturber to cause the observed TTVs of Kepler-160 c. We tested planets near various period commensurabilities (see labels at the top) and find that M_d ranges between about 1 and 100 M_⊕ within the error bars, details depending on the period of the planet. Our results show that various solutions for a nontransiting planet exist and that more data are required before a unique solution can be determined. At this point, our purpose therefore is not a full study of the parameter space, and we accordingly restrict our analysis to the selection of period commensurabilities shown in Fig. 9. Some periods have two mass solutions, where the emcee optimization resulted into two solutions with similar likelihood. Generally speaking, a planet interior to the orbit of Kepler-160 c could have M_d ~ 1 M_⊕, whereas a planet in an outer orbit might be much more massive. These mass estimates of the nontransiting disturber are sensitive to its orbital inclination with respect to Kepler-160 c. Our results are only formally valid for a hypothetical perturber transiting very near our line of sight. That said, the order of magnitude of M_d would still be similar if the mutual inclination of planet c and the perturber were small.

We also tested the hypothesis that the 879 d period in the TTV series of Kepler-160 c is caused by the same planet that is also responsible for the new transit candidate signal detailed above. These attempts, however, did not produce a consistent solution, which prevents us from associating the new transit candidate with the observed TTVs and thus confirming it. In summary, all our simulations suggest that the perturber responsible for the observed TTVs of Kepler-160 c is (a) deep in the planetary mass range and (b) cannot be associated with either the known transiting planet Kepler-160 b or the new transiting candidate. We therefore conclude that there is another nontransiting planet in the system, which we label Kepler-160 d, and which might be described by one of the solutions shown in Fig. 9.

Figure 9b predicts the stellar radial velocity (RV) amplitude that would result from the nontransiting planet Kepler-160 d in a near line-of-sight orbit. These values range between about 1 and 20 m s⁻¹, which could make Kepler-160 d detectable using ground-based stellar spectroscopy.

For the mean P and M_d values shown in Fig. 9 we calculated the transit times of Kepler-160 c over roughly 7000 d, that is, until about March 2028, with our TTV model, the results of which are shown in Fig. 10. Panel a shows the TTVs of all configurations that we considered and panel b shows a zoom into the final 2000 d of the predictions for a selection of configurations. We find that many of our TTV predictions for Kepler-160 c are very similar to the 1:2 MMR (solid black line). Most notably, the TTV divergence of the 2:1 (blue line), 2:3 (green line), and 9:4 (purple line) MMRs exceeds the typical error of the Kepler short-cadence data for a substantial part of our simulations. These configurations might therefore most easily be tested with follow-up transit observations of Kepler-160 c. Alternatively, stellar RVs might be used.

Fig. 7 New transiting candidate KOI-456.04 in the context of the known extrasolar planet population. The position of the open circles indicates the stellar flux received at the top of the atmosphere and the effective temperature of the host star of all known transiting planets as listed in the Exoplanet Encyclopaedia (http://exoplanets.eu) on 29 April 2020. The sizes of the circles scale with the radii of the planets. The green shaded area denotes the habitable zone within the limits of the maximum greenhouse effect (right border) and the runaway greenhouse effect (leftborder). For details, see Sect. 3.1.4. Panel a: the rectangle outlines the margins displayed in panel b, and the arrow points at the new transiting candidate. Panel b: the Earth’s position and size are indicated for reference. The new candidate is labeled with an arrow, and error bars refer to the 1 σ confidence intervals from our MCMC fits. Circles are to scale within each panel but not across panels.

Fig. 8 Panel a: TTVs of Kepler-160 b; panel b: TTVs of Kepler-160 c as derived from our one-planet transit fits to the two transit sequences.

Fig. 9 Results of our MCMC fitting procedure. Panel a: expected mass of the nontransiting planet Kepler-160 d as a function of its orbital period. Dots signify mean values, and error bars illustrate the posterior distribution with 34% of the sample contained in both positive and negative directions from the mean, respectively. The scaling of the ordinate is logarithmic. The orbital periods of Kepler-160 b and c and the orbital MMRs with planet c that we tested are indicated with vertical dashed lines (see labels). Two solutions with very similar probabilities exist for both the 5:2 and 7:2 MMRs. Panel b: expected RV amplitude of the host star due to the orbital motion of Kepler-160 d as a function of P.

4 Discussion

4.1 Follow-up

The current state of technology and the available astronomical facilities make it extremely challenging to validate or reject the new transit candidate KOI-456.04. In principle, such a validation could be achieved through additional transit observations, stellar RV measurements, or by measuring TTVs of the transiting planets that can be explained by a planet with an orbital period of 378 d.

The resulting RV amplitude of its solar-mass host star (Johnson et al. 2017) can be estimated to range between about 0.3 m s⁻¹ for an ocean world and 1.2 m s⁻¹ for a planet with Earth-like composition. Although these RV amplitudes might be accessible with the Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2010) at the Very Large Telescope, Kepler-160 itself is not observable from Chile. This transit candidate could nevertheless be a suitable target to revisit with future spectrographs on the Northern Hemisphere constructed for the RV detection of super-Earths in Earth-like orbits around Sun-like stars. This would be an important stepping stone to confirm even lower-mass Earth-mass planets that will need to be followed-up after their expected discovery with the PLATO mission (Rauer et al. 2014), with an expected launch in 2026.

If Kepler-160 were reobserved with PLATO, then it would not nominally be in any of the stellar samples (P1 to P5), their limiting apparent magnitudes being m_V = 11 and 13, respectively¹². That said, dedicated imagette analysis within the guest observer program might be possible. The potential for a recovery of the candidate signal of KOI-456.04 will depend on the actual location of the star on the PLATO field of view. PLATO will have 24 normal cameras (plus 2 fast cameras), and any star will be observed by either 6, 12, 18 or 24 cameras synchronously. The PLATO performance and noise behavior is expected to be comparable to that of Kepler for the stars that will be observed with all 24 cameras, that is to say, in the center of the field of view. The random noise of an m_V = 11 star is expected to be 34 ppm h⁻¹ for a combination of data from 24 cameras. In the optimal case of a coverage by all 24 cameras, the flux of Kepler-160 would be 10^{(11−14.6)∕2.5} = 3.6% compared to that from that from an m_V = 11 star. This means that the random noise will be 34 ppm h⁻¹/ ppm h⁻¹, which is slightly more than half of the measured transit depth (see Appendix A). Such a configuration or, equivalently,the coverage with 12 cameras and monitoring of 2 transits, might therefore in principle be able to marginally recover the transits of KOI-456.04 if it genuinely is a planet.

In Table 2 we predict the next nine transits of this transit candidate based on our MCMC sampling of the entire light curve. We find very long transit durations (about 13 h) for this candidate, which allows many measurements during a single telescope visit. This fact compensates to some extent for the relative faintness (G = 14.6160, J = 13.460, H = 13.101, K = 13.052) (Cutri et al. 2003; Gaia Collaboration 2018) of this star. The faintness is typical for the Kepler sample. We also find that the uncertainty in the transit times increases by about 2 h between transits, that is, at a rate of roughly 2 h yr⁻¹. Consequently, observations of one of the upcoming putative transits listed in Table 2 will be very important to secure or reject the planet interpretation of this signal.

Fig. 10 Predicted TTVs of Kepler-160 c for various MMRs of the nontransiting planet Kepler-160 d. Times of the transit mid-points are given along the primary x-axis in units of the BKJD and on the secondary x-axis in units of the Coordinated Universal Time (UTC). For reference, BKJD = 7000 d corresponds to 2 March 2028. Panel a: predicted TTV curves for all resonances shown in Fig. 9. The color-coding of the curves is described in the legend. Panel b: zoom into the final ~ 2000 d of our simulations. For clarity, the curves are restricted to some of the most diverging resonances (see labels) and the 5:2, 5:3, 7:2, 7:3, and 7:4 MMRs have been excluded.

4.2 Planet interpretation of the new transit candidate

The radius of of the transit candidate suggests that it might be surrounded by a voluminous layer of volatiles such as H or He (Rogers 2015). Although direct mass measurements of the candidate are currently not available and the mean density is unknown, planetary evolution tracks (Lopez & Fortney 2014) suggest that any potential hydrogen or helium atmosphere would contribute less than about 0.75% to its total planetary mass. Planetary structure models (Fortney et al. 2007) suggest thatfor an Earth-like (or, alternatively, a purely rocky) composition the planetary mass would be about or 13.5 M_⊕ (or 9.9 M_⊕). If the transit candidate, however, were an ocean world composed of equal amounts of rock and water, then its mass would only be about 3.5 M_⊕ for the measured radius.

Although these properties do not seem to make it very Earth-like, its putative super-Earth mass and estimated age of billion years(Johnson et al. 2017) could nevertheless make it a superhabitable world (Heller & Armstrong 2014), potentially with even better prospects for life than truly Earth-like planets.

4.3 Alternative stellar parameterization

The T_eff value of K from the Gaia DR2 that we used for the purpose of consistency throughout our analysis is somewhat lower than the values of K determined from HIRES high-resolution spectroscopy (Johnson et al. 2017) and K derived from RCSpec mid-resolution spectroscopy (Everett et al. 2013). The stellar radius of from Gaia DR2 is well within the uncertainties of as measured based on HIRES spectroscopy (Johnson et al. 2017), although it is somewhat larger than the value of 1.04 (± 0.05) R_⊙ estimated via RCSpec spectroscopy (Everett et al. 2013). The stellar effective temperature of 5857 (± 200) K and radius of 0.884 (± 0.430) R_⊙ given in the Kepler Input Catalog (Brown et al. 2011; Mathur et al. 2017) (available at the NASA Exoplanet Archive)¹³ were derived using five-band photometry without additional stellar spectroscopy or model-independent parallax measurements. Instead, the stellar extinction between Earth and the source was estimated and used as a free parameter, which is why we consider these values as less reliable.

The Gaia DR2 value of the stellar radius that we used is the highest of all the sources discussed above. As a consequence, and because R_p scales linearly with R_s, we derive relatively large planetary radii. If we were to use best-fitting stellar radii of 1.11 R_⊙ (Johnson et al. 2017), 1.04 R_⊙ (Everett et al. 2013), or 0.884 R_⊙ (Brown et al. 2011; Mathur et al. 2017), for comparison, the resulting planetary radius for candidate transit signal would be 1.89 R_⊕, 1.77 R_⊕, or 1.50 R_⊕, respectively.This shows that although the different literature values for the stellar radius all agree within about two standard deviations, the effect on the planetary radius might have dramatic implications for the composition of the planet and its nature as a super-Earth or a mini-Neptune (Rogers 2015).

4.4 Coplanarity of the nontransiting planet Kepler-160 d

In our MCMC optimization procedure of the TTVs of Kepler-160 c (Sect. 2.4) we assumed an orbital inclination of 90° of the test planet, that is to say, a line-of-sight arrangement of its orbit. Strictly speaking, this assumption contradicts our nondetection of transits from our newly discovered planet, Kepler-160 d, which our analysis suggests to have a mass(and we expect thus a radius) in the super-Earth or mini-Neptune regime of maybe between 2 R_⊕ and 4 R_⊕. If the orbital inclination of Kepler-160 d were truly 90°, then we would almost certainly have detected its transits.

A short trigonometric exercise shows that the critical orbital inclination for the disk center of Kepler-160 d to perform a grazing transit is i_g = arccos(R_s∕a). The possible orbital periods of Kepler-160 d between about 5 and 50 d, as suggested by our TTV fitting procedure (see Fig. 9), correspond to orbital semimajor axes between about 0.06 and 0.27 AU. The resulting values for i_g then range between 84°.8 in the innermost orbits and 88°.9 in the outer orbits.

Our results for the mass and semimajor axis of the test planet will only differ by a few percent if its actual inclination differs by only a few degrees from our line-of-sight assumption. We argue that it is very plausible for Kepler-160 d to be on a near-coplanar orbit with planet c given that planet c is essentially coplanar with planet b, given the near-coplanarity of the Solar System planets, and given the many detections of near-coplanar multiplanet transiting systems from Kepler (Rowe et al. 2014).

4.5 Stability of the nontransiting planet Kepler-160d

We confirm the stability of the planetary system Kepler-160 b/c upon inclusion of the nontransiting planet Kepler-160 d near the period commensurabilites to planet c modeled in Sect. 3.2. We integrated the system with the N-body code Mercury (Chambers 1999) over 1 Myr with a time step of 0.2 days (roughly a twentieth of the period of Kepler-160 b). The system is more than a thousand times older than our integration time, but we had to limit ourselves to 1 Myr to achieve results for the numerous configuration within a reasonable computation time. A more detailed stability analysis with longer integration times and over a range of orbital inclinations is reasonable when the number of possible orbital solutions has been narrowed down to a few or when a unique solution for the mass and period of Kepler-160 d has been found.

5 Conclusions

We have analyzed the Kepler light curve of the solar-type star Kepler-160, which has been known to host two transiting planets. We detect a new sequence of three periodic transit-like features that we attributed to a 1.9 R_⊕ planet with virtually Earth-like insolation from a Sun-like star. The S/N of this signal is 9.5, the signal detection efficiency (SDE_TLS) value is 16.3, the astrophysical FPP when the planet multiplicity boost is accounted for (FPP₃) is 1.81 × 10⁻³, the multiple event statistic (MES) is 10.7, and the reliability against Kepler’s rolling-band artifact and other instrumental artifacts is ~ 85%. We also showed that the signal is not caused by any nearby stellar contaminant. All of these metrics point toward a genuine planet as the causeof the observed transit sequence. For now, we remain cautious because there is mild evidence of the rolling-band artifact in the light curve and refer to it as a planet candidate KOI-456.04. Our MCMC fits to the light curve suggest that this planet, if real, could have a global average surface temperature of around + 5 °C assuming anEarth-like Bond albedo and an Earth-like atmospheric greenhouse effect. We present a list of the next nine transit times until December 2028 (Table 2). The next transit is projected to occur on 14.09.2020 (UTC).

We analyzed the previously known TTVs of Kepler-160 c in trying to identify our new transit candidate as the source of the observed TTVs and thereby confirm it. Our MCMC analysis of this system with a three-body (one star plus two planets) TTV model did not result in a solution that involves the transiting candidate as a source of the observed TTVs of planet c, however. Instead, we found that Kepler-160 is orbited by a nontransiting planet, which we refer to as Kepler-160 d. The sinusoidal nature of the observed TTVs suggests that Kepler-160 d is close to a low-integer MMR with Kepler-160 c. Under the assumption of a low mutual orbital inclination, we find that the resulting mass of this new planet Kepler-160 d is always within the planetary mass range, with specific values between about 1 and 100 Earth masses.

Acknowledgements

The authors thank Jeffrey Coughlin for helpful discussions of the planetary nature of the new transit signal. We also thank an anonymous referee for their valuable review of our manuscript. 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 work made use of NASA’s ADS Bibliographic Services. This research has made use of “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute forAstronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, 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. R.H. is supported by the German space agency (Deutsches Zentrum für Luft- und Raumfahrt) under PLATO Data Center grant 50OO1501.

Appendix A Detrending

Here we address the question of whether the new candidate signal that we detected might be due to the use of a particular detrending function. Prior to our nominal transit search and subsequent MCMC fitting of the light curve, we removed any nontransit variability on timescales longer than the expected transit duration (see Sect. 2.1) using the Wōtan implementation (Hippke et al. 2019) of Tukey’s biweight filter (Mosteller & Tukey 1977). We refer to the resulting detrended flux as F_biw, as shown in Figs. 3e–g.

To explore the effect of the detrending algorithm on the transit signal, we also detrended the quarterly normalized PDCSAP flux (see Fig. 3a) using the other detrending algorithms from Wōtan. As examples, we focus on the cosine filter, resulting in a flux series F_cos, the locally weighted or estimated scatter plot smoothing (LOWESS; Cleveland 1979, 1981), providing a flux series F_low, and the running median filter as used in Heller et al. (2019a), producing a flux series F_med. In Fig. A.1 we show the differential flux series between the detrended light curve obtained with the biweight filter and F_cos (top row), F_low (center row), and F_med (bottom row).

The overall variability of the differential flux is up to 200 ppm, but mostly < 100 ppm, on a timescale of about 10 d. On timescales comparable to the duration of our new transit candidate signal of about 0.5 d, however, the variability is typically on the order of a few times 10 ppm. In particular, we find that the differential effect between thedetrending methods is about 20 ppm during the first and second transit both, and it is up to 45 ppm during the third transit. This excursion is labeled in the top right panel in Fig. A.1. All things combined, the overall effect is small, although potentially not negligible, compared to the transit depth of about 300 ppm.

Moreover, we find that the biweight filter that we used to prepare the light curve for our subsequent MCMC characterization of the system always resulted in transit depths that are a few times 10 ppm shallower than obtained with the other detrending algorithms. As an example, we note that the differential in-transit flux during the third transit shown in Fig. A.1 is slightly increased (F_biw − F_x > 0), with F_x being any of the three flux series from the alternative detrending methods, which means that F_biw > F_x. In words, the other detrending algorithms produced slightly (~10 ppm) deeper signals during the third transit than the biweight filter that we used for our nominal system characterization. In other words, the transit depth (and likely also the signal detection efficiency) that we derived can be regarded as conservative estimates.

We conclude that the signal appears repeatedly even when the detrending filter is changed. The effect of the filter is about a few times 10 ppm, which is small compared to the transit depth of about 300 ppm.

Fig. A.1 Differential flux of different detrending methods compared to Tukey’s biweight filter. The three panel columns refer to the three transits, and the in-transit data points, as detected by TLS, are highlighted in green as in Figs. 3a–d. The detected transit width is also indicated with a horizontal bar in the center of each panel. The maximum in-transit differential flux that we found is 45 ppm, which occurs during the third transit and for Fbiw −Fcos in the top right panel (see label). For comparison, the transit depth of our new candidate signal is about 300 ppm.

References

All Tables

All Figures

Fig. 1 Panel a: Kepler CCD photometry of Kepler-160 from the Data Validation Report. The dotted line outlines the aperture mask used to extract the light curve. The cross denotes the target position, and the stars (with annotated KIC numbers) indicate the positions of stars that might act as contaminants of the target light curve. The color scale illustrates the difference between mean flux out-of-transit and in-transit of Kepler-160 b (KOI-456.01) after normalizing by the uncertainty in the difference for each pixel. North and east are indicated with yellow lines. The Kepler CCD pixels are squares and appear stretched in this figure from the Data Validation Report. Panel b: Pan-STARRS image that shows Kepler-160 in the center. The imageis properly oriented so that north coincides with the Kepler image in panel a. The distance of 6.3′′ to one nearby star is indicated with an arrow.

In the text

Fig. 2 Lomb-Scargle power spectrum and signal detection efficiency spectrum. Panel a: Lomb-Scargle power spectrum of the raw PDCSAP Kepler light curve that shows stellar (and possibly systematic) activity mostly at periods of between about 10 d and 50 d. The peak at roughly 22 d might correspond to the stellar rotation period. Panel b: the activity has been removed from the light curve, which was detrended with a running biweight filter. Panel c: signal detection efficiency spectrum produced with the TLS algorithm reveals a new candidate transit signal at a period of about 378 d.

In the text

Fig. 3 Kepler photometry of Kepler-160. Panel a: entire light curve of Kepler-160 with four years of data. Black dots illustrate the PDCSAP flux, and the red line shows our detrending as obtained with the biweight filter. All values have been normalized by the arithmetic mean of the respective Kepler quarter. Green data points highlight the in-transit flux measurements detected with TLS after detrending. Panels b–d: three new candidate transits highlighted in green and labeled with an x. Transits of the previously known planets Kepler-160 b and c are labeled in violet and orange and labeled b and c, respectively. Black dots refer to the Kepler PDCSAP flux without transits, and the red line shows the detrending function obtained with a running biweight window of 2.2 d length. Panels e–g: detrended light curves and best-fit model photometry of the system centered around the candidate transits. Panels h–j: phase-folded light curves of Kepler-160 b, Kepler-160 c, and new candidate, respectively. The red line shows the respective best-fit transit model and the 17-bin running mean is plotted as the blue line.

In the text

	Fig. 4 Mean flux and difference S/N of the new transit candidate on the Kepler pixel map around Kepler-160 (center pixel) for Q5 (a, b), Q9 (c, d), and Q13 (e, f). The positions and Kepler magnitudes of Kepler-160 and two nearby sources are indicated with dots and labels. The difference S/N maps show that the transits are not caused by the nearby sources.
In the text

Fig. 5 Automated candidate vetting with Model-Shift. The MES values for the Pri, Sec, Ter, and Pos events along with other metrics are listed in the table at the top. The plots show (from top to bottom) the phase-folded light curve without the transits of planets b and c (red: data; blue: binned data; black: model), a zoom into the binned data light curve, and six subpanels illustrating the phase-folded primary, odd, even, secondary, tertiary, and positive events.

In the text

Fig. 6 Posterior distribution of the surface temperature of the new transiting candidate as derived from 150 000 MCMC samplings of the Kepler light curve of Kepler-160. We assumed an Earth-like Bond albedo of 0.3. The median value of 244.8 (+2.1, −2.9) K is indicated with a vertical marker, as is the Earth’s effective global mean surface temperature as calculated from its top-of-the atmosphere solar insolation, that is, neglecting the atmospheric greenhouse effect.

In the text

Fig. 7 New transiting candidate KOI-456.04 in the context of the known extrasolar planet population. The position of the open circles indicates the stellar flux received at the top of the atmosphere and the effective temperature of the host star of all known transiting planets as listed in the Exoplanet Encyclopaedia (http://exoplanets.eu) on 29 April 2020. The sizes of the circles scale with the radii of the planets. The green shaded area denotes the habitable zone within the limits of the maximum greenhouse effect (right border) and the runaway greenhouse effect (leftborder). For details, see Sect. 3.1.4. Panel a: the rectangle outlines the margins displayed in panel b, and the arrow points at the new transiting candidate. Panel b: the Earth’s position and size are indicated for reference. The new candidate is labeled with an arrow, and error bars refer to the 1 σ confidence intervals from our MCMC fits. Circles are to scale within each panel but not across panels.

In the text

	Fig. 8 Panel a: TTVs of Kepler-160 b; panel b: TTVs of Kepler-160 c as derived from our one-planet transit fits to the two transit sequences.
In the text

Fig. 9 Results of our MCMC fitting procedure. Panel a: expected mass of the nontransiting planet Kepler-160 d as a function of its orbital period. Dots signify mean values, and error bars illustrate the posterior distribution with 34% of the sample contained in both positive and negative directions from the mean, respectively. The scaling of the ordinate is logarithmic. The orbital periods of Kepler-160 b and c and the orbital MMRs with planet c that we tested are indicated with vertical dashed lines (see labels). Two solutions with very similar probabilities exist for both the 5:2 and 7:2 MMRs. Panel b: expected RV amplitude of the host star due to the orbital motion of Kepler-160 d as a function of P.

In the text

Fig. 10 Predicted TTVs of Kepler-160 c for various MMRs of the nontransiting planet Kepler-160 d. Times of the transit mid-points are given along the primary x-axis in units of the BKJD and on the secondary x-axis in units of the Coordinated Universal Time (UTC). For reference, BKJD = 7000 d corresponds to 2 March 2028. Panel a: predicted TTV curves for all resonances shown in Fig. 9. The color-coding of the curves is described in the legend. Panel b: zoom into the final ~ 2000 d of our simulations. For clarity, the curves are restricted to some of the most diverging resonances (see labels) and the 5:2, 5:3, 7:2, 7:3, and 7:4 MMRs have been excluded.

In the text

Fig. A.1 Differential flux of different detrending methods compared to Tukey’s biweight filter. The three panel columns refer to the three transits, and the in-transit data points, as detected by TLS, are highlighted in green as in Figs. 3a–d. The detected transit width is also indicated with a horizontal bar in the center of each panel. The maximum in-transit differential flux that we found is 45 ppm, which occurs during the third transit and for Fbiw −Fcos in the top right panel (see label). For comparison, the transit depth of our new candidate signal is about 300 ppm.

In the text