© ESO 2020

1 Introduction

Hitherto, studies of the variability of exoplanetary transit timings (i.e., start/end time and duration) have supposed the form of transit light curves (TLCs) to be symmetric, suggesting a spherical exoplanet (e.g., Holczer et al. 2016). As a result, the independent positions of transit borders are unknown. Furthermore, the borders of the TLC are most sensitive to the shape of the transiting object, which is affected by possible presence of exo-rings. (Schlichting & Chang 2011) and dust formations (Arkhypov et al. 2019; Wang & Dai 2019).

In particular, the variable stellar ionizing radiation could vary the effective size of the shadow of an exoplanet with outflowing (Wang & Dai 2019), apparently dusty (Arkhypov et al. 2019) atmosphere. A circumplanetary disk of dusty plasma, formed in some cases (Khodachenko et al. 2015), could also affect the planetary shadow, resulting in variation of its effective size due to reconnection events in the disk. Another mechanism causing variability in transits may be related with shocks and material flows around hot Jupiters (Llama et al. 2013). All these factors, as well as precessing exo-rings, might result in variability in transit borders. Here, we therefore perform independent measurements of transit border timings and TLC minimum time and investigate variations thereof.

In that respect, TLCs provided by the Kepler space telescope appear to be a promising but still superficially studied source of information concerning possible variability in exoplanetary transits. So far, such variability has only been studied for decaying rocky planets (Sanchis-Ojeda et al. 2015; Garai 2018 and therein) and asteroids (Rappaport et al. 2016 and therein). In this paper, we analyze for the first time the TLC variability of different types of transiting Kepler objects of interest (KOIs). The maximal duration of available data records and their high precision make Kepler data the best choice to search for tenuous signatures of variability in TLCs.

For the analysis we apply a specially elaborated method, which is explained in Sect. 2. The obtained results are described in Sect. 3, discussed in Sect. 4, and concluded in Sect. 5.

2 Method and stellar set

2.1 Transit light curve analysis

We use publicly available light curves from the Kepler mission (NASA Exoplanet Archive, thereafter NASA EA¹), after pre-search data conditioning(light flux F_PDC hereinafter),where instrumental drifts, focus changes, and thermal transients are removed or suppressed (Jenkins et al. 2010). Our survey relies on long-cadence data with a photon accumulation (i.e., exposure) period of δt_L = 29.4 min, or 0.02 days. This approach provides the highest precision of the input TLCs, although it cannot resolve the details of ingress and egress parts. The short-cadence data with an exposure period of δt_S = 1 min, in spite of a higher time resolution, are available for a lesser number of stars and not for all quarters;and they also increase the photometry error due to photon noise by a factor of $\sqrt{\delta t_{\text{L}}/\delta t_{\text{S}}}=5.4$.

To prepare the analyzed light curves for determining TLC timing parameters we first remove the residual instrumental drifts as well as the stellar variability at timescales longer than the transit duration Δt_tr, taken from the NASA EA. This procedure consists of the following steps: (a) Extraction of a fragment of light curve, enclosed within a time interval of ± 10Δt_tr centered at the particular transit (see Fig. 1a). (b) Normalization of the extracted fragment: F_PDC∕⟨F_PDC⟩, where ⟨F_PDC⟩ means an average flux over the light-curve fragment. (c) Removal of the transit from the analyzed light-curve fragment, including the half-exposition margins of ± 0.01 days, based on the published transit duration values (NASA EA). (d) Approximation of the remaining light-curve fragment without transit with a sixth-order polynomial: F_b(t) = at⁶ + bt⁵ + ct⁴ + dt³ + et² + ft + g, where t is the time of flux measurement, and a, b, c, d, e, f, g are the fitted coefficients. This approximation is an iterative process (ten iterations), with consequent exclusion of remaining stellar flares, light-curve artifacts, and residual transit effects above the threshold |[F_PDC(t)∕⟨F_PDC(t)⟩] − F_b(t)| > 3σ_b, where σ_b is a standard deviation from the approximation F_b(t). Since the final approximation F_b(t) is insensitive to the transit (Fig. 1a), we consider it as a reference level with which to determine the flux decrease during the transit: ΔF_k = [F_PDC(t)∕⟨F_PDC(t)⟩] − F_b(t), which is used in the analyses below.

In the long-cadence data records with the above-mentioned time step δt_L, there are only a few flux measurement points n = Δt_tr∕δt_L ~ 5 on the time interval of transit duration with a typical value of Δt_tr ~ 0.1 day. This is insufficient for our study of the transit border timing. We therefore increased the number of data points up to K = NΔt_tr∕δt_L ~ 100 by combining the measurements of N = Δt_w∕P_tr transits within a time-window of width Δt_w centered at t_w into one phase-folded TLC ΔF(Δt). Here, Δt ≡ t − t_E is the time counted for each transit with number E = 0, 1, 2, ... relative to its mid-time t_E = t₀ + P_trE. The reference transit time t₀ and the transit period P_tr are taken in their cumulative versions from the NASA EA. The used time-windows are defined as a sequence of adjacent frames, so that t_w = t₀ + Δt_wJ_w, where J_w = 0, 1, 2 ... is the number of a particular time-window.

The creation of a folded TLC by combining individual transits within a window is also an iterative process aimed at exclusion of light-curve artifacts, stellar flares, and starspot eclipses. In order to clean the analyzed TLC of outliers related with the above-mentioned effects, we approximate it with a fourth-order polynomial $\Delta F\approx \Delta F_{\text{ap}}(\Delta t)=a_{\text{ap}}{(\Delta t)}^4+b_{\text{ap}}{(\Delta t)}^3+c_{\text{ap}}{(\Delta t)}^2+d_{\text{ap}}\Delta t+e_{\text{ap}}$, where a_ap, b_ap, c_ap, d_ap and e_ap are the fitted coefficients. All the flux values with |ΔF(Δt) − ΔF_ap(Δt)| > 3σ_ap are excluded from further consideration and the whole approximation procedure is repeated with a new estimate for σ_ap in order to remove another group of outliers in the TLC. After ten rounds of such iterations, the final folded TLC is sufficiently smooth for further analysis. An example of such a folded and iteratively smoothed TLC is shown in Fig. 1b.

For the correlation analysis it is important to provide an acceptable number (≥ 3) of phase-folded TLCs (i.e., the number of time-windows N_w over the whole available Kepler data record of about 4 yr), whilst keeping a sufficient number of flux measurement points in each folded TLC,K, for the transit borders study. We note that, in the case of short transit periods (P_tr ≲ 10 days), the value of K can be easily kept up to 200, whereas for intermediate (10 ≲ P_tr ≲ 20 days) and long (P_tr  ≳ 20 days) periods the number K was decreased to 100 and 60, respectively. These values were adopted in the regular survey of the whole target list in Sects. 3.1 and 3.2, whereas during the analysis of particular objects in Sect. 3.3 the window parameter N has been additionally varied to optimize the time resolution for optimal visualization of the transit variability.

To determine the parameters of a TLC with high accuracy, we separately approximate the start-, end-, and minimum- parts of the obtained folded TLCs with second-order polynomials. In particular, the minimum-part of the TLC is approximated as ΔF ≈ F_a = gΔt² + hΔt + q within the interval |Δt| < Δt_tr∕4 centered at the folded TLC mid-time (Δt = 0). An additional reason for using the second-order polynomial here is to avoid the uncertainty associated with selection among multiple extremes, which are typical for higher-order polynomial approximations. Moreover, the single minimum at Δt_m = −h∕(2g) is considered also as the TLC minimum-time, which depends mainly on the stellar limb darkening. This method reduces the sensitivity of the TLC timing parameter Δt_m to possible local minima caused by starsports. Correspondingly, the depth of the TLC minimum can be found as $\Delta F_{\text{tr}}=g\Delta t_{\text{m}}^2+h\Delta t_{\text{m}}+q$.

For the determination of TLC border timing parameters, that is, its start and end time moments, indicated with the indexes “s” and “e”, respectively, we use second-order polynomial approximations of the reverse dependence Δt(ΔF) = a_s,eΔF² + b_s,eΔF + c_s,e within the interval 0.5ΔF_tr < ΔF < 0 separately for the start- and end- parts of the TLC (i.e., indexes “s” or “e”, respectively). Correspondingly, two different ranges of Δt, namely (i) − Δt_tr∕2 − δt_L∕2 < Δt < 0 and (ii) 0 < Δt < Δt_tr∕2 + δt_L∕2, are considered (see in Fig. 1b). We note that the above specified ranges for the ΔF and Δt are jointly used for the extraction of flux count points involved in the approximations of TLC start- and end-parts. As the transit border must be located at ΔF = 0, the approximation polynomial coefficient c_s,e directly gives the border times Δt_s = c_s and Δt_e = c_e of the folded TLC defined within a particular time-window centered at a certain t_w. To study variability of the parameters ΔF_tr, Δt_m, Δt_s, and Δt_e, we similarly associate them with the corresponding central time t_w of the time-window for which the folded TLC was obtained. The standard errors of ΔF_tr, Δt_m, Δt_s and Δt_e (Fig. 2) are expressed in terms of the definition errors of the coefficients of the approximating polynomials.

