A SYSTEMATIC SEARCH FOR TROJAN PLANETS IN THE KEPLER DATA

A transiting trojan planet will manifest itself in a lightcurve as a set of quasi-periodic transits, with a mean period equal to that of the planet with which it shares a 1:1 resonance. If we define the mid-transit phase of the primary planet as ϕ_p = 0, then the typical trojan transit ephemeris will be at phases near ϕ_t = π/3 or ϕ_l = −π/3, corresponding to the 60^o by which the L5 and L4 Lagrangian points are trailing and leading the primary planet's orbit, respectively. The ephemeris will oscillate around these points due to the epicyclic motion of the trojan. As a result of this quasi-periodic behavior, it is useful for the analysis and for illustrative purposes to structure the lightcurve in a "river diagram" framework, which has been previously used in the analysis of quasi-periodic orbits (Carter et al. 2012; Carter & Agol 2013). A river diagram is a folded lightcurve in matrix form, in which each column represents a specific phase with respect to a fixed orbital period, and each row is one such period. A planet with a constant period will leave a straight vertical trace in a river diagram folded to its period, whereas a quasi-periodic planet will leave a wiggly vertical trace if the diagram is folded to its mean period. A general example of a river diagram is shown in Figure 1. In our case, we create one river diagram per KOI by folding to the KOI period, and impose ϕ_p = 0 in practice by shifting the primary planet mid-transit phase to the first column of the matrix. Hence, in this framework, transiting trojan planets will be represented by wiggly traces near specific columns (corresponding to ϕ_t = π/3 or ϕ_l = −π/3) in the matrices.

Zoom In Zoom Out Reset image size

Before converting the lightcurves into river diagrams, we subject the data to a subtraction of a median moving box filter with a box width equal to four times the transit duration. This eliminates most low-frequency variations, which benefits visual inspection, and does not significantly affect the strength of the transit feature (the specific morphology of the transit is affected, as it becomes slightly elevated with respect to the mean continuum flux, but since the continuum directly before and after the transit is elevated along with it, the strength and visibility of the feature itself is conserved). The lightcurves used were those produced by the standard Kepler pipeline, for quarters Q0–Q12. The KOIs and their values were taken from the official Kepler listing on March of 2013, and include those (and only those) KOIs that were either classified as planets or planet candidates at that time. This gives a total of 2244 KOIs, where the count includes individual KOIs in candidate multi-planet systems.

As mentioned in Section 2.1.1, a typical trojan planet will not merely stay fixed at its Lagrangian point, but will oscillate around it in a manner that can be characterized by two superpositioned epicycles, Ω₁ and Ω₂. As derived in Laughlin & Chambers (2002), these cycles have periods that scale linearly with the orbital period P_orb, and that otherwise depend only on the masses, in the following ways:

$$\begin{equation} P_1 \sim P_{\rm orb} \left(1 + \frac{27(m_1 + m_2)}{8(m_0 + m_1 + m_2)} \right) \end{equation} \tag{ 1 }$$

$$\begin{equation} P_2 \sim P_{\rm orb} \sqrt{\frac{4(m_0 + m_1 + m_2)}{27(m_1 + m_2)}}, \end{equation} \tag{ 2 }$$

where m₀ is the mass of the star, m₁ the mass of the primary planet, and m₂ the mass of the trojan. Since m₀ ≫ m₁ + m₂, it follows that P₁ is close to P_orb, and that P₂ is substantially longer than P_orb, typically by a factor of ∼10–100 for realistic and potentially detectable planetary systems. The amplitude of Ω₁ is significantly smaller than that of Ω₂ (by a factor of ∼25, as we will see later).

Due to these oscillations, the transit ephemeris of the trojan with respect to the primary planet will vary with time, as the trojan is alternately located ahead of and behind the Lagrangian point around which it librates. In order to be able to systematically detect trojans, we must thus simulate the ephemerides for the full range of epicyclic orbits that the trojan may exhibit. We do this numerically, using the acceleration equation for the full three-body problem adapted from Laughlin & Chambers (2002):

$$\begin{equation} \ddot{\boldsymbol {r_2}} = -\frac{G(m_0 + m_2)}{r_2^3} \boldsymbol {r_2} - \frac{G m_1}{\Delta ^3} \boldsymbol {\Delta } - \frac{G m_1}{r_1^3} \boldsymbol {r_1}, \end{equation} \tag{ 3 }$$

