Light-Curve Evolution due to Secular Dynamics
and the Vanishing Transits of KOI 120.01

The Kepler mission (Borucki et al., 2010) has led to a revolution in our knowledge of planetary systems (e.g. Batalha et al., 2013)), providing information on the orbital architectures and the multiplicity distribution of planetary systems. In particular, it was found that a single multiplicity distribution is insufficient to explain the excess of singly transiting planets systems (Ballard & Johnson, 2016; Lissauer et al., 2011), and a second population, consisting of either single, or more inclined planets, may be required to explain the observations (also known as ”the Kepler dichotomy”). A recent statistical study of planetary population found it unlikely that this dichotomy is explained solely by a population of real single-planet systems (He et al., 2019). Statistical analysis of the distribution of transit durations has shown that multiply-transiting systems are consistent with having a roughly co-planar structure (Fabrycky et al., 2014). With accurate determination of stellar densities Xie et al. (2016) used transit duration distribution to fit for distribution of eccentricities (and mutual inclinations, for multiple-transiting systems), and inferred that there is a dichotomy within the singly-transiting planets population, which consists of dynamically hotter and dynamically colder systems.

Non-transiting companions may hide in systems with singly transiting planets particularly in mutually inclined systems (Lai & Pu, 2017), and the former may be probed by their influence on the latter (Mills et al., 2019). Secular effects applied by an exterior companion on an inner transiting planet would yield transit variations, usually described as TTV’s (Agol et al., 2005; Holman & Murray, 2005) and Transit Duration Variation (TDV’s) (Miralda-Escudé, 2002). Secular precession rates due to a planetary perturber were estimated for many systems (e.g., Moro-Martín et al., 2007; Lee & Peale, 2003). The time scale of the precession period of a planet due to a secular external perturber is $\tau=4P_{1}/(\mu_{2}\alpha^{2}b_{3/2}^{(1)}(\alpha))$, where $P$ is the orbital period, $\mu_{2}$ is the perturber/host mass ratio, $\alpha=a_{1}/a_{2}$, where $a$ is the semi-major axes of objects, and 1 and 2 correspond to the inner and outer planets, respectively, and $b_{3/2}^{(1)}$ is a Laplace coefficient. Derivations can be found in Murray & Dermott (1999, chapter 7), and expressions for nodal precession rates to second order in inclination were introduced by Xu & Fabrycky (2019). For comparison, Earth’s orbital precession period due to only Jupiter is of the order $10^{5}$ years. For a close-in planet with an 8 day orbit perturbed by an external Jupiter-mass planet with a 35-day orbit, the precession period is of order 500 years. During Kepler’s primary mission time of $\sim$4 years, this corresponds to an accumulated precession angle of $\sim$3^∘, which may be observable in the current data.

Other sources for secular precession are stellar quadruple moment (e.g., KOI 13.01, a hot Jupiter which exhibits nodal precession of $\sim$0.7 ^∘$\rm yr^{-1}$ (Szabó et al., 2012)), general relativity (e.g, HD 80606 b (Blanchet et al., 2019), and Mercury (Einstein, 1916)), and at specific configurations, even tidal interactions (Wu & Goldreich, 2002). In these examples, different drivers yield a similar dynamical phenomenon of secular precession.

In this work, an analytic method for photo-dynamical fitting of the light-curve is presented. The light-curve calculation relies on the assumption that the orbital motion is dominated by secular precession, a likely condition for systems in which there are no near-resonant interactions. As described above, multiple physical processes result in similar behaviour (secular precession) and depending on the physical process the result of such a fit may inform us on the system’s orbital elements (e.g. planet-planet interaction) or other properties (e.g. stellar $J_{2}$). There are a few advantages for this method. First, it does not require a massive number of n-body integration iterations. Second, it utilizes only the transiting planet information to construct a model for the orbital motion. Hence, the mass and orbital elements of the perturbing object are not required (which are often unknown since long-period planets rarely transit). Only solutions for the orbital elements of the transiting, perturbed, planet are sought. Third, it establishes a direct link between the parameters of the model and the dynamics of the system with respect to its invariable plane.