As the quadratic approximations of TLC border parts are the simplest ones, they also provide an acceptable fit for the transit border regions of interest above the inflection point (Fig. 1b). The cubic and higher order polynomials increase the dispersion of Δt_s and Δt_e times defined for different time-windows, i.e., windows with different J_w. Although the polynomial method could result in systematic errors in the defining of Δt_s and Δt_e, such regular displacements do not affect the detection of transit variability. Further below, the transit shape parameter A_s =(Δt_m − Δt_s)∕(Δt_e − Δt_s) and the Pearson correlation coefficients r_se between Δt_s and Δt_e, as well as analogous cross-correlations r_me (between Δt_m and Δt_e), r_ms (between Δt_m and Δt_s), r_As (between A_s and Δt_s), r_Ae (between A_s and Δt_e) and r_Am (between A_s and Δt_m), are used for the phenomenological analysis of the folded TLCs.

Fig. 1 Example of light-curve processing (the star KIC 11414511). Panel a: approximation of the transit background (red squares) with a sixth-order polynomial Fb(t) (green curve). Panel b: folded profile ΔF(Δt) of 25 adjacent individual transits. Blue curves are the applied polynomial approximations used for measurements of labeled parameters.

2.2 Compiled data set

The considered data set of 98 KOIs and the results obtained for these are presented in Table A.1. Our targets were selected among the most qualitatively observed objects from the list by Aizawa et al. (2018) and supplemented by other Kepler transiting bodies, for which only long-cadence light curves are available, but with high signal-to-noise ratios (S∕N ≳ 1000) according to NASA EA. To avoid the inter-transit interference in multi-planetary systems, only the cases with single transiting object (i.e., transiter) were included in our target list.

The KOIs with current false positive status were not considered, apart from KOI 125.01, KOI 631.01, and KOI 1416.01, which were reclassified as false positives already after the compilation of our data set. The status of these objects is marked in Table A.1 with “n”, that is, probable non-planets. The same mark is assigned to KOI 823.01, KOI 971.01, KOI 1154.01, KOI 6085.01 and KOI 6774.01 which have unknown masses and anomalously high radii, R_p > 24 R_e (here R_e stays for the radius of Earth), exceeding the empirical range of the radius–mass relation (see Fig. 1 in Ulmer-Moll et al. 2019). Altogether, the highest status “p!” of confirmed exoplanets with measured masses below 13 M_J (here M_J stays for the mass of Jupiter) is assigned to 40% of the objects in the considered set of 98 KOIs. One other object, KOI 680.01, was also given the status “p!”. This is a so-called super-puff which has an obviously planetary mass of 0.84 ± 0.15 M_J (Almenara et al. 2015), in spite of a relatively large radius $R_{\text{p}}={26.5}_{-9.4}^{+3.6}$ R_e. Confirmed planets with unknown masses are marked with “p” (35% of cases), whereas unconfirmed candidates were given the status index “c” (17%). To check the status of objects, we additionally used the Extrasolar Planets Encyclopedia² and the SIMBAD Astronomical Database³. The objects for which the previously suggested false-positive or stellar eclipse status has been canceled in NASA EA are marked with the sign “?” in our Table A.1. Finally, the status mark “*” signifies that the KOI formally has the current flag of stellar eclipse false positive, but the measured mass of the object contradicts its assumed stellar nature.

There isa special notation in the KOI cumulative list of NASA EA for KOIs “that are observed to have a significant secondary event (i.e., the secondary TLC minimum due to the eclipse of the planet by the star), transit shape, or out-of-eclipse variability” indicating that the transit-like event is most likely caused by an eclipsing binary. “However, self-luminous, hot Jupiters with a visible secondary eclipse will also have this flag set, but with a disposition of PC (i.e., planetary candidate)”.

For example, in our set, the objects KOI 13.01 and KOI 18.01 (Kepler-13b and -5b) orbiting the stars KIC 9941662 and KIC 8191672, respectively, have the stellar eclipse false positive flag in the KOI cumulative list of NASA EA. Nevertheless, the measurements of their masses by two different methods give non-stellar values (Esteves et al. 2013) corresponding to the classification of objects as confirmed planets in the same NASA EA. As to the secondary eclipses, they are an important instrument of exoplanet meteorology used in many dedicated studies of genuine hot Jupiters (e.g., Jackson et al. 2019; Pass et al. 2019).

The presence in our data set of the admixture (8%) of objects with the supposed non-planetary status “n” is motivatedby the fact that non-planetary KOIs need to be studied as well as exoplanets. Furthermore, the false-positive status is a probabilistic issue (e.g., Morton et al. 2016) and for some objects it changes from one catalog to another. Hence, the current false-positive status of an object does not completely exclude a planetary nature.

Figure 3 demonstrates the ranges of physical parameters (radius R_p of the transiting body and its orbitradius a_p) of the KOIs in the analyzed data set. One can see that our target list includes mainly hot jupiters with an admixture of hot neptunes and super-earths as well as KOI 823.01 and KOI 971.01 of stellar size, classified as planetary “candidates”in the Kepler data of NASA EA. In the case of KOI 1478.01 at KIC 12403119 with maximal P_tr = 76.13 days, which appeared in a kind of resonance with the flux counting period, the parameters Δt_s, Δt_m and Δt_e were found only in two time-windows. This always gives an exact formal correlation, that is, r_se = ±1. Such fictive correlations were excluded from Table A.1.

Fig. 2 Temporal behavior of the transit timing parameters of real TLCs like in Fig. 1b, compiled for a sequence of adjacent time-windows, vs. time in Julian days JD: panel a: transit start-time Δts. Panel b: transit end-time Δte. Panel c: transit duration Δte − Δts. Panel d: transit minimum-time Δtm (maximal flux decrease). Panel e: transit shape parameter As. Panel f: transit flux decrease ΔFtr.

3 Results and their analysis

Table A.1 shows the main results from processing the compiled data set, including the standard deviations σ_m, σ_s and σ_e of transit timings Δt_m, Δt_s, and Δt_e, respectively, as well as the timing cross-correlations r_ms, r_me, and r_se and the average transit shape parameter 〈A_s〉 (see Sect. 2.1 for details). Despite the dispersed errors, Table A.1 contains hidden patterns which enable certain interpretations and conclusions. Below we present some examples of our findings after processing the results, including general statistics, and we highlight several specific individual cases.

Fig. 3 Distribution of selected KOI targets (listed in Table A.1) vs. the radii Rp and orbital semimajor axes ap of transiters according to NASA EA3. Color shows the status of objects: violet – probably non-planets (“n” in Table A.1); blue – candidates (“c” in Table A.1); green – confirmed planets but with unknown masses (“p” in Table A.1); red – confirmed planets with measured planetary masses (“p!” in Table A.1). The cases with “*” and “?” are included. The same color coding is used throughout the paper.

3.1 Transit light curve border timing variability

Figure 2 shows some typical examples of the temporal behavior of the obtained transit timing parameters Δt_s, Δt_m, Δt_e, and ΔF_tr, as well as of their combinations (Δt_e − Δt_s and A_s). We use the following three formal criteria for the detection of the best cases of transit border timing variability:

• (a)

  For Δt_s, the ratio σ_s ∕ε_s > 2.5 holds true, where σ_s is the standard deviation of the estimate of Δt_s over the time-windows, and ε_s is an error of Δt_s measurements averaged over all time-windows.

• (b)

  For Δt_e, the same criterion as for Δt_s is considered, that is, σ_e∕ε_e > 2.5.

• (c)

  For r_se, significant cross-correlation |r_se|∕σ_rse > 2.5 is required, where ${\sigma }_{\text{rse}}=(1-r_{\text{se}}^2)/\sqrt{N_{\text{w}}-1}$ is a standard error of r_se according to Hotelling (1953), and N_w is the number of time-windows in which the phase-folded TLCs are generated.

If any of these criteria hold true, the KOI is considered as potentially interesting. An additional condition of N_w≥ 10 was applied to define the best tracked cases. The results of such a selection are presented in Table A.2.

