Issue 
A&A
Volume 638, June 2020



Article Number  A143  
Number of page(s)  17  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201937303  
Published online  30 June 2020 
Variability of transit light curves of Kepler objects of interest
^{1}
Space Research Institute (IWF), Austrian Academy of Sciences,
Schmiedlstrasse 6,
8042
Graz,
Austria
email: oleksiy.arkhypov@oeaw.ac.at; maxim.khodachenko@oeaw.ac.at
^{2}
Skobeltsyn Institute of Nuclear Physics, Moscow State University,
119992
Moscow,
Russia
^{3}
Institute of Astronomy of the Russian Academy of Sciences,
119017
Moscow,
Russia
^{4}
Institute of Physics, KarlFranzens University of Graz, Universitätsplatz 5,
8010
Graz,
Austria
email: arnold.hanslmeier@unigraz.at
Received:
12
December
2019
Accepted:
26
April
2020
Context. Hitherto, the study of exoplanetary transit timing and duration variability has supposed transit light curves (TLCs) to be symmetric, suggesting a priori a spherical shape for the transiting object, for example, an exoplanet. As a result, the independent positions of transit borders are unknown. However, the borders of TLCs are most sensitive to the presence of exorings and/or dust formations of great interest.
Aims. For the first time we check for a timing variability of independently treated borders of transits of different types of exoplanets.
Methods. Using quadratic approximation for the start, end, and minimum parts of the longcadence TLCs from the Kepler mission archive after their whitening and phase folding, we find the corresponding transit border timings: Δt_{s}, Δt_{e}, respectively, and the TLC minimum time Δt_{m}. These parameters were found separately for each folded TLC constructed in the consequent nonoverlapping timewindows with the respective medium time t_{w}. Temporal and crosscorrelation analysis of the obtained series of Δt_{s}(t_{w}), Δt_{e}(t_{w}), and Δt_{m}(t_{w}) were applied for the detection and diagnostics of variability of transit borders and TLC asymmetry.
Results. Among the considered TLCs of 98 Kepler objects of interest (KOIs), 15 confirmed giant exoplanets and 5 objects with still debatable status (probably nonplanets) show variations in their transit timing parameters at timescales from ≈400 to ≳1500 days. These variations are especially well manifested as an anticorrelation between Δt_{s} and Δt_{e}, indicating variability in the dimensions of transiting shadows, especially along their trajectories. There are also objects with well pronounced oscillations of transit border timing and asymmetry.
Conclusions. The discovered variability of transit timing is important as an indicator of largescale nonstationary processes in the atmospheres of KOIs, as well as dust and aerosol generation in their upper layers and in their close vicinity. These findings highlight the need for a dedicated and detailed study.
Key words: planets and satellites: general / planets and satellites: rings / planetstar interactions
© 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 exorings. (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 exorings, 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 (SanchisOjeda 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^{1}), after presearch 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 longcadence 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 shortcadence 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 .
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 lightcurve fragment. (c) Removal of the transit from the analyzed lightcurve fragment, including the halfexposition margins of ± 0.01 days, based on the published transit duration values (NASA EA). (d) Approximation of the remaining lightcurve fragment without transit with a sixthorder polynomial: F_{b}(t) = at^{6} + bt^{5} + ct^{4} + dt^{3} + et^{2} + 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, lightcurve 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 longcadence data records with the abovementioned 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 timewindow of width Δt_{w} centered at t_{w} into one phasefolded TLC ΔF(Δt). Here, Δt ≡ t − t_{E} is the time counted for each transit with number E = 0, 1, 2, ... relative to its midtime t_{E} = t_{0} + P_{tr}E. The reference transit time t_{0} and the transit period P_{tr} are taken in their cumulative versions from the NASA EA. The used timewindows are defined as a sequence of adjacent frames, so that t_{w} = t_{0} + Δt_{w}J_{w}, where J_{w} = 0, 1, 2 ... is the number of a particular timewindow.
The creation of a folded TLC by combining individual transits within a window is also an iterative process aimed at exclusion of lightcurve artifacts, stellar flares, and starspot eclipses. In order to clean the analyzed TLC of outliers related with the abovementioned effects, we approximate it with a fourthorder polynomial , 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 phasefolded TLCs (i.e., the number of timewindows 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 secondorder polynomials. In particular, the minimumpart of the TLC is approximated as ΔF ≈ F_{a} = gΔt^{2} + hΔt + q within the interval Δt < Δt_{tr}∕4 centered at the folded TLC midtime (Δt = 0). An additional reason for using the secondorder polynomial here is to avoid the uncertainty associated with selection among multiple extremes, which are typical for higherorder polynomial approximations. Moreover, the single minimum at Δt_{m} = −h∕(2g) is considered also as the TLC minimumtime, 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 .
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 secondorder polynomial approximations of the reverse dependence Δt(ΔF) = a_{s,e}ΔF^{2} + 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 endparts. 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 timewindow 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 timewindow 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 timewindows, 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 crosscorrelations 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 lightcurve processing (the star KIC 11414511). Panel a: approximation of the transit background (red squares) with a sixthorder polynomial F_{b}(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 longcadence light curves are available, but with high signaltonoise ratios (S∕N ≳ 1000) according to NASA EA. To avoid the intertransit interference in multiplanetary 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 nonplanets. 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 UlmerMoll 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 socalled superpuff 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_{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^{2} and the SIMBAD Astronomical Database^{3}. The objects for which the previously suggested falsepositive 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 outofeclipse variability” indicating that the transitlike event is most likely caused by an eclipsing binary. “However, selfluminous, 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 (Kepler13b 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 nonstellar 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 nonplanetary status “n” is motivatedby the fact that nonplanetary KOIs need to be studied as well as exoplanets. Furthermore, the falsepositive 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 falsepositive 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 superearths 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 timewindows. 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 timewindows, vs. time in Julian days JD: panel a: transit starttime Δt_{s}. Panel b: transit endtime Δt_{e}. Panel c: transit duration Δt_{e} − Δt_{s}. Panel d: transit minimumtime Δt_{m} (maximal flux decrease). Panel e: transit shape parameter A_{s}. Panel f: transit flux decrease ΔF_{tr}. 
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 crosscorrelations 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 R_{p} and orbital semimajor axes a_{p} of transiters according to NASA EA^{3}. Color shows the status of objects: violet – probably nonplanets (“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 timewindows, and ε_{s} is an error of Δt_{s} measurements averaged over all timewindows.
 (b)
For Δt_{e}, the same criterion as for Δt_{s} is considered, that is, σ_{e}∕ε_{e} > 2.5.
 (c)
For r_{se}, significant crosscorrelation r_{se}∕σ_{rse} > 2.5 is required, where is a standard error of r_{se} according to Hotelling (1953), and N_{w} is the number of timewindows in which the phasefolded 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 anticorrelation r_{se} = −0.73 ± 0.11. All 17 KOIs in Table A.2 show analogous (r_{se} ≤−0.51) negative crosscorrelation. Supposing a stable host star (i.e., with a constant radius), such an anticorrelation 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 semitrajectory TM as , where R_{*} and R_{p}, i.e., the distances cS and cT, are the stellar and transiting planet radii, respectively. This semitrajectory 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: (1)
Let us consider the effect of only the first term in Eq. (1) expressed via standard deviations: . Then Eqs. (1) and (2) give (4)
where 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), (5)
For R_{p} ≪ R_{*}, we obtain from Eq. (5) the standard deviations: (6)
where 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 β: (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 crosscorrelation 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 crosscorrelation 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 celestialmechanic 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 timewindows. 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 nonplanetary 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 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 , one can express and find its fluctuation, which corresponds to the fluctuation Δx, by differentiating . Assuming Δx ≪ x, one can transform the fluctuations ΔR_{p} and Δx into the standard deviations and , respectively.Taking into account the fact that the squared derivative can be expressed as , one can write σ_{Rp} = σ_{s,e}V_{p} cosγ. To approximate the true value of the standard deviation of estimates of the transit borders’ timing, we introduce its corrected version , 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 transitborders in TLC (see Sect. 2.1), averaged over all time windows. Using the corrected value 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 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: (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: (9) (10)
where 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 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 (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 nonplanets and unconfirmed candidates are ignored here. For the normal distribution, the sample skewness has an expected zero value and a standard error (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 Anticorrelated variations of Δt_{s} (red) and Δt_{e} (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 Δt_{s} or Δt_{e}. 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 r_{se}< 0. Panel c: for the KOIs from Table A.1 with β > 0.07 and r_{se}> 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 nonplanets (“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 r_{se} < 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}∕R_{p} 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σ_{0} and 12.4σ_{0} respectively, where σ_{0} = 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 crosscorrelation 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 lightobscuring taillike 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}): (12)
where the partial derivatives of A_{s} are (13) (14) (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 quasisymmetric TLC. Correspondingly, the correlation coefficients are (16) (17) (18)
where , σ_{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 shapeparameter fluctuation in Eq. (12) can be transformed into a standard deviation: (19)
Substituting Eq. (19) into Eq. (17) with r_{ms}= r_{me} = 0, one obtains (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: (21)
As a quasiflat 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} (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, (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 , and the standard deviation is . 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 nonzero 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 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 timewindows 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 ⟨A_{s}⟩ for all KOIs from Table A.1: panel a: distribution histogram of ⟨A_{s} ⟩. Panel b: ⟨A_{s}⟩ versus r_{se} for significantly nonsymmetric TLCs with ⟨A_{s}⟩− 0.5 > 3σ_{A}. The KOI numbers are labeled (no nonstellar 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 A_{s} and TLC timing parameters: panel a: distribution histogram of r_{Am} – the correlation between A_{s} and Δt_{m}. Panel b: distribution histogram of r_{As} – the correlation between A_{s} and Δt_{s}. Panel c: distribution histogram of r_{Ae}  the correlation between A_{s} and Δt_{e}. Color showsthe status of objects within a bin, as in Figs. 3 and 6. 
Fig. 12 Distribution of r_{ms} vs. r_{me} according Table A.1. Color shows the status of objects, as in Fig. 3. 
Fig. 13 Quasisinusoidal oscillations of Δt_{s} 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 Δt_{e} by the sameobjects. 
3.3 Individual peculiarities
In the case of nonvarying 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 quasisinusoidal oscillations of Δt_{s} (Figs. 13a and b), whereas the behavior of Δt_{e} appears chaotic (Figs. 13c and d). Correspondingly, the crosscorrelation 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 wellpronounced anticorrelation 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 anticorrelation 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 (24)
where the planetarytostellar 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 δ Scttype 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 A_{s} 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 Δt_{m} 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 shortcadence 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 HRdiagram 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 sunlike 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 exorings or disks. To model these cyclic phenomena, we consider a transiting spherical planet orbiting a real star KIC 11359879 (host of Kepler15b) 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 pixelbypixel 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 (25)
where we use a coordinate system cocentered on the stellar disk, with the xaxis 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: (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}, semimajor axis a_{p} of the orbit, impact parameter β, midtime t_{0} 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
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 crosscorrelation 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 longliving (~10^{3} 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 quasistable, that is, detectable Δt_{m}≠0 within one timewindow, which at the same time can vary slightly from one timewindow 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 abovementioned 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 Δt_{e} − Δt_{s}; panel b: transit starttime Δt_{s}, and panel c: transit endtime Δt_{e} in the extreme case of grazing transits of KOI 971.01 orbiting KIC 11180361. 
Fig. 17 Stellar gravity (log (g)) vs. effective temperature (T_{eff}) 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 anticorrelation 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 anticorrelation 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 anticorrelated 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 socalled 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 pretransit 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 quasiaccidental (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^{3} 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 intransit probing of exoplanetary exteriors.
Fig. 18 Modeling of TLC variability for an analogue of Kepler15b 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 A_{s} of the simulated synthetic TLC. Panel c: variations of the simulated transit duration Δt_{e} − Δt_{s}. Panel d: variations of the simulated transit depth ΔF_{tr}. Panel e: variations of the simulated TLC minimum time Δt_{m}. Panel f: variations of the simulated transit starttime Δt_{s}. Panel g: variations of the simulated transit endtime Δt_{e}. The abscissascale ΔJD = JD − 2 454 833.0 is in Julian days (JD). 
Acknowledgements
The authors acknowledge the projects I2939N27 and S11606N16 of the Austrian Science Fund (FWF) for the support. M.L.K. is grateful also to the grant No. 181200080 of the Russian Science Foundation and acknowledges the project “Study of stars with exoplanets” within the grant No.0751520191875 from the government of Russian Federation.
Appendix A Tables
Analyzed target set and main processing results.
Best examples of KOIs with a varying transit borders timing.
References
 Aizawa, M., Masuda, K., Kawahara, H., & Suto, Y. 2018, AJ, 155, 206 [NASA ADS] [CrossRef] [Google Scholar]
 Almenara, J. M., Damiani, C., Bouchy, F., et al. 2015, A&A, 575, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arkhypov, O. V., Khodachenko, M. L., & Hanslmeier, A., 2019, A&A, 631, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balona, L. A. 2016, MNRAS, 459, 1097 [CrossRef] [Google Scholar]
 Cauley, P. W., Shkolnik, E. L., Llama, J., & Lanza, A. F. 2019, Nat. Astron., 3, 1128 [CrossRef] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2013, ApJ, 772, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Garai, Z. 2018, A&A, 611, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Holczer, T., Mazeh, T., Nachmani, G., et al. 2016, ApJS, 225, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Hotelling, H. 1953, J. R. Stat. Soc. Ser B, 15, 193 [Google Scholar]
 Jackson, B., Adams, E., Sandidge, W., Kreyche, S., & Briggs, J. 2019, AJ, 157, 239 [CrossRef] [Google Scholar]
 Jenkins, J. M., Douglas, A.C., Chandrasekaran, H., et al. 2010, ApJ, 713, L87 [NASA ADS] [CrossRef] [Google Scholar]
 Kendall, M. G., & Stuart, A. 1969, The Advanced Theory of Statistics (Houston: Griffin), 1 [Google Scholar]
 Kepler Mission Team 2009, VizieR Online Data Catalog: V/133 [Google Scholar]
 Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., & Prokopov, P. A. 2015, ApJ, 813, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Lampens, P., Frémat, Y., Vermeylen, L., et al. 2018, A&A, 610, A17 [CrossRef] [EDP Sciences] [Google Scholar]
 Llama, J., Vidotto, A. A., Jardine, M., et al. 2013, MNRAS, 436, 2179 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, T. D., Bryson, S. T., Coughlin, et al. 2016, ApJ, 822, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Pass, E. K., Cowan, N. B., Cubillos, P. E., & Sklar, J. G. 2019, MNRAS, 489, 941 [CrossRef] [Google Scholar]
 Qu, Z. N., Kong, D. F., Xiang, N. B., & Feng, W. 2015, ApJ, 798, 113 [CrossRef] [Google Scholar]
 Rappaport, S., Gary, B. L., Kaye, T., et al. 2016, MNRAS, 458, 3904 [NASA ADS] [CrossRef] [Google Scholar]
 SanchisOjeda, R., Rappaport, S., Pallè, E., et al. 2015, ApJ, 812, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Schlichting, H. E., & Chang, Ph. 2011, ApJ, 734, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Slawson, R. W., Prsa, A., Welsh, W. F., et al. 2011, AJ, 142, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Szabó, G. M., Pál, A., Derekas, A., et al. 2012, MNRAS, 421, L122 [NASA ADS] [CrossRef] [Google Scholar]
 UlmerMoll, S., Santos, N. C., Figueira, P., Brinchmann, J., & Faria, J. P. 2019, A&A, 630, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, L, & Dai, F. 2019, ApJ, 873, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, R. K., Gastine, T., Christensen, U. R., & Reiners, A. 2015, A&A, 573, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Example of lightcurve processing (the star KIC 11414511). Panel a: approximation of the transit background (red squares) with a sixthorder polynomial F_{b}(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 timewindows, vs. time in Julian days JD: panel a: transit starttime Δt_{s}. Panel b: transit endtime Δt_{e}. Panel c: transit duration Δt_{e} − Δt_{s}. Panel d: transit minimumtime Δt_{m} (maximal flux decrease). Panel e: transit shape parameter A_{s}. Panel f: transit flux decrease ΔF_{tr}. 

In the text 
Fig. 3 Distribution of selected KOI targets (listed in Table A.1) vs. the radii R_{p} and orbital semimajor axes a_{p} of transiters according to NASA EA^{3}. Color shows the status of objects: violet – probably nonplanets (“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 Anticorrelated variations of Δt_{s} (red) and Δt_{e} (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 Δt_{s} or Δt_{e}. 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 r_{se}< 0. Panel c: for the KOIs from Table A.1 with β > 0.07 and r_{se}> 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 nonplanets (“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 r_{se} < 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}∕R_{p} 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 ⟨A_{s}⟩ for all KOIs from Table A.1: panel a: distribution histogram of ⟨A_{s} ⟩. Panel b: ⟨A_{s}⟩ versus r_{se} for significantly nonsymmetric TLCs with ⟨A_{s}⟩− 0.5 > 3σ_{A}. The KOI numbers are labeled (no nonstellar 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 A_{s} and TLC timing parameters: panel a: distribution histogram of r_{Am} – the correlation between A_{s} and Δt_{m}. Panel b: distribution histogram of r_{As} – the correlation between A_{s} and Δt_{s}. Panel c: distribution histogram of r_{Ae}  the correlation between A_{s} and Δt_{e}. Color showsthe status of objects within a bin, as in Figs. 3 and 6. 

In the text 
Fig. 12 Distribution of r_{ms} vs. r_{me} according Table A.1. Color shows the status of objects, as in Fig. 3. 

In the text 
Fig. 13 Quasisinusoidal oscillations of Δt_{s} 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 Δt_{e} by the sameobjects. 

In the text 
Fig. 14 Transit shape parameter A_{s} 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 Δt_{m} 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 shortcadence data are marked by the red and crimson colors, respectively. 

In the text 
Fig. 16 Temporal variation of panel a: transit duration Δt_{e} − Δt_{s}; panel b: transit starttime Δt_{s}, and panel c: transit endtime Δt_{e} 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 (T_{eff}) 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 Kepler15b 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 A_{s} of the simulated synthetic TLC. Panel c: variations of the simulated transit duration Δt_{e} − Δt_{s}. Panel d: variations of the simulated transit depth ΔF_{tr}. Panel e: variations of the simulated TLC minimum time Δt_{m}. Panel f: variations of the simulated transit starttime Δt_{s}. Panel g: variations of the simulated transit endtime Δt_{e}. The abscissascale ΔJD = JD − 2 454 833.0 is in Julian days (JD). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.