Searching for transiting candidates which could be affected by slow, secular, transit variations, we came across the interesting case of KOI 120.01. This object exhibits unusual vanishing transiting signals, which can be explained by this model.

In §2 the mathematical description for transit variations in the presence of secular precession is presented. In §3 the method is applied to the vanishing transit signal of KOI 120.01, and a few possible models for explaining the data are presented. The results are summarized and discussed in §4.

The aim of the following mathematical derivation is to relate the instantaneous orbital elements with respect to the system’s invariable plane (inclination $I$, longitude of the ascending node $\Omega$, and argument of periapse $\omega$) to their counterparts measured with respect to the sky plane ($I_{\rm sky},\omega_{\rm sky}$). The longitude of ascending node with respect to the sky plane ($\Omega_{\rm sky}$) is irrelevant for the transit shape, as long as the stellar disk is radially symmetric. The coordinates of the planet parameterized by the true anomaly $f$, are given by (Murray & Dermott, 1999, page 51):

where $x,y,z$ are the coordinates of the planet with the convention that $x,y$ is the system’s invariable plane, $x$ lies in the sky plane, and $r$ is the distance from the central mass. The equation describing the sky plane is given by $y\cos{\beta}+z\sin{\beta}=0$, where $\beta$ is the angle between the observer’s line of sight and its projection on the invariable plane ($\beta=0$ corresponds to viewing the invariable plane edge-on, such that the sky plane is simply $x-z$). Solving for the intersection of the orbit with the sky plane yields two solutions for the true anomaly $f$ (representing the ascending and descending nodes through the sky plane). This allows for expressing the orbital elements with respect to the sky plane in terms of those with respect to the invariable plane. The sky plane elements are not required for constructing the light-curve, but give a geometrical interpretation revealing the changes in impact parameter and mid-transit velocity. The sky plane elements are

Gradual changes of the angles $\Omega$ and $\omega$ due to nodal regression or apsidal advance will affect the transit shape in a several ways. The varying separation $s$ between the planet and star will affect the impact parameter and hence the depth and duration of the transit and the shape of ingress and egress. The different longitude observed at each transit will yield a different mid-transit tangential velocity $v_{\rm mid}$, changing in turn the transit duration. Lastly, the different time spent by the planet at each longitude will result in Transit Timing Variations (TTV’s).

Making use of the relation $f_{\rm mid}=\pi/2-\omega_{\rm sky}$ (with $f_{\rm mid}$ the true anomaly at mid-transit) the transit parameters can be expressed as

and

where $n=2\pi/P$ is the mean motion, $a$ the semi-major axis and $e$ is the orbital eccentricity. The orbital elements $a,e,\Omega,\omega,I,M$ (where $M=n(t-t_{p})$ is the mean anomaly and $t_{p}$ is the time of periapse) are translated instantaneously to invariable plane coordinates ($x,y,z$) through solution of Kepler’s equation. Using the angle $\beta$, they are then translated to sky plane coordinates ($X,Y$). From these coordinates and the ratio between transiting planet and host star radii $R_{p}/R_{*}$, relative flux is computed, and a model light-curve is constructed (Mandel & Agol, 2002). Instead of using $t_{p}$, here the parametrization is for the first time of mid-transit $T_{\rm mid}$, which is known to a good accuracy.

To briefly summarize: given a full set of orbital parameters (and $R_{p}/R_{*}$) a model flux point can be directly calculated. Importantly, these orbital elements need not be constant in time, and for a given time evolution of these elements (i.e. a list of such sets) a full light curve can be calculated. Below we assume a secular model in which $\Omega,\omega$ are linear in time. Specifically, by assuming all the dynamics is captured by some long term variability in $\Omega,\omega$ a full light-curve can be generated through the parameters $P,T_{\rm mid},R_{p}/R_{*},a/R_{*},e,\omega_{0},\dot{\omega},I,\Omega_{0},\dot{\Omega},\beta$ (where the subscripts ”0” refers to the value at $T_{\rm mid}$) without a full N-body integration.