Figure 4 shows an example of varying transit border timing for KOI 18.01 in the light curve of KIC 8191672 from Table A.2. The plots of Δt_s and Δt_e versus time are combined to demonstrate their clear anti-correlation r_se = −0.73 ± 0.11. All 17 KOIs in Table A.2 show analogous (r_se ≤−0.51) negative cross-correlation. Supposing a stable host star (i.e., with a constant radius), such an anti-correlation effect can have two explanations: (1) a varying impact parameter β of the transiting body (i.e., the minimal distance between the stellar disk center and the transiter’s trajectory projection on the disk, expressed in stellar radii), or (2) a varying effective size of the transiter. To distinguish between these options, let us consider their contribution to the TLC border timing.

Figure 5 shows the transit scheme used to specify the transit border times Δt_s or Δt_e. From the triangle △TMS, one can express the semi-trajectory TM as $x=\sqrt{{(R_{*}+R_{\text{p}})}^2-{(\beta R_{*})}^2}$, where R_* and R_p, i.e., the distances cS and cT, are the stellar and transiting planet radii, respectively. This semi-trajectory is related with the transit border timing as follows: |Δt_s,e| = x∕V_p, where V_p = 2πa_p∕P_tr is the velocity of a visible transiter (across the stellar disk), and a_p is the orbit radius from NASA EA. We note that we assume circular orbits, because all KOIs in our set have zero orbital eccentricities according to the “Kepler candidate overview pages” in NASA EA. Therefore, a small variation Δx of x can be expressed in terms of fluctuations (i.e., deviations from averaged values) Δβ and ΔR_p of the corresponding parameters: $$\Delta x=\frac{\partial x}{\partial \beta }\Delta \beta +\frac{\partial x}{\partial R_{\text{p}}}\Delta R_{\text{p}},$$(1)

where $$\frac{\partial x}{\partial \beta }=\frac{-\beta R_{*}^2}{x},$$(2) $$\frac{\partial x}{\partial R_{\text{p}}}=\frac{R_{*}+R_{\text{p}}}{x}.$$(3)

Let us consider the effect of only the first term in Eq. (1) expressed via standard deviations: ${\sigma }_{\text{s},\text{e}}\equiv \sqrt{\langle {(\Delta x)}^2\rangle }/V_{\text{p}}$. Then Eqs. (1) and (2) give $${\sigma }_{\text{s},\text{e}}=\frac{1}{V_{\text{p}}}\left|\frac{\partial x}{\partial \beta }\right|{\sigma }_{\beta }=\frac{\beta R_{*}^2}{x{V}_{\text{p}}}{\sigma }_{\beta },$$(4)

where ${\sigma }_{\beta }\equiv \sqrt{\langle {(\Delta \beta )}^2\rangle }$ is the standard deviation of the impact parameter β. Correspondingly, for β = 0, σ_s,e= 0, whereas for β → 1 (since R_p ≪ R_*) in view of x→ 0, σ_s,e→∞. Therefore, varying impact parameter β should result in a maximum of σ_s,e at β≈ 1, which means that the KOI objects with most varying TLC border timing should have β close to unity.

Similarly, let us consider now the effect of the second term in Eq. (1), $$\Delta x=\frac{\partial x}{\partial R_{\text{p}}}\Delta R_{\text{p}}=\frac{R_{*}+R_{\text{p}}}{x}\Delta R_{\text{p}}.$$(5)

For R_p ≪ R_*, we obtain from Eq. (5) the standard deviations: $${\sigma }_{\text{s},\text{e}}\equiv \frac{\sqrt{\langle {(\Delta x)}^2\rangle }}{V_{\text{p}}}\approx \frac{R_{*}}{V_{\text{p}}\sqrt{R_{*}^2-{(\beta R_{*})}^2}}{\sigma }_{R\text{p}},$$(6)

where ${\sigma }_{R\text{p}}\equiv \sqrt{\langle {(\Delta R_{\text{p}})}^2\rangle }$ is the standard deviation of R_p. We note that, here, the values R_p and its standard deviation σ_Rp do not concern the regular planetary radius, but its “effective” value at the point where the planetary shadow has its first or last contact with the stellar limb (the point “c” in Fig. 5). This “effective radius” may depend on the orientation of the line of sight with respect to the orbit of the planet. Even for a single planet with constant shape, σ_Rp can depend on the impact parameter β.

Taking into account the above explanations, we differentiate Eq. (6) with respect to the parameter β: $$\frac{\partial {\sigma }_{\text{s},\text{e}}}{\partial \beta }\approx \beta \frac{R_{*}^3}{x^3{V}_{\text{p}}}{\sigma }_{R\text{p}}+\frac{R_{*}}{x{V}_{\text{p}}}\frac{\partial {\sigma }_{R\text{p}}}{\partial \beta }.$$(7)

The first term in Eq. (7) is always positive and at β ≈1 (i.e., x → 0) it reaches its maximum, which, according to Eq. (6), also means the achieved maximal value of σ_s,e. Another maximum of σ_s,e is possible due to the second term. Its sign is controlled by the derivative ∂σ_Rp∕∂β. Furthermore, this derivative characterizes the shape of the transiting body because the parameter β is related with the position angle ∠MTS ≡ γ = arcsin[βR_*∕(R_* + R_p)] of the first (or last) contact point “c” in Fig. 5.

Figure 6a shows a histogram of β estimates for all KOIs from NASA EA. The narrow peak at β ≈0 means a likely biased coplanarity of KOI orbits with a distant observer. Therefore, when preparing the histograms in Figs. 6b and c for our set of targets (Table A.1), we exclude objects with β ≲0.07 from consideration, focusing on more realistic random distribution over β. In particular, the histogram in Fig. 6b constructed for the KOIs with the negative border timing cross-correlation r_se < 0 shows a clear concentration of such objects especially those with the confirmed planet status marked by red and green colors at low values of β, whereas thehistogram in Fig. 6c made for the objects with positive cross-correlation r_se > 0 has a maximum at β≈ 1. The cases with r_se > 0 correspond to the TLC shifting as a whole, which is known as the transit timing variability (TTV in Holczer et al. 2016). This effect, likely connected with the celestial-mechanic perturbations, complicates the analysis by adding a new term in Eq. (1). This complicated case is beyond the scope of our study. At the same time, the cases with negative correlation r_se < 0, corresponding to the negligible TTV, admit further simple diagnostics, which is demonstrated below.

We interpret the histogram maximum in Fig. 6b as a manifestation of increased σ_s,e related with the negative correlation r_se. The latter was a criterion for the inclusion of a case in the histogram. In particular, the histogram maximum (maximal value of n) in Fig. 6b indicates the maximal σ_s,e at β ≈0. This fact contradicts the hypothesis of a varying β and therefore indicates that the neglecting by the variations of R_p, done above,is an inconsistent assumption. Therefore, the effect of varying R_p (i.e., second term in Eq. (7)) dominates. Hence, the corresponding derivatives ∂σ_s,e∕∂β and ∂σ_Rp∕∂β are negative in Eq. (7). Moreover, the negative ∂σ_s,e∕∂β is directly seen in Fig. 7, prepared for the most reliable cases with the standard error of border timing estimates ε_s,e < 0.0007 days, found in terms of the coefficient errors of the approximating polynomials and averaged over all time-windows. We note that only objects with r_se < 0 and β > 0.07 (like in Fig. 6b) were included in Fig. 7. In particular, Figs. 7c and d is constructed only for confirmed planets. The fact that Figs. 7a and b and Figs.7c and d are almost the same suggests that the possible non-planetary objects have a negligible influence on the overall result.