where $\boldsymbol {r_1}$ is the vector from m₀ to m₁, $\boldsymbol {r_2}$ is the vector from m₀ to m₂, and $\boldsymbol {\Delta }$ is the vector from m₁ to m₂. The baseline simulation setup has a time step of 10⁻⁵ yr and contains 2 × 10⁶ steps, in a system where m₀ = 1 M_Sun, m₁ = 1 M_Jup, and m₂ = 1 M_Earth, and with a semi-major axis of 1 AU. A few examples of a trojan orbit are shown in Figure 2, to illustrate the epicyclic motion. We then extract the phase ϕ(t) of the trojan in a frame that co-rotates with the primary planet. This is done for a series of simulations with a different initial perturbation δ, where δ is defined in such a way that if δ = 0.001 AU, then the starting point of the trojan deviates by 0.001 AU from the relevant Lagrangian point in both the x- and y-direction. We vary δ between 0 and 0.004 AU, where the latter corresponds to a total epicyclic amplitude of 23^o, roughly equal to the range of Jupiter's trojans. The resulting ϕ(t) is fit with a series of sine functions. For small initial perturbations, both Ω₁ and Ω₂ are well represented by pure sine curves, but as the amplitude increases, Ω₁ becomes increasingly asymmetric between its peaks and troughs, as the orbit starts acquiring the asymmetrically elongated characteristic that has given it the name "tadpole orbit." We find that this characteristic can be well approximated by adding a co-phased sine wave with Ω₂, with twice the frequency, in addition to a constant offset A₄. That is,

$$\begin{eqnarray} \phi (t) &=& A_1 \sin {(\omega _1 t - B_1)} + A_2 \sin {(\omega _2 t - B_2)}\nonumber\\ &&+ A_3 \sin {(2 \omega _2 t - B_2)} + A_4. \end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Here, B₁ and B₂ are arbitrary phase shifts, and each A_i is a function of δ. We determine these functions by fitting the above relation to the series of simulations with different δ, and fit power laws to the results, since they are well represented by either linear or quadratic fits. We derive the following dependencies:

$$\begin{equation} A_1 = 2.8322 \delta \end{equation} \tag{ 5 }$$

$$\begin{equation} A_2 = 50.3711 \delta \end{equation} \tag{ 6 }$$

$$\begin{equation} A_3 = -525.1390 \delta ^2 - 0.0467 \delta \end{equation} \tag{ 7 }$$

$$\begin{equation} A_4 = 1691.7934 \delta ^2 - 0.1014 \delta + 1.0472. \end{equation} \tag{ 8 }$$

As can be seen in Figure 3, the fits are very good but not perfect, and the standard deviation of the residuals to the fit can be characterized by

$$\begin{equation} e_{\phi } = 273.0478 \delta ^2 - 0.3034 \delta. \end{equation} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

Due to the quadratic dependence, the error is only significant for the largest values of δ. We will discuss the global results of errors on the ephemerides in Section 2.2, but as an illustrative example, if δ = 0.004 AU, the error corresponds to 1 Kepler long cadence (t_step = 0.5 hr) if the orbital period is ∼42 days (∼2000 t_step).

The search for signals from potential trojan planets in the data was conducted among all the planet/candidate KOIs using the two complementary approaches of a visual inspection, and a systematic algorithmic search. The visual inspection allows to efficiently discern any false positives and potentially detect unexpected signals or types of orbits that are not actively searched for in our automatic algorithm, such as horseshoe orbits. The algorithmic search, on the other hand, provides the more objective and quantified answer as to whether signals exist in the data or not. All searches were based on the river diagrams explained in Section 2.1.1. Our visual inspection did not reveal any convincing trojan candidates, but a few interesting cases of other phenomena are discussed in Section 2.3.

The idea behind the automatic procedure is the following: for each KOI, we wish to be able to detect any possible corresponding trojan orbit, so long as it is co-inclined with the KOI and has an amplitude of epicyclic motion within that of the trojan swarms of Jupiter (in addition, of course, the trojan has to be of sufficiently large size to be detectable at all in the data). We therefore wish to search a grid of all possible orbits within these constraints, where for a given orbit, a trace is predicted in the river diagram of the KOI, and the would-be transit signal along this trace is summed up. Any orbit among the family of orbits tested that shows a significant feature can then be examined in detail to assess whether or not it is real and to better characterize its properties.

