THE ORBIT AND MASS OF THE THIRD PLANET IN THE KEPLER-56 SYSTEM

Our Keck/HIRES observations were obtained largely following the standard procedures of the California Planet Survey (CPS) team (Howard et al. 2010), modified slightly for the faint stars of the Kepler field, following the approach of Huber et al. (2013). For all observations, we used the C2 decker (140 × 085), which is a factor of four taller than the B5 decker typically used for observations of bright stars. This setup allows for more background light to enter the spectrograph, allowing for better sky subtraction while maintaining a resolving power of R ≈ 50,000.

Each observation was made with an iodine cell mounted along the light path before the entrance to the spectrograph. The iodine spectrum superposed on the stellar spectrum provides a precise, stable wavelength scale and information on the shape of the instrumental profile of each observation (Valenti et al. 1995; Butler et al. 1996).

The integration times range from 600 to 1800 s. The star-times-iodine spectrum was modeled using the Butler et al. (1996) method, with the instrumental profile removed through numerical deconvolution. The RV of the star at each observation is compared to a template spectrum of the star obtained without iodine, with the instrumental profile removed through numerical deconvolution. The observed RVs are listed in Table 1.

The data set used here includes the 10 observations used by Huber et al. (2013), re-analyzed after all observations were recorded. An improved stellar template spectrum causes the measured RV from these observations to be slightly different than those reported by Huber et al. (2013), although the differences are smaller than the formal uncertainties on each observation.

We also obtained 19 observations of Kepler-56 with HARPS-North, a high-precision echelle spectrograph at the 3.6 m Telescopio Nazionale Galileo (TNG) at the Roque de los Muchachos Observatory, La Palma, Spain. HARPS-N is a fiber-fed high-resolution (R = 115,000) spectrograph optimized for measuring precise RVs.

The exposure times for all observations with HARPS-N were 1800 s, and the data were reduced with version 3.7 of the standard HARPS-N pipeline. RVs were derived with the standard weighted cross-correlation function method (Baranne et al. 1996; Pepe et al. 2002). These data are also listed in Table 1.

Note that the HARPS-N pipeline includes the systemic RV, while the Keck/HIRES pipeline does not, leading to a 54.25 km s⁻¹ apparent shift between the two sets.

First, we perform maximum likelihood estimation via Python's scipy.optimize.minimize routine. For the possible companion, we take all values as unknown. Specifically, we fit for $\sqrt{e}\cos \omega$, $\sqrt{e}\sin \omega$, $m\sin i$, P, ${M}_{\star }$, t_tr, and $\dot{\gamma }$, the acceleration of the entire system over time. Since our measurements come from two instruments, we include independent offset terms for each, γ_HARPS and γ_HIRES, where γ is the systemic RV offset term introduced in Section 3. There are 17 free parameters in total—these, plus $\sqrt{e}\cos \omega$, $\sqrt{e}\cos \omega$, and $m\sin i$ for each planet (as mentioned in Section 3).

Maximum likelihood estimation is a process in which we calculate the logarithm of likelihood (L) by comparing our data (D) to the sum of our generated RV curves through the standard equation

$$\begin{eqnarray}&&\mathrm{ln}(L)=-\displaystyle \frac{1}{2}\displaystyle \sum _{i}^{N-1}{\left(\displaystyle \frac{{D}_{i}-{\mathrm{RV}}_{i}}{{\sigma }_{\mathrm{jitter},i}}\right)}^{2}+\displaystyle \frac{1}{2}N\mathrm{ln}(2\pi {\sigma }_{\mathrm{jitter},i}^{2})\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\sigma }_{\mathrm{jitter},i}^{2}={\sigma }_{i}^{2}+{j}^{2}.\end{eqnarray} \tag{ 3 }$$

We use σ_jitter in order to incorporate jitter. Sources of jitter include uncertainties in measurements beyond photon noise that arise from sources like noise in the detector or stellar activity. For sub-giant stars, typical jitter values are 3–5 m s⁻¹ (Johnson 2008). Given the longer exposures for this star relative to previous studies of planets around relatively bright subgiants, we might expect a lower level of jitter as the integrations will average over the higher-order modes.

We initialize the fit with values from Huber et al. (2013). However, we note a typo in Table 1ajaa4237t1 of the discovery paper: the listed times of transit in that paper are too large by 20 days. They should be 2454958.2556 and 2454958.6560 days for Kepler-56 b and c, respectively, rather than 2454978.2556 and 2454978.6560 days.

We reject trials with nonphysical results such as negative masses and periods. For steps that are not rejected, we apply normal priors with expected values and 1σ uncertainties based on measurements from Huber et al. (2013) for the asteroseismic mass of the host star and the inner planets' photodynamical eccentricity vectors, based on the TTV analysis of the Kepler light curve. The sum of the logarithm of each prior term is saved for each set of parameters that is tested.

Then, we calculate the logarithm of the posterior probability for each model, which is the sum of the log-prior and log-likelihood terms (as maximizing the logarithm of a function is equivalent to maximizing the function itself). Equation (4) illustrates this process:

$$\begin{eqnarray}&&\mathrm{ln}[p(\mathrm{RV}| {\rm{D}})]=\mathrm{ln}[\displaystyle \prod _{N}\,({\boldsymbol{\theta }})]+\mathrm{ln}[L(D| \mathrm{RV})].\end{eqnarray} \tag{ 4 }$$

Equation (4) calculates the logarithm of the posterior probability distribution function for any set of model parameters (${\boldsymbol{\theta }}$) as compared to our RV data (D). The combination of parameters found by this process to make the data most probable then becomes the initial guess for our final fitting process.

We use the result of maximum likelihood estimation from Section 3.1 as the initialization for emcee (Foreman-Mackey et al. 2013), a Markov Chain Monte Carlo (MCMC) implementation for Python of the affine-invariant ensemble sampler of Goodman & Weare (2010).

Our 17 parameter simulation uses 150 walkers and 6000 steps, with an observed burn-in of 1500 steps.

A fourth planet beyond the orbit of Kepler-56 d, if it exists, could be observable through the detection of a long-term trend in the data. Given our three-year baseline of observations, we can place limits on the presence of such an outer companion. From our emcee results, we find a long-term RV acceleration of −0.25 ± 0.32 m s⁻¹ yr⁻¹. The 95th percentile value of the emcee posterior probability distribution for $\dot{\gamma }$ provides an upper limit on acceleration from a fourth planet of 0.80 m s⁻¹ yr⁻¹.

From Montet et al. (2014), we know the maximum trend caused by a planetary companion on a circular orbit is

$$\begin{eqnarray}&&\dot{\gamma }=(6.57\,{\rm{m}}\,{{\rm{s}}}^{-1}\,{\mathrm{yr}}^{-1})\left(\displaystyle \frac{{m}_{p}\sin i}{{M}_{\mathrm{Jup}}}\right){\left(\displaystyle \frac{a}{5\mathrm{au}}\right)}^{-2},\end{eqnarray} \tag{ 5 }$$

where m_p is the mass of the planet, M_Jup the mass of Jupiter, and a the orbital semimajor axis. From this, we can place limits on the presence of outer companions with $m\sin i$ larger than 0.49M_Jup at 10 au and 1.95M_Jup at 20 au; such companions must be at particular points in their orbits or at low inclination in order to evade RV detection.

At V ≃ 13 mag, Kepler-56 falls just within Gaia's bright-star limit (Perryman et al. 2014). A fourth planet's acceleration on Kepler-56 in Gaia astrometry might be detectable at the level of 10–20 μas yr⁻² over the course of the mission. Averaging over flat priors for orbital angles and eccentricity, at the nominal distance of Kepler-56 (d ∼ 850 pc), Gaia could in principle detect curvature due to orbital motion of a companion of ≳20 M_Jup at 10 au or ≳80 M_Jup at 20 au. These values in Equation (5) return, at the lowest, an acceleration of 32.85 m s⁻¹ yr⁻¹. This is much higher than the limits returned by our fit, suggesting that, save for face-on orbits, Gaia will be less helpful than continued RV observation in placing further limits on a fourth planet.

A fourth planet in a near-resonant orbit with Kepler-56 d could masquerade as a single eccentric planet, as described by Anglada-Escudé et al. (2010). However, we find the probability of this scenario to be low. Re-running our emcee fit with the outer planet's eccentricity fixed at 0 leads to decreased likelihoods for the fit as a whole, and we do not detect any long-term structure in the residuals. However, our observations do not have the time resolution necessary to make a definitive assertion on this effect. Complete phase coverage of Kepler-56 d is needed to answer this question.

Our research supports that of Huber et al. (2013) in finding strong evidence for a massive, non-transiting exoplanet in the Kepler-56 system. Now that our observations span a full Kepler-56 d orbit, we can compare our results with the projections from Huber et al. (2013), who predicted that both the planetary obliquity and long-term RV trend could both be broadly explained by a non-transiting companion with a period of 900 days and mass of 3.3 ${M}_{\mathrm{Jup}}$.

Both our minimum mass and period are similar to the representative values listed by Huber et al. (2013). Kepler-56 d's minimum mass could be commensurate with that of a giant planet or a brown dwarf (for inclinations below 30°). This could have implications for the near 2:1 resonance of the inner planets' orbits, as well as for the misalignment of their orbital plane with that of Kepler-56's rotation. Indeed, Li et al. (2014) simulated several scenarios and found a higher probability of the observed misalignment being of a dynamical origin (e.g., Fabrycky & Tremaine 2007) than from migration of the bodies in a tilted protoplanetary disk (e.g., Bate et al. 2010) or through angular momentum transport in the star itself that led to an apparent misalignment, even if the system was originally aligned (Rogers et al. 2012).

While Kepler-56 d is a possible source of dynamical perturbation, Gratia & Fabrycky (2016, submitted) simulate the scattering of two giant outer planets and find scattering between a system of three outer planets is required to excite the two inner planets of the system to inclinations similar to those observed in the data while preserving coplanarity. These additional planets, if real, must be scattered to large orbital separations or ejected entirely to evade detection by our RV observations.