The physical meaning of the shown gradients ∂σ_s,e∕∂β can be demonstrated using transformations β → γ and σ_s,e → σ_Rp, where γ is the positional angle as defined above (see Fig. 5) and ${\sigma }_{R\text{p}}=\sqrt{\langle {(\Delta R_{\text{p}})}^2\rangle }$ is the standard deviation of the local radius of the planet at the start/end point of its ingress/egress on the stellar disk (i.e., point “c” in Fig. 5). The radius fluctuation ΔR_p can be found from the triangle ΔTMS with the side TM = x in Fig. 5. As $x=\sqrt{{(R_{*}+R_{\text{p}})}^2-{(\beta R_{*})}^2}$, one can express $R_{\text{p}}=\sqrt{x^2+{(\beta R_{*})}^2}-R_{*}$ and find its fluctuation, which corresponds to the fluctuation Δx, by differentiating $\Delta R_{\text{p}}=(\partial R_{\text{p}}/\partial x)\Delta x=x\Delta x/\sqrt{x^2+{(\beta R_{*})}^2}$. Assuming Δx ≪ x, one can transform the fluctuations ΔR_p and Δx into the standard deviations ${\sigma }_{R\text{p}}=\sqrt{\langle {(\Delta R_{\text{p}})}^2\rangle }$ and ${\sigma }_{\text{s},\text{e}}=\sqrt{\langle {(\Delta x)}^2\rangle }/V_{\text{p}}$, respectively.Taking into account the fact that the squared derivative can be expressed as ${(\partial R_{\text{p}}/\partial x)}^2=x^2/[x^2+{(\beta R_{*})}^2]=1/[1+{(\tan \gamma )}^2]={(\cos \gamma )}^2$, one can write σ_Rp = σ_s,eV_p cosγ. To approximate the true value of the standard deviation of estimates of the transit borders’ timing, we introduce its corrected version ${\sigma }_{\text{s},\text{e}}^c=\sqrt{{\sigma }_{\text{s},\text{e}}^2-{\epsilon }_{\text{s},\text{e}}^2}$, which takes into account the calculation error ε_s,e found as an error of the coefficient c_s,e in the approximating polynomials for the transit-borders in TLC (see Sect. 2.1), averaged over all time windows. Using the corrected value ${\sigma }_{\text{s},\text{e}}^c$ enables us minimize the artificial increase of σ_Rp during the grazing transits at large β and related positional angles γ. For the objects with σ_s,e < ε_s,e, the calculation of ${\sigma }_{\text{s},\text{e}}^c$ is impossible, and so zero values are adopted in such cases. As a result, one obtains the relative radius deviation σ_Rp ∕R_p of a variable transiting object: $$\frac{{\sigma }_{R\text{p}}}{R_{\text{p}}}=\frac{V_{\text{p}}\cos \gamma }{R_{\text{p}}}{\sigma }_{\text{s},\text{e}}^c.$$(8)

In Figure 8, which shows σ_Rp∕R_p versus γ for the most reliable cases (i.e., confirmed planets with r_se < 0, β > 0.07, and border timing errors ε_s,e < 0.0007 days), one can see the clusters of increased σ_Rp∕R_p values at low γ. Altogether, Figs. 6b, 7, and 8 show results that support the existence of variable obscuring zones (VOZ) around some KOIs at altitudes up to ~ 0.1 R_p near the orbital plane (γ ≲ 50 deg.). In this respect, it is worth noting that the parameter R_p used above should be understood in the sense of effective local radius of the transiter at the moment of its first (or last) contact with the stellar disk. The linear regressions ⟨σ_Rp∕R_p⟩ = Gγ + Q, were found with G = −0.17 ± 0.08%/deg, Q = 8.51 ± 2.24% for the ingress and G = −0.24 ± 0.08%/deg, Q = 11.46 ± 2.17% for the egress, respectively (see the lines in Fig. 8). Using these regressions one can estimate the borders of VOZ in the stellar disk in the corresponding Cartesian coordinates as follows: $$\frac{X}{R_{\text{p}}}=\left(\sqrt{12}\langle \frac{{\sigma }_{R\text{p}}}{R_{\text{p}}}\rangle +1\right)\cos \gamma ,$$(9) $$\frac{Y}{R_{\text{p}}}=\pm \left(\sqrt{12}\langle \frac{{\sigma }_{R\text{p}}}{R_{\text{p}}}\rangle +1\right)\sin \gamma ,$$(10)

where $\sqrt{12}$ is the ratio of the amplitude of randomized ΔR_p to its standard deviation, assuming a homogenous distribution function.

The X axis is directed along the visible KOI movement and is positive in the ingress and negative in the egress scenarios. Since photometry cannot distinguish between positive and negative values of β, or of γ, both options ± Y naturally contribute to the recovered shape of VOZ (i.e., the range of KOI’s effective radius fluctuation) shown in Fig. 8c. This shape confirms the conclusion on ∂σ_Rp∕∂β < 0 in the above comparative analysis of histograms in Fig. 6. Hence, a typical KOI with r_s,e < 0 shows signs of variable size along its orbital path (i.e., along the X axis).

Although Fig. 8 shows some asymmetry in the VOZ between its ingress and egress parts, this effect is not significant because the maximal difference between $\langle \frac{{\sigma }_{R\text{p}}}{R_{\text{p}}}\rangle =Q$ at γ = 0 as well as the corresponding X∕R_p (see Eq. (9)) for the ingress and egress are approximately equal to their standard errors. Figure 9a, constructed only for the confirmed planets, shows that in general σ_s ≈ σ_e. However, the maximum of the histogram in Fig. 9b is shifted toward higher values of σ_s relative to the corresponding maximum in Fig. 9c. The histogram of difference Δσ ≡ σ_s − σ_e in Fig. 9d has the significantly positive skewness S = 1.81 which is defined as $$S=\frac{\sqrt{m(m-1)}}{m-2}\left\{\frac{\frac{1}{m}{\sum }_{i=1}^m{(\Delta {\sigma }_i-\langle \Delta \sigma \rangle )}^3}{{\left[\frac{1}{m}{\sum }_{i=1}^m{(\Delta {\sigma }_i-\langle \Delta \sigma \rangle )}^2\right]}^{3/2}}\right\},$$(11)

where Δσ_i is an estimate of Δσ for an individual KOI with number i; ⟨Δσ⟩ is an averaged value over all Δσ_i estimates; and m = 73 is the total number of considered objects. We note that non-planets and unconfirmed candidates are ignored here. For the normal distribution, the sample skewness has an expected zero value and a standard error ${\sigma }_o=\sqrt{6m(m-1)/[(m-2)(m+1)(m+3)]}=0.28$ (Kendall & Stuart 1969). The null hypothesis (i.e., S = 0) is rather unlikely, because in our case S∕σ_o = 1.81∕0.28 = 6.5. Hence, there is the slight statistical asymmetry σ_s > σ_e.

Fig. 4 Anti-correlated variations of Δts (red) and Δte (blue) for KOI 18.01 revealed in the TLC of KIC 8191672. The measured timing values are shown as crosses with error bars.

Fig. 5 Transit scheme. Labels “S” and “T” mark the centers of a stellar disk (yellow) and transiter’s disk (blue), respectively.The contact point “c” of the disks corresponds to the transit start (or end) time, i.e., the time moment Δts or Δte. The center of the transiting body moves along the visible trajectory from "T" to “M” (or vise versa), where “M” is the middle point of the trajectory.

Fig. 6 Distributions of the impact parameter β estimates: panel a: over the whole list of KOIs in the NASA Exoplanet Archive. Panel b: for the KOIs from Table A.1 with β > 0.07 and rse< 0. Panel c: for the KOIs from Table A.1 with β > 0.07 and rse> 0. The ordinate value n is the number of β-estimates in a bin of histogram. Similarly to Fig. 3, color shows the status of objects within a bin: violet – probably non-planets (“n” in Table A.1); blue – candidates (“c” in Table A.1); green – confirmed planets, but with unknown masses (“p” in Table A.1); red – confirmed planets with measured planetary masses (“p!”in Table A.1). Cases with “*” and “?” are also included.

Fig. 7 Deviation, σs,e, of transit borders’ timing, versus β estimates, confirming the negative derivative ∂σs,e∕∂β in Eq. (7) for the ingress (panel a) and egress (panel b) parts of the TLC. Only objects with rse < 0, β > 0.07, and border timing errors εs,e < 0.0007 days from Table A.1 were included. The solid lines show the corresponding regressions. Panels c and d are the same as panels a and b, but for the confirmed planets with status “p” only, irrespective of extensions (“!”, “*”, “?”).

Fig. 8 Relative radius deviation σRp∕Rp vs. positional angle γ for the ingress (panel a) and egress (panel b) phases of confirmed planets. The range of effective radius variability of a KOI was calculated using Eq. (8) and the timing standard deviations from Figs. 7c and d. The linesshow the corresponding regressions. Using these regressions, a generalized shape of the VOZ (gray color) around a KOI (black disk) is visualized in panel c.

3.2 Transit shape analysis

Figure 10a shows the distribution of ⟨A_s⟩ estimates for all KOIs from Table A.1. One can see that the histogram has a maximum at ⟨A_s ⟩ = 0.5, which is typical for a symmetric transiter. However, there is the significant skewness of this distribution: S = 3.63 ± 0.24 for all objects in the analyzed data set, and S = 3.48 ± 0.28 for the confirmed planets, i.e., objects marked with green and red in the figure. These values of S are 15.1σ₀ and 12.4σ₀ respectively, where σ₀ = 0.24 and 0.28 are the corresponding standard errors for the normal distribution with S = 0 (see details in Sect. 3.1). As the found skewness is significant and positive, there is an excess of transiting bodies with a slight TLC asymmetry A_s > 0.5, hence the increased radius of a shadow on its leading limb. Although this asymmetry is significant as a cumulative effect over both considered sets of objects (i.e., all KOIs from Table A.1 and those with the confirmed planet status), individual estimates of the transit shape parameter rarely show significant deviation of A_s from the value 0.5 of a symmetric TLC. Such deviated estimates are shown in Figure 10b. It is remarkable that the KOIs with significant TLC asymmetry |⟨A_s⟩− 0.5| > 3σ_A, where σ_A is the standard deviation of A_s (given in Table A.1), are also associated with the negative cross-correlation r_se. We note that the distribution of the TLC shape parameter A_s for planets with measured masses (marked red in Fig. 10a) has the negative skewness S =−0.93 ± 0.38 related with the distribution tail at A_s < 0.5. This fact could be interpreted in terms of light-obscuring tail-like features of some planets. Such tailed bodies slightly affect overall statistics with positive S.