In practice, the family size of orbits for a given KOI varies in extent, and becomes quite substantial for certain types of objects. There are four parameters that need to be scanned through. The first parameter is the perturbation amplitude, δ. This sets all of the A_i factors in Equation (4). We set the range of this parameter such that it covers the oscillation amplitude range of the Jupiter trojans, and we set the grid steps so that the peak-to-peak amplitude is sampled once for every time it increases by 1 t_step. This is one of the steps to ensure that the positional error never gets larger than 1 t_step for any single transit event.

The second parameter is a mass parameter, which is essentially determined by m₁. The two aggregate masses m₀ + m₁ + m₂ ≈ m₀ and m₁ + m₂ together set the two oscillation periods P₁ and P₂, as discussed at the top of this section. While the mass of the primary m₀ has some considerable uncertainty, it is dwarfed by the uncertainty in m₁, since we only know the radius of the (candidate) planet in general, and it can have a wide range of densities. We have furthermore assumed that m₂ < m₁, and so m₁ provides the most important parameter uncertainty. We therefore fix m₀ to the value given in the Kepler catalog, set m₂ to zero for this purpose, and examine a wide possible range of values for m₁. We estimate this range individually for each KOI, based on the estimated planetary radius from Kepler, and the mass–radius relationship of those transiting planets for which both quantities have been measured (Weiss et al. 2013). In doing so, we note that while the scatter in the distribution is large, it is not infinite. All objects with a certain radius span a mass range where the lowest and highest mass are within a factor of 10 of each other, apart from in the super-Jovian mass regime where the mass of a ∼1 R_Jup object can in principle even be in the brown dwarf range, but it is known that such objects are very rare (e.g., Grether & Lineweaver 2006). Hence, we adopt the mean relationships in Weiss et al. (2013) and establish that for objects with r_p < 1 R_Jup, the lower mass limit is $m_{\rm low} = 0.178 r_{\rm p}^{1.89}$ and the upper mass limit is $m_{\rm high} = 1.78 r_{\rm p}^{1.89}$. For objects with r_p > 1 R_Jup, we set a lower limit of 0.3 M_Jup and an upper limit of 3 M_Jup. In other words, we try to generously encompass most properties that the planet could have, while simultaneously excluding some parameter space where nothing could realistically reside. With a mass range in hand, we translate it to a distribution in the pair of periods P₁ and P₂. Within the resulting range, we step through the periods in units of P_orb, again to ensure that the positional error is smaller than 1 t_step.

Finally, the third and fourth parameters are simply the phase shifts B₁ and B₂, which determine which phase each oscillation is in at a given time. We step through B₂ in steps of P_orb up to P₂. B₁ is insignificant when A₁P_orb/2π is smaller than t_step which is the case for a rather large range of values of δ, hence its sampling is set to depend on the amplitude in such a way that 4N_samp uniformly spaced values of B₁ are stepped through, where N_samp = A₁P_orb/2πt_step, rounded downward to the nearest integer.

Looping through all of these parameter ranges for each set of parameters we can generate one specific orbit, which yields a predicted trace in the river diagram, where the central location of each transit is rounded off to the nearest t_step. Summing up the putative transit signals along this trace, we can make a first assessment of whether or not a trojan planets exists in this particular orbit. In practice, in order to save substantial computational effort, we do this by creating signal-to-noise (S/N) versions of each river diagram prior to running the automatic trojan search. In the S/N diagram, each pixel represents the S/N of an individual hypothetical transit centered on each given t_step. That is, for each pixel in the original river diagram, a sum is taken of the signal from t_step − t_dur/2 to t_step + t_dur/2. This sum is divided by the standard deviation of the signal in the local continuum to produce a S/N value at that t_step. The search algorithm then creates a sum along the trace in the S/N diagram. The resulting value, divided by $\sqrt{N_{\rm per}}$, where N_per is the number of rows in the river diagram, gives an estimate of the ensemble S/N for a given trojan orbit. This S/N calculation provides a more robust estimation of the real S/N than a standard root-mean-square estimation, because a strong source of contamination for a given KOI is the presence of transits in the data from other KOIs in multi-planet systems. Such contaminants are (for practical purposes) arbitrarily distributed in phase in a river diagram, and so occasionally one or a few contaminant transits will overlap with the trace of a trojan orbit under investigation. This contamination has a much larger impact on a root-mean-square estimation than our estimation, which favors detection of uniformly distributed events across the trace over strong single-instance events. Nonetheless, our mean-based estimator is still somewhat susceptible to such contaminants, hence we have also used a median-based estimator for the S/N, which is calculated as the median of the trace multiplied by $\sqrt{N_{\rm per}}$.