Huber et al. (2013) inferred masses of the system's inner planets by dynamically modeling their transits, ignoring possible perturbations from the third, outer planet. We verify that this is a reasonable assumption by checking two effects that may be significant: a tidal term corresponding to the change in the gravitational potential as Kepler-56 d completes its orbit, and a Roemer delay as the distance to the inner planets and host star varies over the orbit of the outer planet.

Following Equations (25)–(27) of Agol et al. (2005), the tidal effects would cause, over a long time baseline, the transits of an inner planet with mass m₁ and period P₁ to be perturbed with a standard deviation

$$\begin{eqnarray}\sigma & = & \displaystyle \frac{3\beta {e}_{2}}{\sqrt{2}{(1-{e}_{2}^{2})}^{3/2}}\\ & & \times \,{\left[1-\displaystyle \frac{3}{16}{e}_{2}^{2}-\displaystyle \frac{47}{1296}{e}_{2}^{4}-\displaystyle \frac{413}{27648}{e}_{2}^{6}\right]}^{1/2},\end{eqnarray} \tag{ 6 }$$

where e₂ is the eccentricity of the outer planet and

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{m}_{2}}{2\pi ({m}_{0}+{m}_{1})}\displaystyle \frac{{P}_{1}^{2}}{{P}_{2}}.\end{eqnarray} \tag{ 7 }$$

Here, m₂ is the mass of the outer planet with orbital period P₂, and m₀ the mass of the host star.

For the values in Table 2 for our system, we find perturbations in the time of transit on the order of 4 s for Kepler-56 b and 16 s for Kepler-56 c. Given that the precision in the measurement of times of transit of these planets is typically tens of minutes, we do not expect these perturbations to affect, or be noticeable in, the measured times of transit.

The light travel time, or Roemer, delay is the result of changes in the physical distance between the observer and the host star due to the orbit of the outer body. Following Equations (6) and (7) of Rappaport et al. (2013), its magnitude is bounded such that

$$\begin{eqnarray}&&{A}_{R}\leqslant \displaystyle \frac{{G}^{1/3}}{c{(2\pi )}^{2/3}}{P}_{2}^{2/3}\left[\displaystyle \frac{{m}_{2}\sin {i}_{2}}{{({m}_{0}+{m}_{1}+{m}_{2})}^{2/3}}\right],\end{eqnarray} \tag{ 8 }$$

where G is Newton's constant, c the speed of light, and all other terms retain their meaning from the previous equation. Inserting values from Table 2 again, we find the expected light travel time signal not to exceed 5 s, significantly smaller than the observed uncertainties, so we do not expect Kepler-56 d to affect the orbits of the inner planets in any observable way.

From our model, we measure a time of transit for Kepler-56 d of BJD-2,400,000 = 56,450 ± 7 days. Thus, if the planet transits the star, we would expect to detect a single transit in the Kepler data set that is visible by eye, but do not observe one in this window. As we know the posterior distribution of allowed times of transit, we can determine the probability the planet transited during a data gap. In Quarters 6 and 7, there are four data gaps larger than 12 hr in which a transit could reside. Together, these gaps represent 1.7% of the mass of the posterior distribution of the time of central transit. The transit duration allows us to place even tighter constraints. If Kepler-56 d transited with an impact parameter b = 0, the transit would have a duration of 3.1 days. As none of the gaps are longer than 20 hr in duration, we can additionally rule out any transits with b < 0.95. By again integrating over the posterior distribution but accounting for the nonzero transit duration, assuming a flat distribution in impact parameter, we find that only 0.07% of allowed transits fall fully inside a data gap. If Kepler-56 d were to transit, there is a 99.93% proabability it would be observable in the Kepler data. Given this low probability and the a priori small transit probability for a companion on a ∼1000 day period, it is likely this companion is non-transiting. We note that non-transiting does not necessarily imply non-coplanarity with the inner planets, as the transit probability decreases with increasing semimajor axis (Borucki & Summers 1984).

Having measured the minimum mass ($m\sin i$) and orbital semimajor axis (a) of Kepler-56 d, we can consider the possibility that the Gaia astrometric mission would be able to constrain its inclination. For lower (more face-on) inclinations, the planet will have a higher mass and the center of mass of the system will move closer to the planet. Additionally, the astrometric orbit will change shape on the sky, with more face-on inclinations appearing more circular throughout an orbit.

Given the distance to the system (d ∼ 850 pc) and the inferred semimajor axis 2.13 ± 0.07 au, the orbit of Kepler-56 d has a projected semimajor axis on the sky of ∼2.5 mas. From the mass ratio between the planet and star, we then expect an astrometric signal with a semi-amplitude of $11{\sin }^{-1}$ μas. Perryman et al. (2014) determined that Gaia will detect planets with astrometric signatures larger than 68 μas for stars as bright as Kepler-56, meaning this planet would evade detection at all except the lowest inclinations. However, given that the Gaia data can be combined with the prior information about the orbit of Kepler-56 d from RVs, it may be possible that the planet will be detected at slightly lower inclinations. Regardless, the prospects of a robust determination of the outer planet's complete set of orbital parameters from Gaia appear unlikely.