Figure 11 shows the histograms of the estimated Pearson correlation coefficients between the timing parameters Δt_s, Δt_e, and Δt_m of TLCs and the transit shape parameter A_s. One can see that the transit asymmetry is controlled mainly by Δt_m, which is indicated by the maximum in the histogram in Fig. 11a at r_Am = 0.9 ± 0.1, including 62% of cases. At the same time, the transit borders according to the maxima in histograms in Fig. 11b and c show lower correlations of r_As = −0.5 ± 0.1 and r_Ae = −0.5 ± 0.1.

To interpret Fig. 11, one should consider the influence of small fluctuations of the transit timing δ(Δt_s), δ(Δt_m), and δ(Δt_e) on the variation of the transit shape parameter A_s = (Δt_m − Δt_s)∕(Δt_e − Δt_s): $$\delta A_{\text{s}}\approx \frac{\partial (A_{\text{s}})}{\partial (\Delta t_{\text{s}})}\delta (\Delta t_{\text{s}})+\frac{\partial (A_{\text{s}})}{\partial (\Delta t_{\text{m}})}\delta (\Delta t_{\text{m}})+\frac{\partial (A_{\text{s}})}{\partial (\Delta t_{\text{e}})}\delta (\Delta t_{\text{e}}),$$(12)

where the partial derivatives of A_s are $$\frac{\partial (A_{\text{s}})}{\partial (\Delta t_{\text{s}})}=-\frac{\Delta t_{\text{e}}-\Delta t_{\text{m}}}{{(\Delta t_{\text{e}}-\Delta t_{\text{s}})}^2}\approx -\frac{1}{2D},$$(13) $$\frac{\partial (A_{\text{s}})}{\partial (\Delta t_{\text{m}})}=\frac{\Delta t_{\text{e}}-\Delta t_{\text{s}}}{{(\Delta t_{\text{e}}-\Delta t_{\text{s}})}^2}\approx \frac{1}{D},$$(14) $$\frac{\partial (A_{\text{s}})}{\partial (\Delta t_{\text{e}})}=-\frac{\Delta t_{\text{m}}-\Delta t_{\text{s}}}{{(\Delta t_{\text{e}}-\Delta t_{\text{s}})}^2}\approx -\frac{1}{2D}.$$(15)

In Eqs. (13)–(15) we introduce the notation D ≡ Δt_e − Δt_s ≈ Δt_tr + δt_L and express Δt_m − Δt_s ≈ Δt_e − Δt_m ≈ 0.5D assuming a quasi-symmetric TLC. Correspondingly, the correlation coefficients are $$r_{\text{As}}=\frac{\langle \delta A_{\text{s}}\delta (\Delta t_{\text{s}})\rangle }{{\sigma }_{\text{A}}{\sigma }_{\text{s}}}\approx -\frac{{\sigma }_{\text{s}}}{2D{\sigma }_{\text{A}}}+\frac{r_{\text{ms}}{\sigma }_{\text{m}}}{D{\sigma }_{\text{A}}}-\frac{r_{\text{se}}{\sigma }_{\text{e}}}{2D{\sigma }_{\text{A}}},$$(16) $$r_{\text{Am}}=\frac{\langle \delta A_{\text{s}}\delta (\Delta t_{\text{m}})\rangle }{{\sigma }_{\text{A}}{\sigma }_{\text{m}}}\approx -\frac{r_{\text{ms}}{\sigma }_{\text{s}}}{2D{\sigma }_{\text{A}}}+\frac{{\sigma }_{\text{m}}}{D{\sigma }_{\text{A}}}-\frac{r_{\text{me}}{\sigma }_{\text{e}}}{2D{\sigma }_{\text{A}}},$$(17) $$r_{\text{Ae}}=\frac{\langle \delta A_{\text{s}}\delta (\Delta t_{\text{e}})\rangle }{{\sigma }_{\text{A}}{\sigma }_{\text{e}}}\approx -\frac{r_{\text{se}}{\sigma }_{\text{s}}}{2D{\sigma }_{\text{A}}}+\frac{r_{\text{me}}{\sigma }_{\text{m}}}{D{\sigma }_{\text{A}}}-\frac{{\sigma }_{\text{e}}}{2D{\sigma }_{\text{A}}},$$(18)

where ${\sigma }_{\text{A}}\equiv \sqrt{\langle \delta A_{\text{s}}^2\rangle }$, σ_s, σ_m, and σ_e are the standard deviations of A_s, Δt_s, Δt_m, and Δt_e, respectively.

Let us begin with a hypothesis that the timing changes are only due to statistical error without any connection some physical effect or process. In this case one might expect uncorrelated fluctuations of timing (i.e., r_se = r_ms = r_me = 0). Correspondingly, the shape-parameter fluctuation in Eq. (12) can be transformed into a standard deviation: $${\sigma }_{\text{A}}=\sqrt{\langle {(\delta A_{\text{s}})}^2\rangle }=\frac{1}{D}\sqrt{\frac{{\sigma }_{\text{s}}^2}{4}+{\sigma }_{\text{m}}^2+\frac{{\sigma }_{\text{e}}^2}{4}}.$$(19)

Substituting Eq. (19) into Eq. (17) with r_ms= r_me = 0, one obtains $$r_{\text{Am}}=\frac{1}{\sqrt{1+{\left(\frac{{\sigma }_{\text{s}}}{2{\sigma }_{\text{m}}}\right)}^2+{\left(\frac{{\sigma }_{\text{e}}}{2{\sigma }_{\text{m}}}\right)}^2}}.$$(20)

Equation (20) shows that r_Am is always positive (0 < r_Am < 1). The most probable correlation r_Am > 0.8 in Fig. 11a is possibly provided by the following condition: $$\frac{\sqrt{{\sigma }_{\text{s}}^2+{\sigma }_{\text{e}}^2}}{{\sigma }_{\text{m}}}<2\sqrt{\frac{1}{{(0.8)}^2}-1}=\mathrm{1.5.}$$(21)

As a quasi-flat minimum of TLC is harder to localize than the border timing, the realistic condition σ_s,e ≲ σ_m is reasonable. However, an obviously unrealistic condition σ_s,e ≫ σ_m is needed for r_Am ≈ 0. Nevertheless, one can see that the histogram in Fig. 11a continues even further, in the region r_Am ≲ 0, which contradicts Eq. (20).

Using Eqs. (16) and (19) for r_ms= r_se = 0, one can derive the expression for r_As $$r_{\text{As}}=-\frac{1}{\sqrt{1+{\left(2\frac{{\sigma }_{\text{m}}}{{\sigma }_{\text{s}}}\right)}^2+{\left(\frac{{\sigma }_{\text{e}}}{{\sigma }_{\text{s}}}\right)}^2}}.$$(22)

Equation (22) shows that r_As is always negative. Nevertheless, the histogram in Fig. 11b has a continuation in the positive region of r_As> 0, which is true for 19% of the considered KOIs. The analogous continuation of the histogram of r_Ae in the region of positive values r_Ae > 0 is even moreprominent and takes place for 35% of cases in Fig. 11c, although Eq. (18) predicts negative r_Ae for r_se = r_me = 0, $$r_{\text{Ae}}=-\frac{1}{\sqrt{1+{\left(\frac{{\sigma }_{\text{s}}}{{\sigma }_{\text{e}}}\right)}^2+{\left(2\frac{{\sigma }_{\text{m}}}{{\sigma }_{\text{e}}}\right)}^2}}.$$(23)

Now, let us test the above assumption of r_ms = r_me = 0 with the diagram r_ms versus r_me. In Fig. 12 one can see a clustering of most estimate points around r_ms = r_me = 0 in accordance with the typical error bars. However, there are highly deviating points and even an outlying cluster at r_ms < −0.5 and r_me > 0.5. This fact together with the values r_As > 0 and r_Ae > 0 forbidden for r_ms = r_me = 0 leads us to suspect that the uncorrelated timing changes are not the only factor contributing to r_ms and r_me.