Fitting for the parameters listed above (or any combination of them) can reveal the precession rate of the transiting planet, imprinted in small variations in the photometry of the transits. As the ratio $\dot{\omega}/\dot{\Omega}\sim-2$ (or, equivalently, $\dot{\varpi}/\dot{\Omega}\sim-1$, where $\varpi=\Omega+\omega$ is the longitude of periapse), instead of fitting for $\dot{\omega}$, the fit is performed on the ratio $\dot{\omega}/\dot{\Omega}$. The exact ratio between the secular precession rates depends on the Laplace coefficients, and derivations can be found in Murray & Dermott (1999, chapter 7). Instead of using the angle $\beta$ which is unconstrained, the angle $I_{\rm sky0}$ was used (i.e., the inclination with respect to the sky plane at first mid-transit time). This angle is constrained to be up to a few degrees away from 90^∘ by the existence of observable transits.

In the first part of this work a photodynamical analytic approach is presented. Such an approach is useful when the temporal functional dependence of the orbital elements is known and only the magnitude of these effects is sought. Specifically, the model described here is relevant to a scenario in which a distant object gravitationally perturbs an inner transiting planet in a secular manner, and thus the orbital motion is dominated by precession of the periapse and regression of the nodes over the observed epoch. In the second part of the work, this analytic approach is utilized to fit for the vanishing transit signals of KOI 120.01, assumed to arise due to such a configuration.

It was found that a planet candidate exhibiting both large apsidal precession and large nodal regression rates can account for the observed transit signal of KOI 120.01, though the phase space is large and multiple solutions are possible given the current data. Four scenarios were treated separately: in each the transit is diluted by a different third light source, corresponding to possible arrangements of the relevant objects within the Kepler pixel. Though various sets of parameters were found to fit the data with similar quality, the obtained precession rate magnitude was always in no less than about $1^{\circ}{\rm yr}^{-1}$, and some solutions were found with precession rates as high as $\sim 20^{\circ}{\rm yr}^{-1}$. If indeed the precession scenario describes this transit signal, then the precession rates found here are unusual in their magnitude. For example, the rapidly precessing orbit of KOI 13.01 exhibits a nodal regression rate of about 0.7^∘/year (Szabó et al., 2012).

We examine the physical inferences that can be drawn for KOI 120.01 from a scenario of a rapid precession induced by an external unseen companion. Importantly, we note that even when good knowledge of the precession rate is available, the physical properties of the perturber cannot be uniquely constrained, as the observed rate depends on the perturber semi-major axis and mass, and on the mutual inclination between the orbits, in a degenerate manner. The following analysis, therefore, aims to rule out scenarios inconsistent with the data, rather than finding a specific solution.

In order to account for such a precession rate by a low-inclination planetary perturber, a few Jupiter-masses object is required to be located at a semi-major axis such that $\alpha$ is approximately 0.6-0.75, close either to the 2:1 MMR or to the 3:2 MMR. Because the TTV shape is parabolic over the $\sim$780 days during which transits are observed, the super-period of such an interaction must exceed $\sim$ 1500 days. This corresponds to a distance from resonance $\Delta\lesssim 0.01$. Such a close proximity to resonance of two giant planets is not common, but possible, for example in the system Kepler-9 (Holman et al., 2010; Freudenthal et al., 2018). However, a close, massive planetary perturber would also induce a chopping signal of order $(\mu^{\prime}P/2\pi)(P/(P-P^{\prime}))^{2}$ (Hadden & Lithwick, 2016), where $\mu^{\prime}$ and $P^{\prime}$ are the planet-to-star mass ratio and orbital period of the perturber, respectively. In this case, the chopping signal would be of order of ten minutes, which is not seen in the transiting object’s TTVs (Figure 3). Therefore, a near-MMR interaction with a giant planet is ruled out by the data.

Another possibility is that the motion of the planet out of transit does not result from secular interactions, but from other dynamical phenomena which are not covered by our secular model. For example, the planets orbiting the star K2-146 exhibit rapid precession combined with an orbital mean motion resonance, that together yield rapid impact parameter variations (Hamann et al., 2019; Lam et al., 2019). Analysis of such scenarios requires adding yet more objects to the model, and is beyond the scope of this study and is left for future work. We focused on scenarios with the minimal number of objects.