The most reasonable physical interpretation of the dominating influence of Δt_m on A_s is the strong effect of starspots which may slightly shift Δt_m in the central part of the TLC. However, this effect should be weakened towards the transit borders due to the decrease of the visible area of a spot near the stellar limb (projection effect). Correspondingly, one can expect r_ms≈ 0 and r_me≈ 0. In the cases of a dominating influence of spots (i.e., σ_m ≫ σ_s,e, hence, statistically |δ(Δt_m)|≫|δ(Δt_s,e))|), it follows from Eqs. (12) and (14) that δA_s ≈ δ(Δt_m)∕D. In this case, the variance of A_s is $\langle \delta A_{\text{s}}^2\rangle =\langle \delta {(\Delta t_{\text{m}})}^2\rangle /D^2$, and the standard deviation is ${\sigma }_{\text{A}}={(\langle \delta A_{\text{s}}^2\rangle )}^{1/2}={\sigma }_{\text{m}}/D$. Substitution of the latter in Eq. (17) with r_ms = r_me = 0 gives r_Am = 1, which agrees with the histogram in Fig. 11a. The cases with significantly non-zero r_ms and r_me in Fig. 12 are not a problem in cases with strong spot effect, i.e., σ_s,e≪ σ_m, meaning that the terms containing r_ms and r_me in Eq. (17) are negligible. Hence, the prediction r_Am ≈ 1 is still valid for any values of r_ms and r_me in the approximation of a strong spot effect.

For a moderate spot effect (σ_m ~ σ_s,e), the first and third terms in Eq. (17) could give significantly negative deposits in r_Am when r_ms> 0 and r_me> 0, respectively. Such an effect is possible in the case of TTVs (Holczer et al. 2016). When the TTV effect shifts the TLC as a whole, the positive correlations r_ms> 0 and r_me> 0 take place. The corresponding decrease of r_Am could explain the fact that many objects (38%) in Fig. 11a show a lower correlation r_Am < 0.8 (up to negative values).

The positive correlations r_As and r_Ae, forbidden for the uncorrelated border timing fluctuations, may still take place due to statistical errors in the correlation coefficient estimates. The scale of the leakage of these estimates into regions r_As> 0 and r_Ae> 0 is of the same order as the statistical errors ${\sigma }_{\text{rAs},\text{e}}=(1-r_{\text{As},\text{e}}^2)/\sqrt{N_{\text{w}}-1}$ of Pearson coefficients r_As or r_Ae for the correlations between the TLC shape parameter A_s and the border timing Δt_s or Δt_e. As the number N_w of the used time-windows is the same for r_As and r_Ae for each particular KOI, and the histograms of Figs. 11b,c are similar for negative correlations, one can predict the statistical proximity σ_rAs ≈ σ_rAe and the similarity between these histograms in the regions of positive correlations r_As > 0 and r_Ae > 0. Nevertheless, there is a noticeable asymmetry. Specifically, in the histogram of r_Ae (Fig. 11c), 35% of the considered objects have r_Ae > 0, while in the histogram of r_As (Fig. 11b) only 19% of the objects show a positive correlation r_As > 0. Therefore, there should be an additional source of positive correlation besides statistical errors.

The uncorrelated estimates of Δt_s,m,e as well as the spot effect justify the guess that r_ms = r_me = 0, making negligible the related terms in Eqs. (16) and (18), which could lead to positive deposits in r_As and r_Ae. Another source of the positive correlations r_As > 0 and r_Ae > 0 are the terms containing r_se in the same equations. If r_se < 0, as was demonstrated in Sect. 3.1, the third term in Eq. (16) and the first term in Eq. (18) give the positive deposits in r_As and r_Ae, respectively. The more pronounced (35% of cases) continuation of the histogram to the region r_Ae > 0 in Fig. 11c, as compared with 19% of cases with r_As > 0 in Fig. 11b, means that there is a statistical inequality σ_s > σ_e in the symmetric terms − r_seσ_s∕(2Dσ_A) and − r_seσ_e∕(2Dσ_A). This results in the unequal positive deposits in r_Ae (see Eq. (18)) and in r_As (see Eq. (16)), which is consistent with the conclusions in Sect. 3.1 on detected statistical inequalities σ_s > σ_e and r_se < 0 in many cases.

Fig. 9 Comparison of standard deviations of the transit border timing: panel a: σs vs. σe. Panel b: distribution histogram of σs. Panel c: distribution histogram of σe. Panel d: distribution histogram of σs − σe. Only confirmed planets from Table A.1 were considered.

Fig. 10 Average TLC shape parameter ⟨As⟩ for all KOIs from Table A.1: panel a: distribution histogram of ⟨As ⟩. Panel b: ⟨As⟩ versus rse for significantly non-symmetric TLCs with |⟨As⟩− 0.5| > 3σA. The KOI numbers are labeled (no non-stellar objects among them). The color coding of the KOI status is similar to one used in Figs. 3 and 6.

Fig. 11 Pearson correlation coefficients between the transit shape parameter As and TLC timing parameters: panel a: distribution histogram of rAm – the correlation between As and Δtm. Panel b: distribution histogram of rAs – the correlation between As and Δts. Panel c: distribution histogram of rAe - the correlation between As and Δte. Color showsthe status of objects within a bin, as in Figs. 3 and 6.

Fig. 12 Distribution of rms vs. rme according Table A.1. Color shows the status of objects, as in Fig. 3.

Fig. 13 Quasi-sinusoidal oscillations of Δts in TLCs of KOI 840.01 at KIC 5651104 (panel a) and KOI 908.01 at KIC 8255887 (panel b). For comparison, panels c and d show the chaotic behavior of Δte by the sameobjects.

3.3 Individual peculiarities

In the case of non-varying transits, the estimates of timing parameters Δt_s, Δt_e, and Δt_m are dispersed irregularly within the range of about plus or minus one standard error (Figs. 2a, b, and d). However, TLCs of KOI 840.01 and KOI 908.01, orbiting KIC 5651104 and KIC 8255887, respectively, show quasi-sinusoidal oscillations of Δt_s (Figs. 13a and b), whereas the behavior of Δt_e appears chaotic (Figs. 13c and d). Correspondingly, the cross-correlation r_se in Table A.1 is insignificant: 0.37 ± 0.31 (for KOI 840.01) and 0.37 ± 0.31 (for KOI 908.01).

Moreover, in some cases (see Fig. 14) there is a clear transit shape asymmetry (i.e., A_s ≠ 0.5), which can also vary. The known case of KOI 13.01 (KIC 9941662) in Fig. 14a has been interpreted as an example of gravity darkening effect appearing during a tilted transit in front of a fast rotating (the period is 1.06 days) star (Szabó et al. 2012 and therein). At the same time, it is worth mentioning that besides the typical gravity darkening phenomenology, the object KOI 13.01 also demonstrates a well-pronounced anti-correlation r_se = −0.55 ± 0.14 at the level beyond three standard errors (r_se∕ε_se = 3.93, see in Table A.2). This attribute constitutes an additional enigma, which cannot be interpreted as merely a darkening effect and therefore suggests an additional manifestation of VOZ. Moreover, gravitational darkening is not valid at all in the case of KOI 3.01 (at KIC 10748390) in Fig. 14b, because its host star is a slow rotator with a period of ≈ 30 days (see Fig. 15).

Figure 16 shows clear cyclic variation of the visible transit duration D = Δt_e − Δt_s, as well as anti-correlation r_se = −0.57 ± 0.11 between Δt_s and Δt_e for the KOI 971.01 orbiting KIC 11180361. This system exhibits grazing transits, where less than half of the transiting shadow is projected onto the stellar disk. Such geometry, in combination with an extremely short transit period of P_tr = 0.533 days, is verysensitive to variations in the size of the transiter. Based on the geometry treatment in Sect. 3.1, it can be shown that $$\frac{1}{D}\frac{\partial D}{\partial z}=\frac{1+z}{{(1+z)}^2-{\beta }^2}=13.0,$$(24)

where the planetary-to-stellar radius ratio z = R_p∕R_* = 0.302 and the impact parameter β = 1.263 are defined with the cumulative data from NASA EA. Therefore, in order to provide the measured ± 11% of the relative variations of D (as seen in Fig. 16a), the ratio R_p∕R_* according to Eq. (24) has to change within the range of ~ ± 11%∕13 = ±0.8%. We note that the star KIC 11180361 shows δ Sct-type pulsations (Balona 2016), the period of which is nevertheless too short (~ 0.1 day) to be related with the cycles in Fig. 16a of ~ 500 days. Although this system was considered as an eclipsing star (Slawson et al. 2011), no variations of radial velocity were detected, nor were signs of spectral contamination from the companion (Lampens et al. 2018). The current status of the KIC 11180361/KOI 971.01 in NASA EA is a candidate. Even if KIC 11180361/KOI 971.01 is a star, in light of its rather unique features, this object deserves a dedicated study, which we attempt to initiate here.

Fig. 14 Transit shape parameter As of KOI 13.01 orbiting KIC 9941662 (panel a) and KOI 3.01 at KIC 10748390 (panel b). For comparison, panels c and d show the synchronous variations of Δtm by the sameobjects.

Fig. 15 Rotational modulation with a 30 day period in the light curve of KOI 3.01 hosted by KIC 10748390. The long- and short-cadence data are marked by the red and crimson colors, respectively.

4 Discussion

The host stars of all the KOIs considered here with varying transit border timings (Table A.2) are main sequence objects. This can be seen in an analog of the HR-diagram in Fig. 17 prepared for these objects, i.e., the logarithm of gravity log (g) versus effective temperature T_eff in the stellar photospheres. A cluster of stars around the solar position can be clearly seen there. For example, KIC 12019440 has almost the same parameters as the Sun. The numerous attempts (Qu et al. 2015 and therein) to detect variability in the solar radius on long timescales from 1 day to decades give a negligible upper limit of the effect at the level of < 0.5 arcseconds or <0.05%. At the same time, according to the analysis in Sect. 3.1, the distance TS = R_p + R_* in Fig. 5, which defines the transit border timing variations, changes much more than just only the varying stellar radius might contribute. In particular, the regressions in Figs. 8a and b reveal the variations of R_p with the scale of σ_Rp ∕R_p ≈ 10%, which yields an estimate of σ_Rp∕R_* = (σ_Rp∕R_p)(R_p∕R_*) ~ 1% for a typical hot jupiter ratio R_p∕R_*~ 0.1. We note that the sun-like stars show variability in their radii with an order of magnitude smaller amplitude than the observed scale of variations of TS = R_p + R_*. Therefore, the main source of TS variability is connected with the varying R_p rather than with R_s. This conclusion supports theconstant R_s assumption,which is applied above during the derivation and analysis of Eqs. (1)–(10).

Let us verify whether or not the cyclic variability of Δt_s, Δt_e, and Δt_m found above and the corresponding A_s could indicate the presence of precessing exo-rings or disks. To model these cyclic phenomena, we consider a transiting spherical planet orbiting a real star KIC 11359879 (host of Kepler-15b) with a planetocentric opaque disk tilted at 30 degrees relative the orbital plane and precessing with a period of 2400 days (see in Fig. 18a). For the reconstruction of the corresponding TLSs, we use a pixel-by-pixel integration, which can be used for any type or geometry of transiter. The dimming of stellar flux during the transit is characterized by the part of starlight blocked by the transiting object $$\Delta F=\frac{{\displaystyle \int {\displaystyle \int I}}(x,y)\text{d}x\text{d}y}{{\displaystyle \int {\displaystyle \int I_{\text{s}}}}(x,y)\text{d}x\text{d}y},$$(25)

where we use a coordinate system co-centered on the stellar disk, with the x-axis along the planet orbit projection onto the stellar disk. By this, I_s is radiation intensity at a given position (x,y) on the visible stellar disk, and I is the same intensity but disturbed by the transiter. The integrals in Eq. (25) can be replaced by sums over N_p pixels with serial number i: $$\Delta F=\frac{{\sum }_{i=1}^{N_{\text{p}}}I(x_i,y_i)}{{\sum }_{i=1}^{N_{\text{p}}}{I}_{\text{s}}(x_i,y_i)}.$$(26)

The stellar limb darkening is taken into account according to the best (four coefficients) approximation by Claret & Bloemen (2011), depending on the particular stellar effective temperature and gravity adopted from NASA EA. The planetary data (radius R_p, semi-major axis a_p of the orbit, impact parameter β, mid-time t₀ of the first observed transit, and transit period P_tr) are also taken from the NASA Exoplanet Archive. In the reference system used here, the moving center of the exoplanetary shadow has coordinates $$\begin{array}{lll}{x}_{\text{p}}\hfill & =\hfill & a_{\text{p}}\sin \left[\frac{2\pi }{P_{\text{tr}}}(t-t_0)\right],\hfill \\ *y_{\text{p}}\hfill & =\hfill & \beta R_{*},\hfill \end{array}$$

with the value for stellar radius R_* taken from NASA EA. Using this approach, we calculate a synthetic light curve, which is processed analogously to the real photometricdata. Figure 18 summarizes the obtained results.

In particular, Fig. 18 reveals a clear variability in A_s, Δt_s, and Δt_e, but an almost complete lack of noticeablevariability in Δt_m. These results, to a certain degree, resemble the measured behavior of the TLC border timing and shape parameters of real objects shown in Figs. 13 and 14b and confirm that although the disk effectively influences the transit borders, it cannot noticeably displace Δt_m because of the transiter’s symmetry. However, as shown in the analysis in Sect. 3.2, the varying minimum time Δt_m of the TLCis the main driver of the variations of its shape parameter A_s. This, in particular, results in the cross-correlation r_Am values closeto unity obtained for the real objects with oscillating Δt_s, for example, r_Am = 0.98 ± 0.02 and 0.94 ± 0.05 for KIC 8255887 and KIC 10748390 respectively. Altogether, this discrepancy, with regard to the relation between Δt_m and A_s as predicted bythe model and revealed from observations, rules out the hypothesis of a precessing circumplanetary disk as the main driver of thevarying TLC shape parameter, whereas the variation of only border timing in some cases could still be connected with such a disk.

As the TTV effect cannot influence the shape parameter A_s, as explained in Sect. 3.2, the measured (in some cases) phenomenon of a variable A_s remains to be related with the effect of starspots. If this is the case, then the regular variability of Δt_m and correspondingly A_s seen in Figs. 14b and d suggests some long-living (~10³ days) starspot or activity complex in the stellar disk of KIC 10748390, the position of which appears to have a certain relation to the transiting exoplanet KOI 3.01. In other words, in order to have a quasi-stable, that is, detectable Δt_m≠0 within one time-window, which at the same time can vary slightly from one time-window to another, the starspot(s) must occupy approximately the same position on the stellar disk during different transits in the window. Since in the considered case of KIC 10748390 the stellar rotation period is much longer (≈ 30 days according to Fig. 15) than the transiter’s orbital period (P_tr = 1.24 days), the above-mentioned condition could be fulfilled if the starspot appearance is related with the planet. This conclusion does not seem impossible in light of recent reports on planetary related activity regions in host stars (e.g., Cauley et al. 2019 and therein). An alternative scenario with a transit over a polar starspot requires fast stellar rotation (Yadav et al. 2015) which is not the case in the stellar system considered here.

Fig. 16 Temporal variation of panel a: transit duration Δte − Δts; panel b: transit start-time Δts, and panel c: transit end-time Δte in the extreme case of grazing transits of KOI 971.01 orbiting KIC 11180361.

Fig. 17 Stellar gravity (log (g)) vs. effective temperature (Teff) distribution of the systems with the noticeably varying transit borders (see Table A.2) according to the NASA EA. The red cross marks the position of the Sun.

5 Conclusions

The results reported in this paper lead to several conclusions which can be summarized as follows.

• 1.

  The transit border timings of some KOIs (exclusively hot jupiters in Table A.2) are variable on timescales from ≈ 400 to ≲ 1500 days (see Figs. 4, 13 and 16).

• 2.

  Among the most typical features of the TLC timing variability is the significant anti-correlation r_se < 0 between the transit border timing parameters Δt_s and Δt_e (see Table A.2 as well as Figs. 4, 16). This anti-correlation is likely a manifestation of the variability of the dimensions of KOIs. A hypothetical variability of the impact parameter or the stellar radius, assumed as an alternative mechanism for the anti-correlated TLC border timings were shown to contradict the revealed facts (see the argumentation in Sects. 3.1 and 4).

• 3.

  The range of variability in the dimensions of KOIs, introduced as a so-called varying obscuring zone (VOZ), extends mainly along the orbit, up to a maximum of ~10% of the transiter’s size on average, and disappears in the perpendicular direction as shown in Fig. 8. This feature suggests generation of dust or aerosol in the upper layers of the atmospheres of KOIs and above, especially in the equatorial regions (i.e.,close to the orbital plane), which are the subject of higher stellar radiative impact as compared to the poles.

• 4.

  The aerosol or dust particle clouds in the VOZ, detected in the present study, may also be connected with the detected dusty obscuring matter (DOM) ahead of hot jupiters at altitudes ~ 2R_p (Arkhypov et al. 2019). This assumption is in particular supported by the fact that the signatures of pre-transit DOM are detected for six KOIs (in Table 1 in Arkhypov et al. 2019) that also have the varying transit borders timing (i.e., appear in Table A.2). These objects are the confirmed exoplanets KOI 13.01, KOI 17.01, KOI 18.01, KOI 20.01, and KOI 186.01, and the candidate exoplanet KOI 6085.01 (according to NASA EA). Indeed, as seen in the histograms in Figs. 10a and 9, there is a perceptible statistical excess of the ahead radius of the KOIs, estimated from the ingress part of the TLC, as compared to the opposite side values, which manifests itself as a significant positive skewness of the distribution of estimates for A_s and statistic inequality σ_s > σ_e (see Sect. 3).