In order to explain such a precession rate by a distant perturber located at a few AU (which would yield secular interactions, far from any MMR), its mass would need to be comparable to the host star, i.e. it would be a binary member. This scenario, if correct, would be notable because of the difficulty in forming giant planets between close binary pair of Sun-like stars. Indeed, planets in S-type orbits were previously found only in binaries separated by more than $\sim 5$ AU (Gong & Ji, 2018). For an S-type planet in a binary system, the question of dynamical stability arises. Holman & Wiegert (1999) estimated $a_{c}$, the critical semi-major axis below which a planet can stably orbit one of the binary members. For a circular orbit of the binary and equal-mass stars, their expression reduces to $a_{c}=(0.274\pm 0.008)a_{b}$, where $a_{b}$ is the binary semi-major axis. In our case, if the host star is Sun-like, the orbital period of the presumed planet is translated to a semi-major axis of $\sim 0.147$ AU, so for a binary star located at a few AU such a planet is within the stable region. A recent study examined a larger number of orbital configurations for S-type orbits of a planet in a binary system, and reached a similar $a_{c}=0.26a_{b}$ for equal-mass binary members at a circular orbit (Quarles et al., 2019, fig. 1).

In order to check the feasibility of a configuration of an S-type planet, we ran a several N-body integrations of a scenario in which the perturber is a binary of the same mass as the host star, located at a few AU. We find that the precession rate matches the theoretical secular values well, but the TTV’s might be larger than the observed ones. We note that the secular model used above takes into account only the TTV part that arises from the orbital precession, but not other dynamical effects which become more and more significant at higher eccentricity values of the transiting planet. For this reason, the secular model underestimates the TTV, and is prone to overestimate the eccentricity, or overestimate the precession rate required to explain the data. This may explain the high eccentricity and precession rates values of the presumed transiting planet obtained by using the secular precession model. The explanation above is illustrated in Fig. 5, which shows sample TTV patterns and nodal regression rates for a few specific scenarios. Specifically, it shows that the scenario of secular precession combined with a parabolic TTV patterns can occur for an S-type planet perturbed by the binary member, and hence qualitatively, the configuration of an S-type planet perturbed by a binary member is possible for KOI 120.01, if the binary separation is large enough.

We conclude that the simple precession model utilized here yields a decent fit to the photometry (§3.3), but the physical inferences derived from the obtained solutions, along with experimenting with N-body simulations, imply that the dynamical nature of the system is likely more complex. Full N-body models (and probably more observational data) would be required in order to constrain the characteristics of the perturber and the transiting object.

The transit signal may still arise from an eclipsing binary, itself perturbed by a tertiary massive star. If this is the case, the orbit should be inclined and eccentric as no secondary eclipse is visible.

A possible clue that can guide future observations is available from the detection of the motion of the center of light on the Kepler CCD. For two objects of a comparable apparent magnitude, and one partially occulted with a relative transit depth of $\sim 2\times 10^{-3}$, the expected motion of the center of light due to the transit is $\sim 1\times 10^{-3}$ of their separation. Such a motion is indeed reported in the Kepler database³³3https://exoplanetarchive.ipac.caltech.edu/data/ KeplerData/011/011869/011869052/dv/kplr011869052-20160209194854_dvr.pdf. According to the Kepler centroid analysis, examining the in-transit points, the centroid RA increases by $\sim$2 mas and the DEC decreases by $\sim$1 mas. Gaia data places object 2135199457417898240 at a separation from 2135199461713461120 of -1.884 as in RA and 0.9978 as in DEC, matching the Kepler observed centroid shift if the transits are of 2135199457417898240. To emphasize, both the magnitude and direction of the centroid shift are in agreement with predictions based on Gaia astrometry, lending significant support to Gaia 2135199457417898240 as the transits host, and to estimating that the perturber itself does not contribute significant amount of light.

Followup RV measurements dedicated to finding this deduced companion star and the observed transiting object would need to be able to resolve two objects separated by 1.56 as. As the orbital period of the presumed planet is about 20.5 days, the RV semi-amplitude arising from it should be of order a few to a few tens of m/s, achievable with the observational capabilities existing today. The RV semi-amplitude arising from the binary member would be much larger. Follow-up observations, including RV measurements, are encouraged in order to shed more light on the nature of this unique Kepler object, and determine if indeed it harbours an unusual planetary system.