• 5.

  Starspots appear to be the main factor disturbing the generally symmetric TLC shape (i.e., giving A_s ≠ 0.5). At the same time, the starspots usually have a limited lifetime and a quasi-accidental (random) appearance in longitude. Thus, they should result in a disordered, irregular variability of A_s. Nevertheless, the system KIC 10748390 / KOI 3.01 demonstrates the regular variation of A_s on a timescale of ~10³ days (Fig. 14). This object is a candidate for having a starspot, which is associated (synchronized) with the transiter.

• 6.

  Altogether, the independent measurement of the transit border timing performed here opens a way to study yet unexplored phenomena, such as for example the oscillations in Δt_s (see Fig. 13).

In summary, we conclude that the transit border timing is a new and effective but still underused instrument for the in-transit probing of exoplanetary exteriors.

Fig. 18 Modeling of TLC variability for an analogue of Kepler-15b with a prescribed precessing ring: panel a: model image of the transiting exoplanet with a precessing ring (solid line shows the trajectory of planet). Panel b: variations of the shape parameter As of the simulated synthetic TLC. Panel c: variations of the simulated transit duration Δte − Δts. Panel d: variations of the simulated transit depth ΔFtr. Panel e: variations of the simulated TLC minimum time Δtm. Panel f: variations of the simulated transit start-time Δts. Panel g: variations of the simulated transit end-time Δte. The abscissascale ΔJD = JD − 2 454 833.0 is in Julian days (JD).

Acknowledgements

The authors acknowledge the projects I2939-N27 and S11606-N16 of the Austrian Science Fund (FWF) for the support. M.L.K. is grateful also to the grant No. 18-12-00080 of the Russian Science Foundation and acknowledges the project “Study of stars with exoplanets” within the grant No.075-15-2019-1875 from the government of Russian Federation.

Appendix A Tables

References

All Tables

All Figures

	Fig. 1 Example of light-curve processing (the star KIC 11414511). Panel a: approximation of the transit background (red squares) with a sixth-order polynomial Fb(t) (green curve). Panel b: folded profile ΔF(Δt) of 25 adjacent individual transits. Blue curves are the applied polynomial approximations used for measurements of labeled parameters.
In the text

	Fig. 2 Temporal behavior of the transit timing parameters of real TLCs like in Fig. 1b, compiled for a sequence of adjacent time-windows, vs. time in Julian days JD: panel a: transit start-time Δts. Panel b: transit end-time Δte. Panel c: transit duration Δte − Δts. Panel d: transit minimum-time Δtm (maximal flux decrease). Panel e: transit shape parameter As. Panel f: transit flux decrease ΔFtr.
In the text

Fig. 3 Distribution of selected KOI targets (listed in Table A.1) vs. the radii Rp and orbital semimajor axes ap of transiters according to NASA EA3. Color shows the status of objects: violet – probably non-planets (“n” in Table A.1); blue – candidates (“c” in Table A.1); green – confirmed planets but with unknown masses (“p” in Table A.1); red – confirmed planets with measured planetary masses (“p!” in Table A.1). The cases with “*” and “?” are included. The same color coding is used throughout the paper.

In the text

	Fig. 4 Anti-correlated variations of Δts (red) and Δte (blue) for KOI 18.01 revealed in the TLC of KIC 8191672. The measured timing values are shown as crosses with error bars.
In the text

	Fig. 5 Transit scheme. Labels “S” and “T” mark the centers of a stellar disk (yellow) and transiter’s disk (blue), respectively.The contact point “c” of the disks corresponds to the transit start (or end) time, i.e., the time moment Δts or Δte. The center of the transiting body moves along the visible trajectory from "T" to “M” (or vise versa), where “M” is the middle point of the trajectory.
In the text

Fig. 6 Distributions of the impact parameter β estimates: panel a: over the whole list of KOIs in the NASA Exoplanet Archive. Panel b: for the KOIs from Table A.1 with β > 0.07 and rse< 0. Panel c: for the KOIs from Table A.1 with β > 0.07 and rse> 0. The ordinate value n is the number of β-estimates in a bin of histogram. Similarly to Fig. 3, color shows the status of objects within a bin: violet – probably non-planets (“n” in Table A.1); blue – candidates (“c” in Table A.1); green – confirmed planets, but with unknown masses (“p” in Table A.1); red – confirmed planets with measured planetary masses (“p!”in Table A.1). Cases with “*” and “?” are also included.

In the text

Fig. 7 Deviation, σs,e, of transit borders’ timing, versus β estimates, confirming the negative derivative ∂σs,e∕∂β in Eq. (7) for the ingress (panel a) and egress (panel b) parts of the TLC. Only objects with rse < 0, β > 0.07, and border timing errors εs,e < 0.0007 days from Table A.1 were included. The solid lines show the corresponding regressions. Panels c and d are the same as panels a and b, but for the confirmed planets with status “p” only, irrespective of extensions (“!”, “*”, “?”).

In the text

Fig. 8 Relative radius deviation σRp∕Rp vs. positional angle γ for the ingress (panel a) and egress (panel b) phases of confirmed planets. The range of effective radius variability of a KOI was calculated using Eq. (8) and the timing standard deviations from Figs. 7c and d. The linesshow the corresponding regressions. Using these regressions, a generalized shape of the VOZ (gray color) around a KOI (black disk) is visualized in panel c.

In the text

	Fig. 9 Comparison of standard deviations of the transit border timing: panel a: σs vs. σe. Panel b: distribution histogram of σs. Panel c: distribution histogram of σe. Panel d: distribution histogram of σs − σe. Only confirmed planets from Table A.1 were considered.
In the text

	Fig. 10 Average TLC shape parameter ⟨As⟩ for all KOIs from Table A.1: panel a: distribution histogram of ⟨As ⟩. Panel b: ⟨As⟩ versus rse for significantly non-symmetric TLCs with |⟨As⟩− 0.5| > 3σA. The KOI numbers are labeled (no non-stellar objects among them). The color coding of the KOI status is similar to one used in Figs. 3 and 6.
In the text

Fig. 11 Pearson correlation coefficients between the transit shape parameter As and TLC timing parameters: panel a: distribution histogram of rAm – the correlation between As and Δtm. Panel b: distribution histogram of rAs – the correlation between As and Δts. Panel c: distribution histogram of rAe - the correlation between As and Δte. Color showsthe status of objects within a bin, as in Figs. 3 and 6.

In the text

	Fig. 12 Distribution of rms vs. rme according Table A.1. Color shows the status of objects, as in Fig. 3.
In the text

	Fig. 13 Quasi-sinusoidal oscillations of Δts in TLCs of KOI 840.01 at KIC 5651104 (panel a) and KOI 908.01 at KIC 8255887 (panel b). For comparison, panels c and d show the chaotic behavior of Δte by the sameobjects.
In the text

	Fig. 14 Transit shape parameter As of KOI 13.01 orbiting KIC 9941662 (panel a) and KOI 3.01 at KIC 10748390 (panel b). For comparison, panels c and d show the synchronous variations of Δtm by the sameobjects.
In the text

	Fig. 15 Rotational modulation with a 30 day period in the light curve of KOI 3.01 hosted by KIC 10748390. The long- and short-cadence data are marked by the red and crimson colors, respectively.
In the text

	Fig. 16 Temporal variation of panel a: transit duration Δte − Δts; panel b: transit start-time Δts, and panel c: transit end-time Δte in the extreme case of grazing transits of KOI 971.01 orbiting KIC 11180361.
In the text

	Fig. 17 Stellar gravity (log (g)) vs. effective temperature (Teff) distribution of the systems with the noticeably varying transit borders (see Table A.2) according to the NASA EA. The red cross marks the position of the Sun.
In the text

Fig. 18 Modeling of TLC variability for an analogue of Kepler-15b with a prescribed precessing ring: panel a: model image of the transiting exoplanet with a precessing ring (solid line shows the trajectory of planet). Panel b: variations of the shape parameter As of the simulated synthetic TLC. Panel c: variations of the simulated transit duration Δte − Δts. Panel d: variations of the simulated transit depth ΔFtr. Panel e: variations of the simulated TLC minimum time Δtm. Panel f: variations of the simulated transit start-time Δts. Panel g: variations of the simulated transit end-time Δte. The abscissascale ΔJD = JD − 2 454 833.0 is in Julian days (JD).

In the text